方腔驱动流源代码(SIMPLE算法)
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
通量 F
2
1 M 1 M 1 M 1 1 M 1 u i 1, j u iM u iM , j u i 1, j u i , j 1, j u i , j 4 Re x 1 1 M M 1 M 1 1 M 1 viM u iM 1, j vi , j u i 1, j u i , j 1, j u i , j 4 Re y
F 1 u 2 F 2 v 2
通量 F
1
1 u , Re x 1 v , Re x
G 1 uv G 2 uv
1 u Re y 1 v Re y
, G 1 分别定义在主控制体的中心和角点，如图所示，并按照如下方法离散
F 1 G 1
M 1
p ，代入修正公式得到 M+1 时刻的速
, v M 1 和压力场 p M 1 。将 M+1 时刻的速度场 u M 1 , v M 1 和压力场 p M 1 作为新的
猜测的速度场和猜测的压力场估计值， 采用上述方法计算下一个时刻的速度场和压力场， 直 到满足收敛条件。 收敛判据
Max bip ,j
1 viM vi*, j ,j
将修正后的速度分量代入离散后的连续性方程，得到压力修正方程
p p aip , j pi , j a p ,q p p ,q bi , j
其中
a a
p i 1, j
y 2 y 2 x 2 x 2 p p p u , ai 1, j u , ai , j 1 v , ai , j 1 v ai , j ai 1, j ai , j ai , j 1
xy M 1 M 2 2 2 M 1 M 1 vi , j vi , j Fi ,2 x 0 j Fi 1, j x Gi , j Gi , j 1 y pi , j 1 pi , j t
其中，数值通量
p
pi, j
x 0
编号为（i,j）的速度修正量 u , v 不仅与压力修正量
p 有关，还与邻近点的速度源自文库正
量有关。SIMPLE 算法的重要假定：速度的改变只与压力的改变有关，忽略邻近点对速度修 正的影响。因而得到如下速度修正量
u i, j vi, j
, G 1 和 F 2 , G 2 的某些项冻结于 M 时间层，使离散化之后的方程对
u M 1 , v M 1 是线性的。将离散化之后的 F 1 , G 1 和 F 2 , G 2 代入离散后的 x 方向和 y 方
向的动量方程，整理之后得离散后的动量方程如下
1 M 1 u M 1 M 1 aiu, j uiM au y 0 ,j p , q u p , q bi , j pi 1, j pi , j 1 M 1 v aiv, j viM av ,j p , q v p , q bi , j M 1 i , j 1 1 piM ,j
以上是 SIMPLE 算法中离散化的动量方程
三、SIMPLE 算法基本思想
SIMPLE 算法是一种解决压力-速度耦合问题的“半隐式”算法。首先给定 M 时刻猜测 的速度场 u 估计值 值u
* M
, v M ，用于计算离散动量方程中的系数和常数项。给定 M+1 时刻猜测的压力场
p * ，迭代求解离散动量方程，得到 M+1 时刻速度场的估计值 u * , v * ，速度场的估计
y pi1, j pi, j aiu, j x pi, j 1 pi, j aiv, j
修正后的速度分量
1 u iM u i*, j ,j
y pi1, j pi, j aiu, j x pi, j 1 pi, j aiv, j
1 2 2 xy 1 M M M M aiu, j y uiM x viM 1, j u i 1, j , j vi 1, j vi , j 1 vi 1, j 1 Re x Re y t 4 4
图 2 交错网格示意图
图 3 求解域离散示意图 图 3 中阴影部分代表方腔内的流动区域，阴影区域的边界代表方腔的上、下、左、右 壁面， 阴影区域外面的网格节点是为边界处理需要而设定的虚拟网格节点， 后面介绍边界处 理方法时详细论述。
方程离散
无量纲化的守恒型不可压缩 N S 方程为
U 0 U 1 2 UU P U 0 t Re
本算例采用求解不可压缩流动的经典算法， 即 SIMPLE 算法， 求解方腔内粘性不可压缩 流体运动的定常解。 SIMPLE 算法的全称为 Semi-Implicit Method for Pressure-Linked Equations， 即求解压力关联方程的半隐式算法。 采用 SIMPLE 算法时， 为了避免中心差分格式将 “棋盘” 型参量分布误认为是均匀分布， 需要用交错网格对计算域进行离散。
交错网格
交错网格如图 2 所示，压力、密度等物理量存储在控制体 i, j 的中心，这个控制体称 为主控制体。速度分量
u, v 分别存储在主控制体的 i 1/ 2, j 和 i, j 1/ 2
位置处，
标记为 i, j 位置，再分别以此为中心，划分速度分量 u、v 的控制体。采用空间均匀网格， 等间距离散整个求解域，如图 3 所示。
, G 2 分别定义在主控制体的中心和角点，如图所示，并按照如下方式离散
F 2 G 2
通量 F
1
1 M 1 M 1 M 1 1 M 1 vi 1, j viM viM , j vi 1, j vi , j 1, j vi , j 4 Re y 1 1 M M 1 M 1 1 M 1 u iM viM 1, j u i , j vi 1, j vi , j 1, j vi , j 4 Re x
biv, j
xy M vi , j t
1 1 1 1 M M aiv, j 1 x viM , aiv, j 1 x viM , j vi , j 1 , j 1 vi , j Re y Re y 4 4 1 1 1 1 M M aiv1, j y u iM , aiv1, j y u iM , j 1 u i , j 1, j 1 u i 1, j Re x Re x 4 4
其中， u , v 和
p 分别速度和压力的修正量，修正量亦满足离散的动量方程
u aiu, j ui, j a u p ,q u p ,q bi , j pi 1, j pi , j y 0 v aiv, j vi, j a v p ,q v p ,q bi , j i , j 1
SIMPLE 算法求解方腔内粘性不可压流动
一、问题描述
假设 0 x, y 1 的方腔内充满粘性不可压缩流体，左、右、下壁固定，上壁以
u 16 x 2 1 x 2 运动，试求 Re 100, 200, 400 时的定常解，方腔如图 1 所示。
图 1 方腔内流动示意图
二、离散格式
p i, j
y 2 y 2 x 2 x 2 u u v v ai , j ai 1, j ai , j ai , j 1
* * * * bip , j y u i , j u i 1, j x vi , j vi , j 1
采用迭代法求解压力修正方程，得到压力修正量 度场 u
, v * 满足如下离散方程。
* u * * aiu, j ui*, j a u p ,q u p ,q bi , j pi 1, j pi , j y 0 * v aiv, j vi*, j a v p ,q v p ,q bi , j * i , j 1
p
pi*, j
x 0
一般地， 速度场 u
*
* * * , v * 不满足离散的连续性方程， 因而需要对速度场 u , v 和压力场 p
进行修正。M+1 时刻的修正值和估计值有如下关系
u M 1 u * u v M 1 v * v p M 1 p * p
边界条件处理
首先对计算区域离散，并流动边界之外扩充一个虚拟网格，将真实流动的离散域包围， 如图 7 所示。
图 7 虚拟网格扩充示意图
为很小的正实数，视计算的精度要求而定。本算例中取 1e 8 。
若 bi , j 0 ，则 pi, j 0 ，此时 U i , j
p M 1
U i*, j ，从而来自于离散动量的 U i*, j 满足离散的连
续性方程。因此 Max
b
p i, j
可以作为收敛判据。
xy M 1 M 1 1 1 M 1 M 1 ui , j ui , j Fi ,1j Fi y 0 1, j y Gi , j Gi , j 1 x pi 1, j pi , j t
Y 方向动量方程在速度 v 控制体上离散，时间采用前差
1 2 2 xy 1 M M M M aiv, j x viM y uiM , j 1 u i , j 1 , j 1 u i , j u i 1, j 1 u i 1, j Re y Re x t 4 4
其积分形式为
n UdS 0
S
V
V
u 1 dV n UudS pn x dS n udS 0 S S t Re S v 1 n vdS 0 dV n UvdS pn y dS S S t Re S
p
x 0
其中
biu, j
xy M ui, j t
1 u 1 1 1 M M aiu1, j y u iM , ai -1, j y u iM 1, j u i , j , j u i -1, j Re x Re x 4 4 1 1 1 u 1 M M aiu, j 1 x viM , ai , j 1 x viM , j vi 1, j , j 1 vi 1, j 1 Re y Re y 4 4
图 4 主控制体
图 5 速度 u 控制体
图 6 速度 v 控制体
采用有限体积法离散 N S 方程，连续性方程在主控制体上离散
u
M 1 i, j
1 M 1 1 uiM viM 1, j y vi , j , j 1 x 0
X 方向动量方程在速度 u 控制体上离散，时间采用前差
2
1 M 1 M 1 M 1 1 M 1 u i 1, j u iM u iM , j u i 1, j u i , j 1, j u i , j 4 Re x 1 1 M M 1 M 1 1 M 1 viM u iM 1, j vi , j u i 1, j u i , j 1, j u i , j 4 Re y
F 1 u 2 F 2 v 2
通量 F
1
1 u , Re x 1 v , Re x
G 1 uv G 2 uv
1 u Re y 1 v Re y
, G 1 分别定义在主控制体的中心和角点，如图所示，并按照如下方法离散
F 1 G 1
M 1
p ，代入修正公式得到 M+1 时刻的速
, v M 1 和压力场 p M 1 。将 M+1 时刻的速度场 u M 1 , v M 1 和压力场 p M 1 作为新的
猜测的速度场和猜测的压力场估计值， 采用上述方法计算下一个时刻的速度场和压力场， 直 到满足收敛条件。 收敛判据
Max bip ,j
1 viM vi*, j ,j
将修正后的速度分量代入离散后的连续性方程，得到压力修正方程
p p aip , j pi , j a p ,q p p ,q bi , j
其中
a a
p i 1, j
y 2 y 2 x 2 x 2 p p p u , ai 1, j u , ai , j 1 v , ai , j 1 v ai , j ai 1, j ai , j ai , j 1
xy M 1 M 2 2 2 M 1 M 1 vi , j vi , j Fi ,2 x 0 j Fi 1, j x Gi , j Gi , j 1 y pi , j 1 pi , j t
其中，数值通量
p
pi, j
x 0
编号为（i,j）的速度修正量 u , v 不仅与压力修正量
p 有关，还与邻近点的速度源自文库正
量有关。SIMPLE 算法的重要假定：速度的改变只与压力的改变有关，忽略邻近点对速度修 正的影响。因而得到如下速度修正量
u i, j vi, j
, G 1 和 F 2 , G 2 的某些项冻结于 M 时间层，使离散化之后的方程对
u M 1 , v M 1 是线性的。将离散化之后的 F 1 , G 1 和 F 2 , G 2 代入离散后的 x 方向和 y 方
向的动量方程，整理之后得离散后的动量方程如下
1 M 1 u M 1 M 1 aiu, j uiM au y 0 ,j p , q u p , q bi , j pi 1, j pi , j 1 M 1 v aiv, j viM av ,j p , q v p , q bi , j M 1 i , j 1 1 piM ,j
以上是 SIMPLE 算法中离散化的动量方程
三、SIMPLE 算法基本思想
SIMPLE 算法是一种解决压力-速度耦合问题的“半隐式”算法。首先给定 M 时刻猜测 的速度场 u 估计值 值u
* M
, v M ，用于计算离散动量方程中的系数和常数项。给定 M+1 时刻猜测的压力场
p * ，迭代求解离散动量方程，得到 M+1 时刻速度场的估计值 u * , v * ，速度场的估计
y pi1, j pi, j aiu, j x pi, j 1 pi, j aiv, j
修正后的速度分量
1 u iM u i*, j ,j
y pi1, j pi, j aiu, j x pi, j 1 pi, j aiv, j
1 2 2 xy 1 M M M M aiu, j y uiM x viM 1, j u i 1, j , j vi 1, j vi , j 1 vi 1, j 1 Re x Re y t 4 4
图 2 交错网格示意图
图 3 求解域离散示意图 图 3 中阴影部分代表方腔内的流动区域，阴影区域的边界代表方腔的上、下、左、右 壁面， 阴影区域外面的网格节点是为边界处理需要而设定的虚拟网格节点， 后面介绍边界处 理方法时详细论述。
方程离散
无量纲化的守恒型不可压缩 N S 方程为
U 0 U 1 2 UU P U 0 t Re
本算例采用求解不可压缩流动的经典算法， 即 SIMPLE 算法， 求解方腔内粘性不可压缩 流体运动的定常解。 SIMPLE 算法的全称为 Semi-Implicit Method for Pressure-Linked Equations， 即求解压力关联方程的半隐式算法。 采用 SIMPLE 算法时， 为了避免中心差分格式将 “棋盘” 型参量分布误认为是均匀分布， 需要用交错网格对计算域进行离散。
交错网格
交错网格如图 2 所示，压力、密度等物理量存储在控制体 i, j 的中心，这个控制体称 为主控制体。速度分量
u, v 分别存储在主控制体的 i 1/ 2, j 和 i, j 1/ 2
位置处，
标记为 i, j 位置，再分别以此为中心，划分速度分量 u、v 的控制体。采用空间均匀网格， 等间距离散整个求解域，如图 3 所示。
, G 2 分别定义在主控制体的中心和角点，如图所示，并按照如下方式离散
F 2 G 2
通量 F
1
1 M 1 M 1 M 1 1 M 1 vi 1, j viM viM , j vi 1, j vi , j 1, j vi , j 4 Re y 1 1 M M 1 M 1 1 M 1 u iM viM 1, j u i , j vi 1, j vi , j 1, j vi , j 4 Re x
biv, j
xy M vi , j t
1 1 1 1 M M aiv, j 1 x viM , aiv, j 1 x viM , j vi , j 1 , j 1 vi , j Re y Re y 4 4 1 1 1 1 M M aiv1, j y u iM , aiv1, j y u iM , j 1 u i , j 1, j 1 u i 1, j Re x Re x 4 4
其中， u , v 和
p 分别速度和压力的修正量，修正量亦满足离散的动量方程
u aiu, j ui, j a u p ,q u p ,q bi , j pi 1, j pi , j y 0 v aiv, j vi, j a v p ,q v p ,q bi , j i , j 1
SIMPLE 算法求解方腔内粘性不可压流动
一、问题描述
假设 0 x, y 1 的方腔内充满粘性不可压缩流体，左、右、下壁固定，上壁以
u 16 x 2 1 x 2 运动，试求 Re 100, 200, 400 时的定常解，方腔如图 1 所示。
图 1 方腔内流动示意图
二、离散格式
p i, j
y 2 y 2 x 2 x 2 u u v v ai , j ai 1, j ai , j ai , j 1
* * * * bip , j y u i , j u i 1, j x vi , j vi , j 1
采用迭代法求解压力修正方程，得到压力修正量 度场 u
, v * 满足如下离散方程。
* u * * aiu, j ui*, j a u p ,q u p ,q bi , j pi 1, j pi , j y 0 * v aiv, j vi*, j a v p ,q v p ,q bi , j * i , j 1
p
pi*, j
x 0
一般地， 速度场 u
*
* * * , v * 不满足离散的连续性方程， 因而需要对速度场 u , v 和压力场 p
进行修正。M+1 时刻的修正值和估计值有如下关系
u M 1 u * u v M 1 v * v p M 1 p * p
边界条件处理
首先对计算区域离散，并流动边界之外扩充一个虚拟网格，将真实流动的离散域包围， 如图 7 所示。
图 7 虚拟网格扩充示意图
为很小的正实数，视计算的精度要求而定。本算例中取 1e 8 。
若 bi , j 0 ，则 pi, j 0 ，此时 U i , j
p M 1
U i*, j ，从而来自于离散动量的 U i*, j 满足离散的连
续性方程。因此 Max
b
p i, j
可以作为收敛判据。
xy M 1 M 1 1 1 M 1 M 1 ui , j ui , j Fi ,1j Fi y 0 1, j y Gi , j Gi , j 1 x pi 1, j pi , j t
Y 方向动量方程在速度 v 控制体上离散，时间采用前差
1 2 2 xy 1 M M M M aiv, j x viM y uiM , j 1 u i , j 1 , j 1 u i , j u i 1, j 1 u i 1, j Re y Re x t 4 4
其积分形式为
n UdS 0
S
V
V
u 1 dV n UudS pn x dS n udS 0 S S t Re S v 1 n vdS 0 dV n UvdS pn y dS S S t Re S
p
x 0
其中
biu, j
xy M ui, j t
1 u 1 1 1 M M aiu1, j y u iM , ai -1, j y u iM 1, j u i , j , j u i -1, j Re x Re x 4 4 1 1 1 u 1 M M aiu, j 1 x viM , ai , j 1 x viM , j vi 1, j , j 1 vi 1, j 1 Re y Re y 4 4
图 4 主控制体
图 5 速度 u 控制体
图 6 速度 v 控制体
采用有限体积法离散 N S 方程，连续性方程在主控制体上离散
u
M 1 i, j
1 M 1 1 uiM viM 1, j y vi , j , j 1 x 0
X 方向动量方程在速度 u 控制体上离散，时间采用前差