方腔驱动流源代码(SIMPLE算法)
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
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
其积分形式为
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
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 控制体上离散，时间采用前差
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
, 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
以上是 SIMPLE 算法中离散化的动量方程
三、SIMPLE 算法基本思想
SIMPLE 算法是一种解决压力-速度耦合问题的“半隐式”算法。首先给定 M 时刻猜测 的速度场 u 估计值 值u
* M
, v M ，用于计算离散动量方程中的系数和常数项。给定 M+1 时刻猜测的压力场
p * ，迭代求解离散动量方程，得到 M+1 时刻速度场的估计值 u * , v * ，速度场的估计
交错网格
交错网格如图 2 所示，压力、密度等物理量存储在控制体 i, j 的中心，这个控制体称 为主控制体。速度分量
u, v 分别存储在主控制体的 i 1/ 2, j 和 i, j 1/ 2
位置处，
标记为 i, j 位置，再分别以此为中心，划分速度分量 u、v 的控制体。采用空间均匀网格， 等间距离散整个求解域，如图 3 所示。
p
pi, j
x 0
编号为（i,j）的速度修正量 u , v 不仅与压力修正量
p 有关，还与邻近点的速度修正
量有关。SIMPLE 算法的重要假定：速度的改变只与压力的改变有关，忽略邻近点对速度修 正的影响。因而得到如下速度修正量
u i, j vi, j
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
本算例采用求解不可压缩流动的经典算法， 即 SIMPLE 算法， 求解方腔内粘性不可压缩 流体运动的定常解。 SIMPLE 算法的全称为 Semi-Implicit Method for Pressure-Linked Equations， 即求解压力关联方程的半隐式算法。 采用 SIMPLE 算法时， 为了避免中心差分格式将 “棋盘” 型参量分布误认为是均匀分布， 需要用交错网格对计算域进行离散。
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
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
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
边界条件处理
首先对计算区域离散，并流动边界之外扩充一个虚拟网格，将真实流动的离散域包围， 如图 7 所示。
图 7 虚拟网格扩充示意图
图 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 控制体上离散，时间采用前差
, 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
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
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
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
其中ቤተ መጻሕፍቲ ባይዱ 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
为很小的正实数，视计算的精度要求而定。本算例中取 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
SIMPLE算法及计算例子
动量方程：
—ａ（ｐ＝＿ｕ）＋乳（ｐｕＵ）：ｖ．（印“）一字
（｛２）
Ｏｌ
ｏｘ
＿ａ（ｐ『＿ｖ）＋Ｖ．（ｐｖＵ）：ｖ．（ｐｖｖ）一挈
（１３）
ｕｆ
ｏｙ
—ａ（ｐ＿ｗ）＋ｖ．（户ｗ可）＝ｖ．（胛ｗ）一害
（１．４）
０ｌΒιβλιοθήκη Ｏｚ能量方稃：—ａ（ｐ－Ｔ）＋ｖ．（ｐ功）：ｖ．（生ｖ ７＿）＋ｓｙ
（１ ５）
ｏｌ
Ｃ
若考虑湍流运动，采ＪＨ两方程模式，则还有湍流动能ｋ疗科和湍流动能耗散率占方科：
Ｔｈｅ ｅｆｆｅｃｔ ｏｆ ｓｈａｌｌｏｗ ｗａｔｅｒ ｄｅｐｔｈ ｏｎ ａ ｂｉｇ ＦＰＳＯ（Ｆｌｏａｔｉｎｇ Ｐｒｏｄｕｃｔｉｏｎ，Ｓｔｏｒａｇｅ ａｎｄ Ｏｆｆｉｏａｄｉｎｇ Ｓｙｓｔｅｍ）ｉｓ ｎｕｍｅｒｉｃａｌｌｙ ｓｉｍｕｌａｔｅｄ ｉｎ ｃａｓｅ ｔｈｒｅｅ．Ｉｔ ｈａｓ ｂｅｅｎ ｓｈｏｗｎ ｔｈａｔ ｔｈｅ ｖｉｓｃｏｕｓ ｐｌａｙｓ ａ ｖｅｒｙ＇ｓｍａｌｌ ｒｏｌｅ ａｎｄ ｓｏ ｉｔ ｃａｒｌ ｂｅ ｎｅｇｌｅｃｔｅｄ ｉｎ ｔｈｉｓ ｃｏｎｄｉｔｉｏｎ．Ｔｈｅ ｖａｒｙｉｎｇ ｌａｗ ｏｆ ｐｒｅｓｓｕｒｅ ｐｒｏｄｕｃｅｄ ｂｙ ｓｈａｌｌｏｗ ｗａｔｅｒ ｔｏ ＦＰＳＯ ｉｓ ｇｉｖｅｎ．
Ｔｈｅ ｎｕｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ ｏｆ ｃｕｒｒｅｎｔ ｇｅｎｅｒａｔｉｎｇ ｓｙｓｔｅｍ ｏｆ ｔｈｅ ｏｃｅａｎ ｅｎｇｉｎｅｅｒｉｎｇ ｂａｓｉｎ ｉｓ ｇｉｖｅｎ ｉｎ ｃａｓｅ ｔｏｕｒ ａｎｄ ｓｏｍｅ ｕｓｅｆｕｌ ａｄｖｉｃｅ ｉｓ ｇｉｖｅｎ ｔｏ ａｉｄ ｔｈｅ ｄｅｓｉｇｎ ｏｆ ｃｕｒｒｅｎｔ ｇｅｎｅｒａｔｉｎｇ ｓｙｓｔｅｍ
上海交通大学 博士后士学位论文 SIMPLE算法及应用
4.2 顶盖驱动流(8)
4.2顶盖驱动流4.2.1物理模型在一个正方形的二维空腔中充满等密度的空气，方腔每边长为0.1m，其顶板以0.1m/s 的速度向右移动，同时带动方腔内流体的流动，流场内的流体为层流。
计算区域示意图如图4-2-1所示。
图4-2-1 计算区域示意图4.2.2在Gambit中建立模型Step1：启动Gambit并选择求解器为Fluent5/6。
Step2：创建面操作：→→打开对话框如图4-2-2所示。
输入长度和宽度10，在Direction中选择XY Centered。
图4-2-2 创建面设置对话框Step3：划分面网格操作：→→打开对话框如图4-2-3所示，Shift+鼠标左键选择正方形面，Internal size=0.5，其它保留默认，点击Apply确认。
划分后的网格如图4-2-4所示。
图4-2-3 网格划分设置对话框图4-2-4 计算区域网格图Step4：设置边界类型操作：→●在Name栏输入边界名称wall-1，将Type栏选为Wall，在Entity栏选取Edges，并选中方腔顶部边线。
●在Name栏输入边界名称wall-2，将Type栏选为Wall，在Entity栏选取Edges，并选中方腔其它三条边线。
Step5：输出网格文件操作：Fil m→export→mesh打开对话框如图4-2-5所示，选中Export 2-D mesh 前面的复选框，输出网格文件。
图4-2-5 网格文件输出对话框4.2.3求解计算Step1：启动Fluent选择2d单精度求解器，点击Run，如图4-2-6所示。
图4-2-6 启动求解器图4-2-7 网格尺寸设置对话框Step2：导入并检查网格1．读入网格文件操作：Fil e→Read→Case...找到文件后，单击OK按键确认。
2．检查网格操作：Grid→Check网格读入后，一定要进行网格检查，注意最小体积不能为负值。
3．网格比例设置操作：Grid→Scale...在Gambit中，生成网格使用的单位是cm，在Grid Was Created In下拉菜单中，选取cm，如图4-2-7所示，然后单击Scale，关闭对话框。
顶盖驱动流的格子Boltzmann模拟代码
for(i=1;i<NX;i++) for(j=1;j<NY;j++) for(k=0;k<Q;k++) { ip=i-e[k][0]; jp=j-e[k][1]; F[i][j][k]=f[ip][jp][k]+(feq(k,rho[ip][jp], u[ip][jp])-f[ip][jp][k])/tau_f; } for(i=1;i<NX;i++) for(j=1;j<NY;j++) { u0[i][j][0]=u[i][j][0]; u0[i][j][1]=u[i][j][1]; rho[i][j]=0; u[i][j][0]=0; u[i][j][1]=0; for(k=0;k<Q;k++) { f[i][j][k]=F[i][j][k]; rho[i][j]+=f[i][j][k]; u[i][j][0]+=e[k][0]*f[i][j][k]; u[i][j][1]+=e[k][1]*f[i][j][k]; } u[i][j][0]/=rho[i][j]; u[i][j][1]/=rho[i][j]; } for(j=1;j<NY;j++) for(k=0;k<Q;k++) { rho[NX][j]=rho[NX-1][j]; f[NX][j][k]=feq(k,rho[NX][j],u[NX][j])+(f[NX-1] [j][k]-feq(k,rho[NX-1][j],u[NX-1][j]))*(1-1/tau_f); rho[0][j]=rho[1][j]; f[0][j][k]=feq(k,rho[0][j],u[0][j])+(f[1][j][k] -feq(k,rho[1][j],u[1][j]))*(1-1/tau_f); } for(i=0;i<=NX;i++) for(k=0;k<Q;k++) { rho[i][0]=rho[i][1]; f[i][0][k]=feq(k,rho[i][0],u[i][0])+(f[i][1][k] -feq(k,rho[i][1],u[i][1]))*(1-1/tau_f); rho[i][NY]=rho[i][NY-1]; u[i][NY][0]=U;
模糊控制技术在SIMPLER算法中的应用及求解性能分析
模糊控制技术在SIMPLER算法中的应用及求解性能分析王艳宁;孙东亮;苗政;陈家庆;蔡晓君【摘要】为了提高SIMPLER算法在三维流动问题上的求解性能,引入模糊控制方法来自动调控速度亚松弛因子的大小.在数值计算过程中,将相邻两个迭代层次上的最大动量残差比值作为模糊控制输入量,速度亚松弛因子的变化量作为模糊控制输出量,基于最大动量残差的变化趋势可实现速度亚松弛因子的自动调控,从而达到加快收敛的目的.最后,通过3个经典的流动问题验证了模糊控制方法的优越性.研究表明:当初始亚松弛因子为最不利值时,模糊控制方法的收敛速度约是固定松弛因子方法的5～30倍;当初始亚松弛因子为最佳值时,模糊控制方法迭代次数与固定松弛因子方法迭代次数之比为0.7～2.0,收敛速度相差不大;采用模糊控制方法后,SIMPLER 算法在不同初始亚松弛因子下均能得到高速收敛的解,同时健壮性也显著提高.研究工作将为大幅提升SIMPLER算法在三维流动问题上的求解性能起到重要作用.【期刊名称】《西安交通大学学报》【年(卷),期】2016(050)001【总页数】7页(P78-84)【关键词】模糊控制;SIMPLER算法;三维流动问题;亚松弛因子;收敛速度;健壮性【作者】王艳宁;孙东亮;苗政;陈家庆;蔡晓君【作者单位】华北电力大学可再生能源学院,102206,北京;北京石油化工学院机械工程学院,102617,北京;华北电力大学可再生能源学院,102206,北京;北京石油化工学院机械工程学院,102617,北京;北京石油化工学院机械工程学院,102617,北京【正文语种】中文【中图分类】TK124求解流动与传热问题的压力修正算法首次由著名学者Patankar和Spalding提出,并被命名为SIMPLE算法[1]。
为了克服SIMPLE算法中初始压力和速度场单独进行设定的缺点,随后Patankar提出了改进算法SIMPLER[2]。
目前,SIMPLE系列算法已被广泛应用于求解流动与传热问题[3-7]。
基于MPI+FreeFem++的有限元并行计算
基于MPI+FreeFem++的有限元并行计算摘要：有限元方法是一种灵活而高效的数值求解偏微分方程的计算方法，是工程分析和计算中不可缺少的重要工具之一。
在计算机技术的快速发展使得并行机的价格日益下降的今天，并行有限元计算方法受到了学术界和工程界的普遍关注。
讨论了基于MPI+FreeFem++的有限元并行计算环境的构建，阐述了在该环境下有限元并行程序的编写、编译及运行等过程，并通过具体编程实例，说明了MPI+FreeFem++环境下的有限元并行编程的简单和高效。
关键词：有限元方法；并行计算；MPI；FreeFem++0 引言有限元方法是20世纪50年代伴随电子计算机的诞生，在计算数学和计算工程领域里诞生的一种高效而灵活的计算方法，它将古典变分法与分片多项式插值相结合，易于处理复杂的边值问题，具有有限差分法无可比拟的优越性，广泛应用于求解热传导、电磁场、流体力学等相关问题，已成为当今工程分析和计算中不可缺少的最重要的工具之一。
有限元方法的“化整为零、积零为整”的基本思想与并行处理技术的基本原则“分而治之”基本一致，因而具有高度的内在并行性。
在计算机技术快速发展使得并行机的价格日益下降的今天，有限元并行计算引起了学术界和工程界的普遍关注，吸引了众多科研与工程技术人员。
但要实现有限元并行编程，并不是一件容易的事，特别是对于复杂区域问题，若从网格生成、任务的划分、单元刚度矩阵的计算、总刚度矩阵的组装，到有限元方程组的求解以及后处理，都需要程序员自己编写代码的话，将是一件十分繁琐的事情。
本文探讨了构建基于MPI+FreeFem++的有限元并行计算环境，在该环境下，程序员可避免冗长代码的编写，进而轻松、快速、高效地实现复杂问题的有限元并行计算。
1 FreeFem++简介FreeFem++ 是一款免费的偏微分方程有限元计算软件，它集成网格生成器、线性方程组的求解器、后处理及计算结果可视化于一体，能快速而高效地实现复杂区域问题的有限元数值计算。
方腔驱动流源代码(SIMPLE算法)
p
pi, j
x 0
编号为（i,j）的速度修正量 u , v 不仅与压力修正量
p 有关，还与邻近点的速度修正
量有关。SIMPLE 算法的重要假定：速度的改变只与压力的改变有关，忽略邻近点对速度修 正的影响。因而得到如下速度修正量
u i, j vi, j
本算例采用求解不可压缩流动的经典算法， 即 SIMPLE 算法， 求解方腔内粘性不可压缩 流体运动的定常解。 SIMPLE 算法的全称为 Semi-Implicit Method for Pressure-Linked Equations， 即求解压力关联方程的半隐式算法。 采用 SIMPLE 算法时， 为了避免中心差分格式将 “棋盘” 型参量分布误认为是均匀分布， 需要用交错网格对计算域进行离散。
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 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
以上是 SIMPLE 算法中离散化的动量方程
三、SIMPLE 算法基本思想
SIMPLE 算法是一种解决压力-速度耦合问题的“半隐式”算法。首先给定 M 时刻猜测 的速度场 u 估计值 值u
力学流体simple算法完成
力学流体simple算法完成一、简介力学流体（Computational Fluid Dynamics，CFD）是一种通过数值模拟来研究流体运动的方法。
它可以帮助工程师和科学家了解流体如何在各种条件下移动和相互作用，以及如何在不同的环境中优化设计。
其中simple算法是一种常用的求解Navier-Stokes方程的方法，本文将对其进行详细介绍。
二、simple算法基础1. Navier-Stokes方程Navier-Stokes方程是描述流体运动的基本方程之一，它包含了三个主要部分：质量守恒方程、动量守恒方程和能量守恒方程。
其中，质量守恒方程描述了质量在空间和时间上的变化；动量守恒方程描述了力对流体速度的影响；能量守恒方程描述了温度在空间和时间上的变化。
2. 离散化方法为了求解Navier-Stokes方程，需要将其离散化为有限个离散点上的代数形式。
这可以通过有限差分、有限元或谱方法等技术实现。
其中，有限差分法是最常用的离散化方法之一。
3. simple算法原理simple算法是一种求解Navier-Stokes方程的迭代方法，它采用了分步求解的策略。
simple算法包括三个主要步骤：预处理、压力修正和速度修正。
其中，预处理步骤用于确定初始速度场和压力场；压力修正步骤用于计算压力场的变化；速度修正步骤用于更新速度场。
通过迭代这三个步骤，可以逐渐收敛到稳定的解。
三、simple算法实现1. 离散化首先需要将Navier-Stokes方程离散化为有限差分形式。
例如，在二维情况下，可以将动量守恒方程离散化为：$ \frac{u_{i+1,j}^{n+1}-u_{i,j}^{n}}{\Delta t}+\frac{u_{i,j+1}^{n}-u_{i,j-1}^{n}}{2\Delta y}-\frac{\nu}{\Delta x^2}(u_{i+1,j}^{n+1}-2u_{i,j}^{n+1}+u_{i-1,j}^{n+1})=-\frac{p_{i+1,j}-p_{i,j}}{\Delta x}$其中，$u$表示水平方向的速度，$v$表示垂直方向的速度，$p$表示压力，$\nu$表示粘性系数。
SIMPLE算法及计算例子
第二个问题：压力的一阶导数是以源项的形式出现在动簧方程中。采用分离式求解各 变量的离散方程时，由于压力没有独立的方程，压力与速度的关系隐含于连续性方程中，如 果压力场是正确的，则据此压力场求得的速度场一定满足连续性方程。如何构造求解压力场 的方程，或者说在假定初始乐力分布＿｜亓如何构造计算压力改进值得方程，就成了分离式求解 法中的一个关键问题。
在交错网格系统中，关『．“、ｖ的离散方程可通过对Ｕ、ｖ番白的控制容积作积分而 得出。这时压力梯度的离散形式对虬为（ＰＦ—Ｐ，，）／ｄｘ，对ｋ为（Ｐ．。一Ｐ，，）／ａｙ，亦即相邻 两点间的压力羞构成了印／彘、印／砂，这就从根本上解决了采川一般网格系统时所遇到 的困难，这也是交错网格的成功之处。
Ｋｅｙｗｏｒｄｓ：ＳＩＭＰＬＥ ａｌｇｏｒｉｔｈｍ；ｔｕｒｂｕｌｅｎｃｅ ｍｏｄｅｌｉｎｇ；ｇｒａｖｉｔｙ ｃｕｒｒｅｎｔ； ｎｕｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ；ｃｏｍｐｕｔａｔｉｏｎａｌ ｆｌｕｉｄ ｍｅｃｈａｎｉｃｓ
上海交通大学博十后出站报告
第一章ＳＩＭＰＬＥ算法概述
SIMPLEC 二维、结构网格单相流SIMPLE 算法的C 语言实现的说明
v ∂ ( ρφ ) + div ( ρ uφ ) = div( Γgradφ ) + S （传热 ２．１３） ∂t
（热图 ２．１） 在一个控制容积上的流量平衡 此微元内的传递行为可以用下面的式子统一地表示：
∂ ( ρφ ) + divJ φ = S （１） ∂t ∂ ∂J ∂J y ∂J z ( ρφ ) 代表单位容积内 φ 的变化率。divJ φ = x + + ∂ ∂x ∂y ∂z
关于文中的公式与图表的引用为了方便读者与传热与流体流动的数值计算的对照凡是从原著上摘抄的公式与图表均保持原著中的编号
二维、结构网格单相流 ＳＩＭＰＬＥ 算法的 Ｃ 语言实现的说明 ｑｚｈｅｎ＠ｆｌｏｔｕ．ｏｒｇ Ｐａｇｅ １／４９
第一节 前言
ＳＩＭＰＬＥ 算法的全称是“Ｓｅｍｉ－Ｉｍｐｌｉｃｉｔ Ｍｅｔｈｏｄ ｆｏｒ Ｐｒｅｓｓｕｒｅ－Ｌｉｎｋｅｄ Ｅｑｕａｔｉｏｎｓ”，意思是求解压力耦合 的质量／动量／能量传递方程的半隐式方法。自从 Ｓ．Ｖ．Ｐａｔａｎｋａｒ（帕坦卡）和 Ｄ．Ｂ．Ｓｐａｌｄｉｎｇ（斯伯丁）于 １９７２ 年提出 ＳＩＭＰＬＥ 算法后，它实际上已经成为许多工程流动、传热以及反应体系的数值模拟的最重要的方法。 许多商业 ＣＦＤ 软件，如 ｃｆｘ 与 ｆｌｕｅｎｔ，其核心也都基于 ＳＩＭＰＬＥ 算法。 了解 ＳＩＭＰＬＥ 算法是学习差分数值计算以及 ＣＦＤ（计算流体力学）一个非常好的起点。而且，ＳＩＭＰＬＥ 算法对 于很多工程问题也足够用了。介绍 ＳＩＭＰＬＥ 算法最好的书自然是帕坦卡自己所著的《传热与流体流动的数值 计算》［１］。不过，这本优秀、简洁的著作没有附带源代码。 我在学习 ＳＩＭＰＬＥ 的过程中，了解到一份可靠、简单、说明文档完整、可方便地编译执行的源代码对于初学 者的重要性。于是我用 Ｃ 语言改写了 ＳＩＭＰＬＥ 算法（称之为 ＳＩＭＰＬＥ＿Ｃ），并给出了详细的文档，希望能够缩 短算法与代码之间的距离。算法源自水木清华 ＢＢＳ 上 ｏｒａｎｇｅ０１ 与 ｍｏｒａｎ 两位同学所给予的 ＦＯＲＴＲＡＮ ＳＩＭＰＬＥ 代码。同时清华大学化工系反应工程教研组刘飞同学的笔记也给了我很大的帮助，在此一并谢过。 目前，ＳＩＭＰＬＥ＿Ｃ 只能解决二维直角坐标、结构网格的问题。取决于表征问题的微分方程，ＳＩＭＰＬＥ＿Ｃ 在这个 几何框架下能够模拟瞬态或者时均的流动、传热以及反应耦合的单相流体系。应该说明，下面的一些问题 可能或者根本不在 ＳＩＭＰＬＥ＿Ｃ 的能力范围之内：高精度的流场预测；多相流；三维流动；复杂流场或者几何 结构；边界层行为预测；马赫数接近或超过 １ 的高速流动。 关于 ｏｒａｎｇｅ０１ 与 ｍｏｒａｎ 两位同学的 ＳＩＭＰＬＥ ＦＯＲＴＲＡＮ 代码的版权问题，我在源代码上看见这么一段，在此 原封不动地摘录如下： Ｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃ Ｃｃｃｃ Ｔｈｉｓ ｃｏｍｐｕｔｅｒ ｐｒｏｇｒａｍ ｗａｓ ｃｏｐｉｅｄ ｆｒｏｍ ｔｈｅ ｇｒａｄｕａｔｅ ｓｔｕｄｅｎｔ ｃｏｕｒｓｅ ｐｒｏｇｒａｍ Ｃｃｃｃ ｏｆ ｔｈｅ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｍｉｎｎｅｓｏｔａ． Ｐａｒｔ ｏｆ ｉｔ ｗａｓ ｒｅ－ｆｏｒｍｕｌａｔｅｄ ｔｏ ｍｅｅｔ ｔｈｅ Ｃｃｃｃ ｐｅｒｓｏｎａｌ ｃｏｍｐｕｔｅｒ ｅｎｖｉｒｏｎｍｅｎｔ． Ｓｏｍｅ ｉｎａｐｐｒｏｐｒｉａｔｅ ｅｘｐｒｅｓｓｉｏｎｓ ｗｅｒｅ ａｌｓｏ Ｃｃｃｃ ｃｏｒｒｅｃｔｅｄ． Ｔｈｅ ｐｒｏｇｒａｍ ｉｓ ｕｓｅｄ ｏｎｌｙ ｆｏｒ ｔｅａｃｈｉｎｇ ｐｕｒｐｏｓｅ． Ｎｏ ｐａｒｔ ｏｆ ｉｔ Ｃｃｃｃ ｍａｙ ｂｅ ｐｕｂｌｉｓｈｅｄ． Ｙｏｕ ｍａｙ ｕｓｅ ｉｔ ａｓ ａ ｆｒａｍｅ ｔｏ ｒｅ－ｄｅｖｅｌｏｐ ｙｏｕｒ ｏｗｎ ｃｏｄｅ Ｃｃｃｃ ｆｏｒ ｒｅｓｅａｒｃｈ ｐｕｒｐｏｓｅ． ＸＪＴＵ Ｉｎｓｔｒｕｃｔｏｒ， １９９５．１１ Ｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃｃ 本文的工作是非常粗浅的。本质上它是帕坦卡在７０年代就已解决的问题的一个简单的复述。不过，希望它 能起到抛砖引玉的作用――期望更多的研究者能够将手头基于ＳＩＭＰＬＥ的先进算法的源代码整理并共享出 来。下面一些特性是特别欢迎的： l l 三维网格； 非结构／自适应网格； 高精度差分格式； 欧拉＋欧拉模式的多相流； 欧拉＋拉格朗日模式的多相流； 基于ＰＣ ＣＬＵＳＴＥＲ（个人电脑集群，ＭＩＭＤ体系）的并行算法。
以SIMPLE算法求解瞬态二维顶板驱动方腔流动问题
P旳
(lb)
也+也 =0
dx dy
(lc)
式(la) ~ (lc)中,t一时间,s；p—流体压力,Pa；
u、v一 ％,y方向上的速度,m/s；a—运动粘性系数,
m2/s.
2问题描述
如图1所示的方腔内充满粘性不可压缩流体, 运动粘性系数U在方腔内为常数，其左、右、下壁固 定,上壁以U = S运动.求在这样一种驱动力作用 下，方腔内的流体动力特性•在图1中，方腔的高度
1基本控制方程
对于二维顶板驱动方腔不可压均质粘性瞬态
流动问题,其流场的流体控制方程由式(la) ~(lc)
给出•
=£(u +il du + d( uu) +。(他)
dt dx
du\
丄型
~dx]
i dy) p dx
(la)
dv dt
uv) dx
1
d(仞)
=— d / V dx\
dx]
+斗 v^\ dy\ dyl
(9)
5方程求解
对于本算例的瞬态问题，在每一时间步采用
SIMPLE算法来求解时，需要计算在t + △力时刻的:
1）速度分量叫八％的初算值；2）皱八 %•的修正值此广、巧：阳；3）压力刃J的修正值
刃广•在交错网格的编码系统中，节点GJ）和（/J）
点处的速度分量初算值＜/+1、吩小和压力修正方
程的隐式格式离散方程可表达为：
\-VI,N + 2 -^VI,N+1 ~VI,N
(1 = 1,...⑻
4)压力和压力修正量物理边界条件的插值
为配合采用SIMPLE算法求解，需要将压力P
和压力修正值P'扩充到相应的虚拟网格节点中去. 对于本算例，由于所有垂宜于固壁的速度分量均为