说明：程序旨在熟悉VOF求解流程，投影法求解Matlab编程主要参考Mark Owkes的Guide To CFD资料，对于VOF模型的对流项需要采取特殊的离散方式，如：施主-受主格式，但本文中仅采取中心差分格式进行离散，会带来较大误差，其次对于每个时间步都对体积分数值做矫正，限制其在0-1之间。后续如果使用了更精确的格式，再来更新。
摘要
本文介绍基于matlab平台的VOF编程实例，采用投影法求解二维不可压缩N-S方程，交错网格结合中心差分法离散控制方程。并对Re=100方腔流进行计算验证，证明投影法求解程序的正确性，随后对二维溃坝算例进行数值计算。
关键词：数值模拟，VOF，Matlab，溃坝
第一章 数值离散
1.1投影法求解原理
投影法是求解不可压缩N-S方程的常用方法之一，对于二维不可压缩N-S方程以及连续性方程如公式(1-1)展开形式如公式(1-2)
{ ∂ u ∂ t + u ⋅ ∇ u = − 1 ρ ∇ p + υ ∇ 2 u ∇ ⋅ u = 0 (1 - 1) \begin{matrix} \left\{ \begin{matrix} \frac{\partial\mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla\mathbf{u =} - \frac{1}{\rho}\nabla p + \upsilon\nabla^{2}\mathbf{u} \\ \nabla \cdot \mathbf{u = 0} \\ \end{matrix} \right.\ \\ \end{matrix}\tag{1 - 1} {∂t∂u+u⋅∇u=−ρ1∇p+υ∇2u∇⋅u=0 (1 - 1)
展开形式
∂ u ∂ t + u ∂ u ∂ x + v ∂ u ∂ y = − 1 ρ ∇ p + υ ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 ) ∂ v ∂ t + u ∂ v ∂ x + v ∂ v ∂ y = − 1 ρ ∇ p + υ ( ∂ 2 v ∂ x 2 + ∂ 2 v ∂ y 2 ) ∂ u ∂ x + ∂ v ∂ y = 0 (1 - 2) \begin{matrix} \begin{matrix} \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = - \frac{1}{\rho}\nabla p + \upsilon\left( \frac{\partial^{2}u}{\partial x^{2}} + \frac{\partial^{2}u}{\partial y^{2}} \right) \\ \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = - \frac{1}{\rho}\nabla p + \upsilon\left( \frac{\partial^{2}v}{\partial x^{2}} + \frac{\partial^{2}v}{\partial y^{2}} \right) \\ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \\ \end{matrix}\\ \end{matrix}\tag{1 - 2} ∂t∂u+u∂x∂u+v∂y∂u=−ρ1∇p+υ(∂x2∂2u+∂y2∂2u)∂t∂v+u∂x∂v+v∂y∂v=−ρ1∇p+υ(∂x2∂2v+∂y2∂2v)∂x∂u+∂y∂v=0(1 - 2)
投影法原理可以分为以下三步，将时间推进分成三个子步，在中间步解出压力，首先需要确定时间离散格式，这里时间离散化采用显式欧拉格式，则离散方程如（1- 3）所示。
u n + 1 − u n Δ t = − 1 ρ ∇ p n + 1 − u n ⋅ ∇ u n + υ ∇ 2 u n (1 - 3) \begin{matrix} \frac{\mathbf{u}^{n + 1} - \mathbf{u}^{n}}{\mathrm{\Delta}t}\mathbf{=} - \frac{1}{\rho}\nabla p^{n + 1} - \mathbf{u}^{n} \cdot \nabla\mathbf{u}^{n} + \upsilon\nabla^{2}\mathbf{u}^{n}\ \\ \end{matrix}\tag{1 - 3} Δtun+1−un=−ρ1∇pn+1−un⋅∇un+υ∇2un (1 - 3)
将方程（1-3）分解为个方程（1-4）与（1-5）
预算步 Step1：
u ∗ − u n Δ t = − u n ⋅ ∇ u n + υ ∇ 2 u n (1 - 4) \begin{matrix} \frac{\mathbf{u}^{\mathbf{*}} - \mathbf{u}^{n}}{\mathrm{\Delta}t}\mathbf{=} - \mathbf{u}^{n} \cdot \nabla\mathbf{u}^{n} + \upsilon\nabla^{2}\mathbf{u}^{n}\\ \end{matrix}\tag{1 - 4} Δtu∗−un=−un⋅∇un+υ∇2un(1 - 4)
压力修正步 Step2：
u n + 1 − u ∗ Δ t = − 1 ρ ∇ p n + 1 (1 - 5) \begin{matrix} \frac{{\mathbf{u}^{n + 1}\mathbf{- u}}^{\mathbf{*}}}{\mathrm{\Delta}t}\mathbf{=} - \frac{1}{\rho}\nabla p^{n + 1}\ \\ \end{matrix}\tag{1 - 5} Δtun+1−u∗=−ρ1∇pn+1 (1 - 5)
对于方程（1-5），并不满足连续性方程，所以需要对两边求散度，并且需要满足条件： ∇ ⋅ u n + 1 = 0 \nabla\cdot \mathbf{u}^{n + 1} = 0 ∇⋅un+1=0则方程（1-5）变化为：
∇ 2 p n + 1 = ρ ∇ ⋅ u ∗ Δ t (1 - 6) \begin{matrix} \nabla^{2}p^{n + 1}\mathbf{=}\frac{\rho{\nabla \cdot \mathbf{u}}^{\mathbf{*}}}{\mathrm{\Delta}t}\ \\ \end{matrix}\tag{1 - 6} ∇2pn+1=Δtρ∇⋅u∗ (1 - 6)
最终步
Step3：由中间速度
u
∗
\mathbf{u}^{\mathbf{*}}
u∗与压力计算得到下一刻时间的速度场。至此完成了求解随着时间推进的过程。
u n + 1 − u ∗ Δ t = − 1 ρ ∇ p n + 1 (1 - 7) \begin{matrix} \frac{{\mathbf{u}^{n + 1}\mathbf{- u}}^{\mathbf{*}}}{\mathrm{\Delta}t}\mathbf{=} - \frac{1}{\rho}\nabla p^{n + 1}\\ \end{matrix}\tag{1 - 7} Δtun+1−u∗=−ρ1∇pn+1(1 - 7)
1.2空间离散：
本程序采用有限差分法求解，为了避免棋盘式的压力分布，采用交错网格，压力储存在网格点中心，速度储存在网格左边界与底边界，如图1所示，其中imin与jmin均为2，相当于在边界外层预留有一层网格，用于设置边界条件，后续边界条件一节将会介绍。
1.3方程离散
1.3.1 动量方程离散：
对于方程（1-4）离散，以x方向动量方程为例，展开式如下
u ∗ = u n + Δ t ( − ( u n ∂ u n ∂ x + v n ∂ u n ∂ y ) + υ ( ∂ 2 u n ∂ x 2 + ∂ 2 u n ∂ y 2 ) ) (1 - 8) \begin{matrix} u^{\mathbf{*}}\mathbf{=}u^{n} + \mathrm{\Delta}t\left( - \left( u^{n}\frac{\partial u^{n}}{\partial x} + v^{n}\frac{\partial u^{n}}{\partial y} \right) + \upsilon\left( \frac{\partial^{2}u^{n}}{\partial x^{2}} + \frac{\partial^{2}u^{n}}{\partial y^{2}} \right) \right) \\ \end{matrix}\tag{1 - 8} u∗=un+Δt(−(un∂x∂un+vn∂y∂un)+υ(∂x2∂2un+∂y2∂2un))(1 - 8)
其中，一阶，二阶空间导数项采用中心差分格式离散，离散形式如下，空间位置关系图如图2所示。
1.3.2 压力Poisson 方程离散
∇ 2 p n + 1 = ρ ∇ ⋅ u ∗ Δ t \nabla^{2}p^{n + 1}\mathbf{=}\frac{\rho{\nabla \cdot \mathbf{u}}^{\mathbf{*}}}{\mathrm{\Delta}t} ∇2pn+1=Δtρ∇⋅u∗
压力空间导数项采用中心差分格式
∇ 2 p n + 1 = ∂ 2 p n + 1 ∂ x 2 + ∂ 2 p n + 1 ∂ y 2 = p n + 1 ( i − 1 , j ) − 2 p n + 1 ( i , j ) + p n + 1 ( i + 1 , j ) Δ x 2 + p n + 1 ( i , j − 1 ) − 2 p n + 1 ( i , j ) + p n + 1 ( i , j + 1 ) Δ y 2 \nabla^{2}p^{n + 1} = \frac{\partial^{2}p^{n + 1}}{\partial x^{2}} + \frac{\partial^{2}p^{n + 1}}{\partial y^{2}} = \frac{p^{n + 1}\left( i - 1,j \right) - {2p}^{n + 1}\left( i,j \right) + p^{n + 1}\left( i + 1,j \right)}{\mathrm{\Delta}x^{2}} + \frac{p^{n + 1}\left( i,j - 1 \right) - {2p}^{n + 1}\left( i,j \right) + p^{n + 1}\left( i,j + 1 \right)}{\mathrm{\Delta}y^{2}} ∇2pn+1=∂x2∂2pn+1+∂y2∂2pn+1=Δx2pn+1(i−1,j)−2pn+1(i,j)+pn+1(i+1,j)+Δy2pn+1(i,j−1)−2pn+1(i,j)+pn+1(i,j+1)
∇ ⋅ u ∗ = ∂ u ∗ ∂ x + ∂ v ∗ ∂ y = u ∗ ( i + 1 , j ) − u ∗ ( i , j ) Δ x + v ∗ ( i , j + 1 ) − v ∗ ( i , j ) Δ y {\nabla \cdot \mathbf{u}}^{\mathbf{*}}\mathbf{=}\frac{\partial u^{*}}{\partial x} + \frac{\partial v^{*}}{\partial y} = \frac{u^{*}\left( i + 1,j \right) - u^{*}\left( i,j \right)}{\mathrm{\Delta}x} + \frac{v^{*}\left( i,j + 1 \right) - v^{*}\left( i,j \right)}{\mathrm{\Delta}y} ∇⋅u∗=∂x∂u∗+∂y∂v∗=Δxu∗(i+1,j)−u∗(i,j)+Δyv∗(i,j+1)−v∗(i,j)
将离散后的压力possion方程组写成矩阵形式 L p n + 1 = R \mathbf{L}\mathbf{p}^{n + 1} = \mathbf{R} Lpn+1=R
[ L 11 ⋯ L 1 n ⋮ ⋱ ⋮ L n 1 ⋯ L nn ] ⋅ [ p 1 ⋮ p n ] = [ R 1 ⋮ R n ] (1 - 9 ) \begin{matrix} \begin{bmatrix} L_{11} & \cdots & L_{1n} \\ \vdots & \ddots & \vdots \\ L_{n1} & \cdots & L_{\text{nn}} \\ \end{bmatrix} \cdot \begin{bmatrix} p_{1} \\ \vdots \\ p_{n} \\ \end{bmatrix} = \begin{bmatrix} R_{1} \\ \vdots \\ R_{n} \\ \end{bmatrix}\\ \end{matrix}\tag{1 - 9 } ⎣⎢⎡L11⋮Ln1⋯⋱⋯L1n⋮Lnn⎦⎥⎤⋅⎣⎢⎡p1⋮pn⎦⎥⎤=⎣⎢⎡R1⋮Rn⎦⎥⎤(1 - 9 )
L p n + 1 \mathbf{L}\mathbf{p}^{n +1} Lpn+1的系数矩阵，系数矩阵只与网格划分相关可以在网格划分完成之后的时候建立 L \mathbf{L} L 矩阵。应用matlab编程时网格数量较小，可以采用矩阵除法进行求解，如方程（1-10）所示。
p n + 1 = L 左 除 R (1 - 10) \begin{matrix} \mathbf{p}^{n + 1} = \mathbf{L{左除}R} \\ \end{matrix}\tag{1 - 10} pn+1=L左除R(1 - 10)
最终步，根据求得的压力与中间速度求出下一时刻的速度场，方程与离散格式如下：
u n + 1 − u ∗ Δ t = − 1 ρ ∇ p n + 1 \frac{{\mathbf{u}^{n + 1}\mathbf{- u}}^{\mathbf{*}}}{\mathrm{\Delta}t}\mathbf{=} - \frac{1}{\rho}\nabla p^{n + 1} Δtun+1−u∗=−ρ1∇pn+1
u n + 1 = u ∗ − Δ t ρ ∂ p ∂ x , v n + 1 = v ∗ − Δ t ρ ∂ p ∂ y u^{n + 1}\mathbf{=}u^{\mathbf{*}} - \frac{\mathrm{\Delta}t}{\rho}\frac{\partial p}{\partial x},v^{n + 1}\mathbf{=}v^{\mathbf{*}} - \frac{\mathrm{\Delta}t}{\rho}\frac{\partial p}{\partial y} un+1=u∗−ρΔt∂x∂p,vn+1=v∗−ρΔt∂y∂p
∂ p ∂ x = p ( i , j ) − p ( i − 1 , j ) Δ x , ∂ p ∂ y = p ( i , j ) − p ( i , j − 1 ) Δ x \frac{\partial p}{\partial x} = \frac{p\left( i,j \right) - p\left( i - 1,j \right)}{\mathrm{\Delta}x},\frac{\partial p}{\partial y} = \frac{p\left( i,j \right) - p\left( i,j - 1 \right)}{\mathrm{\Delta}x} ∂x∂p=Δxp(i,j)−p(i−1,j),∂y∂p=Δxp(i,j)−p(i,j−1)
1.4 边界条件设置
由于速度
u
,
v
u,v
u,v存储在网格单元left，bottom边界，压力值存储在网格中心，为了确定边界上部分速度值需要做向外扩充一层网格，通过插值来确定边界上的速度值。以左边界为例。边界上的水平速度可以直接设定，垂向速度需要通过插值确定。
v left = ( v 1 + v 2 ) / 2 v_{\text{left}} = \left( v1 + v2 \right)/2 vleft=(v1+v2)/2
∂ ∂ n v left = 0 : ( v 1 = v 2 ) {\frac{\partial}{\partial n}v}_{\text{left}} = 0\ \ \ \ \ \ \ \ :\ \left( v1 = v2 \right) ∂n∂vleft=0 : (v1=v2)
v left = 0 : ( v 1 + v 2 ) = 0 v_{\text{left}} = 0\ \ \ \ \ \ \ \ :\ \left( v1 + v2 \right) = 0 vleft=0 : (v1+v2)=0
1.5 求解流程
整体求解流程如图5所示，其中红色方框中便是对应的程序函数。代码见附录
1.6 VOF求解原理以及流程
溃坝算例中需要考虑重力的影响
∂ u ∂ t + u ⋅ ∇ u = − 1 ρ ∇ p + υ ∇ 2 u + g ∇ ⋅ u = 0 ( 1 - 11 ) \begin{matrix} \begin{matrix} \frac{\partial\mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla\mathbf{u =} - \frac{1}{\rho}\nabla p + \upsilon\nabla^{2}\mathbf{u + g} \\ \nabla \cdot \mathbf{u = 0} \\ \end{matrix} \\ \end{matrix}\tag{ 1 - 11 } ∂t∂u+u⋅∇u=−ρ1∇p+υ∇2u+g∇⋅u=0( 1 - 11 )
展开形式
{ ∂ u ∂ t + u ∂ u ∂ x + v ∂ u ∂ y = − 1 ρ ∇ p + υ ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 ) + g x ∂ v ∂ t + u ∂ v ∂ x + v ∂ v ∂ y = − 1 ρ ∇ p + υ ( ∂ 2 v ∂ x 2 + ∂ 2 v ∂ y 2 ) + g y ∂ u ∂ x + ∂ v ∂ y = 0 ( 1 - 12 ) \begin{matrix} \left\{ \begin{matrix} \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = - \frac{1}{\rho}\nabla p + \upsilon\left( \frac{\partial^{2}u}{\partial x^{2}} + \frac{\partial^{2}u}{\partial y^{2}} \right) + g_{x} \\ \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = - \frac{1}{\rho}\nabla p + \upsilon\left( \frac{\partial^{2}v}{\partial x^{2}} + \frac{\partial^{2}v}{\partial y^{2}} \right) + g_{y} \\ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \\ \end{matrix} \right.\ \\ \end{matrix}\tag{ 1 - 12 } ⎩⎪⎪⎨⎪⎪⎧∂t∂u+u∂x∂u+v∂y∂u=−ρ1∇p+υ(∂x2∂2u+∂y2∂2u)+gx∂t∂v+u∂x∂v+v∂y∂v=−ρ1∇p+υ(∂x2∂2v+∂y2∂2v)+gy∂x∂u+∂y∂v=0 ( 1 - 12 )
VOF模型带来的额外方程：体积分数输运方程（1-13）
∂ F ∂ t + u ∂ F ∂ x + v ∂ F ∂ y = 0 ( 1 - 13 ) \begin{matrix} \frac{\partial F}{\partial t} + u\frac{\partial F}{\partial x} + v\frac{\partial F}{\partial y} = 0\\ \end{matrix}\tag{ 1 - 13 } ∂t∂F+u∂x∂F+v∂y∂F=0( 1 - 13 )
密度与粘度采用体积分数构造得到的混合密度表示
ρ = ρ 1 ⋅ F + ( 1 − F ) ρ 2 \rho = \rho_{1} \cdot F + \left( 1 - F \right)\rho_{2} ρ=ρ1⋅F+(1−F)ρ2
μ = μ 1 ⋅ F + ( 1 − F ) μ 2 \mu = \mu_{1} \cdot F + \left( 1 - F \right)\mu_{2} μ=μ1⋅F+(1−F)μ2
υ = μ / ρ \upsilon = \mu/\rho υ=μ/ρ
体积分数F与压力相同，定义在网格中心点，有限差分法，时间离散采用显示欧拉格式，对流项采用中心差分格式进行离散。
∂ F ∂ t = F n + 1 − F n Δ t \frac{\partial F}{\partial t} = \frac{F^{n + 1} - F^{n}}{\mathrm{\Delta}t} ∂t∂F=ΔtFn+1−Fn
∂ F ∂ x = F ( i + 1 , j ) − F ( i − 1 , j ) 2 Δ x \frac{\partial F}{\partial x} = \frac{F\left( i + 1,j \right) - F\left( i - 1,j \right)}{2\mathrm{\Delta}x} ∂x∂F=2ΔxF(i+1,j)−F(i−1,j)
∂ F ∂ y = F ( i , j + 1 ) − F ( i , j − 1 ) 2 Δ y \frac{\partial F}{\partial y} = \frac{F\left( i,j + 1 \right) - F\left( i,j - 1 \right)}{2\mathrm{\Delta}y} ∂y∂F=2ΔyF(i,j+1)−F(i,j−1)
加入VOF模型后的整体流程图如图6所示
图6 VOF求解流程
第二章 计算结果及求解程序
2.1 算例1 方腔流
计算Re=100的方腔流，方形区域边长为1m，速度为1m/s，动力粘度0.01，
图8 水平方向速度场（左）残差图（右）
计算300个时间步，残差值低于e-6，可以认为计算收敛达稳定值。
2.2 算例2 溃坝模型：
区域： 1 × 1 \ 1 \times 1 1×1
边界：壁面
网格点： 64 × 64 \ 64 \times 64 64×64
时间步：0.001
g x = 0 , g y = − 1 g_{x} = 0,\ \ g_{y} = - 1 gx=0, gy=−1
相1： ρ 1 = 1 , μ 1 = 1 × 1 0 − 3 \rho_{1} = 1\ ,\mu_{1} = 1 \times 10^{- 3} ρ1=1 ,μ1=1×10−3
相2： ρ 2 = 0.5 , μ 2 = 0.5 × 1 0 − 3 \rho_{2} = 0.5\ ,\mu_{2} = 0.5 \times 10^{- 3} ρ2=0.5 ,μ2=0.5×10−3
图9 0. 分别在0.1s，1s，2s，3s，4s，5s，6s，7s，8s时刻下体积分数分布图
第三章 Matlab求解程序使用说明
程序总共分为8个部分，Main，M-possion，cal_mu_rho，set_BC，solve_F，var，plot_Fra，post_plot，在main文件中输入模拟参数，M-possion，cal_mu_rho，set_BC，solve_F，var为求解过程中函数，plot_Fra为输出数据函数，post_plot实现后处理绘图功能。
3.1 输入配置参数
在主程序main中设定以下4类参数
1、网格参数
nx， ny：分别为x方向，y方向均分网格数目
Lx，Ly：分别为方形区域大小，左下角为原点
2、物理参数
两相流密度，运动粘度，water代表重流体属性，air为轻流体属性
Rho_water，rho_air，nu_water，nu_air
3、初始体积分数区域
以方形区域为例， x 1 ≤ x ≤ x 2 y 1 ≤ y ≤ y 2 x1 \leq x \leq x2\ \ \ \ \ \ \ \ \ \ \ \ \ y1 \leq y \leq y2 x1≤x≤x2 y1≤y≤y2
4、求解参数
时间步dt
3.2 输出结果
使用Plot_Fra来绘制体积分数分布图，可以设置输出数据类型，输出位置，也可自行编写其他数据输出文件代码。
附录
1. Main主体程序
%输入网格参数
global nx ny Lx Ly
global rho_water rho_air nu_water nu_air gx gy
global x1 y1 x2 y2
global dt
nx=32; %x方向网格点数
ny=32; %y方向网格点数
Lx=1; %区域长度
Ly=1; %区域高度
%输入物理参数
rho_water=1;
rho_air=0.5;
nu_water=0.01; %运动粘性系数
nu_air=0.005; %运动粘性系数
gx=0;
gy=-1;
%初始体积分数区域（方形区域）
x1=0;
x2=1/3;
y1=0;
y2=0.5;
%求解参数
dt=0.001;
%mesh=========================================================================
global imax imin jmax jmin dxi dyi
imin=2; imax=imin+nx-1;
jmin=2; jmax=jmin+ny-1;
x(imin : imax+1)=linspace (0 ,Lx, nx+1);
y(jmin : jmax+1)=linspace (0 ,Ly, ny+1);
xm(imin : imax)=0.5*(x(imin : imax)+x(imin+1:imax+1));
ym(jmin : jmax)=0.5*(y(jmin : jmax)+y(jmin+1:jmax+1));
dx=x(imin+1)-x(imin );
dy=y(jmin+1)-y(jmin );
dxi=1/dx ;
dyi=1/dy ;
%参数设置
%初始条件
u=zeros(imax+1,jmax+1);
v=zeros(imax+1,jmax+1);
p=zeros(imax,jmax);
F=zeros(imax+1,jmax+1);
%设置初始体积分数
for j=1:1:jmax
for i=1:1:imax
if(xm(i)>=x1)&&(xm(i)<=x2)&&(ym(j)>=y1)&&(ym(j)<=y2)
F(i,j)=1;
end
end
end
%初始化rho，mu
rho=zeros(imax+1,jmax+1);
mu=zeros(imax+1,jmax+1);
%=================================================
% 创建Laplace算子
L=Laplace_operator(nx,ny,dxi,dyi);
%=====================================================
istep=0;
istep_max=500;
Residual=zeros(istep_max/100,3);
check_mass=zeros(istep_max/100,1);%检查质量
R_limit=15;
count=0;
while (istep<istep_max)
%设置边界条件
[u,v,F]=set_BC(u,v,F);
%由F更新rho，mu
[rho,mu] = cal_mu_rho( F,rho,mu);
%投影法求解压力possion方程
[un,vn,pn] = M_Possion(L,u,v,p,mu,rho);
%显示，差分法求解F
[Fn] = solve_F( un,vn,F);
for j=1:1:jmax
for i=1:1:imax
Fn(i,j)=var(0,1,Fn(i,j));
end
end
istep=istep+1;%时间步+1
u_R=norm(un-u)/(nx*ny);
v_R=norm(vn-v)/(nx*ny);
p_R=norm(pn-p)/(nx*ny);
u=un;
v=vn;
p=pn;
F=Fn;
[u,v,Fn]=set_BC(u,v,Fn);
if (u_R>R_limit)||(v_R>R_limit)
fprintf('残差过大不收敛');
break;
end
if mod(istep,100)==0%每隔100步输出数据
count=count+1;
Residual(count,1)=u_R;
Residual(count,2)=v_R;
Residual(count,3)=p_R;
check_mass(count)=sum(sum(abs(Fn(imin:imax,jmin:jmax))));
fprintf('迭代次数%s\n u_R=%s\n v_R=%s\n check mass：%s\n',num2str(istep), ...
num2str(u_R),num2str(v_R),num2str(check_mass(count)));
%输出图片
plot_Fra( Fn,istep,xm,ym );
end
end
2. M-Possion
function [ u,v,p ] = M_Possion(L,u,v,p,mu,rho)
%投影法求解压力泊松方程
global gx gy
global imax imin jmax jmin dxi dyi dt
u_star=zeros(imax+1,jmax+1);
v_star=zeros(imax+1,jmax+1);
for j=jmin:jmax
for i=imin+1:imax
v_here =0.25*(v(i-1,j)+v(i-1,j+1)+v(i,j)+v(i,j +1));
u_star(i,j)=u(i,j)+dt*...
(mu(i,j)*(u(i-1,j)-2*u(i,j)+u(i+1,j))*dxi^2 ...
+mu(i,j)*(u(i,j-1)-2*u(i,j)+u(i,j+1))*dyi^2 ...
-u(i,j)*(u(i+1,j)-u(i-1,j))*0.5*dxi ...
-v_here*(u(i,j+1)-u(i,j-1))*0.5*dyi ...
+gx);
end
end
for j=jmin+1:jmax
for i=imin : imax
u_here =0.25*(u(i,j-1)+u(i,j)+u(i+1,j-1)+u(i+1,j));
v_star(i,j)=v(i,j)+dt* ...
(mu(i,j)*(v(i-1,j)-2*v(i,j)+v(i+1,j))*dxi^2 ...
+mu(i,j)*(v(i,j-1)-2*v(i,j)+v(i,j+1))*dyi^2 ...
-u_here*(v(i+1,j)-v(i-1,j))*0.5*dxi ...
-v(i,j)*(v(i,j+1)-v(i,j-1))*0.5*dyi ...
+gy);
end
end
n=0;
R=zeros(imax-imin+1,1);
for j=jmin : jmax
for i=imin : imax
n=n+1;
R(n)=-rho(i,j)/dt*((u_star(i+1,j)-u_star(i,j))*dxi ...
+(v_star(i,j+1)-v_star(i,j))*dyi);
end
end
pv=L\R;
n=0;
%p=zeros(imax,jmax);
for j=jmin : jmax
for i=imin : imax
n=n+1;
p(i,j)=pv(n);
end
end
for j=jmin : jmax
for i=imin+1:imax
u(i,j)=u_star(i,j)-dt/rho(i,j)*(p(i,j)-p(i-1,j))*dxi;
end
end
for j=jmin+1:jmax
for i=imin : imax
v(i,j)=v_star(i,j)-dt/rho(i,j)*(p(i,j)-p(i,j-1))*dyi;
end
end
end
4. set_BC
function [u,v,F] = set_BC(u,v,F)
global imax imin jmax jmin
%设置边界条件
%
%边界条件
% u_bottom=0;
% v_bottom=0;
% u_top=0;
% v_top=0;
% u_left=0;
% v_left=0;
% u_right=0;
% v_right=0;
%
% u(:,jmin-1)=-u(:,jmin)+2*u_bottom;
% v(:,jmin-1)=-v(:,jmin)+2*v_bottom;
%
% u(:,jmax+1)=-u(:,jmax)+2*u_top;
% v(:,jmax+1)=-v(:,jmax)+2*v_top;
%
% u(imin-1,:)=-u(imin,:)+2*u_left;
% v(imin-1,:)=-v(imin,:)+2*v_left;
%
% u(imax+1,:)=-u(imax,:)+2*u_right;
% v(imax+1,:)=-v(imax,:)+2*v_right;
%设置壁面为不可穿透条件，滑移
%bottom
%u(:,jmin-1)=-u(:,jmin)+2*u_bottom;
v(:,jmin-1)=v(:,jmin);
%top
%u(:,jmax+1)=-u(:,jmax)+2*u_top;
v(:,jmax+1)=v(:,jmax);
%left
u(imin-1,:)=u(imin,:);
%v(imin-1,:)=-v(imin,:)+2*v_left;
%right
u(imax+1,:)=u(imax,:);
%v(imax+1,:)=-v(imax,:)+2*v_right;
% fprintf('设置边界条件完成\n');
%为F设置不可穿透条件
F(imin-1,:)=F(imin,:);%left
F(imax+1,:)=F(imax,:);%right
F(:,jmin-1)=F(:,jmin);%bottom;
F(:,jmax+1)=F(:,jmax);%top
end
5. solve_F
function [Fn] = solve_F( u,v,F)
%差分法求解体积分数F
global imax imin jmax jmin dxi dyi
global dt
Fn=zeros(imax+1,jmax+1);
for j=jmin:1:jmax
for i=imin:1:imax
u_loc=3/8*(u(i,j)+u(i+1,j))+1/16* ...
(u(i,j+1)+u(i+1,j+1)+u(i,j-1)+u(i+1,j-1));
v_loc=3/8*(v(i,j)+v(i,j+1))+1/16* ...
(v(i-1,j)+v(i-1,j+1)+v(i+1,j)+v(i+1,j+1));
Fn(i,j)=F(i,j)-dt* ...
(u_loc*(F(i+1,j)-F(i-1,j))*dxi/2+ ...
v_loc*(F(i,j+1)-F(i,j-1))*dyi/2);
end
end
end
6. var
function [ center ] = var( a,b,c)
%返回三个数的中间值
center=a+b+c-(max([a,b,c])+min([a,b,c]));
end
7. plot_Fra
function [] = plot_Fra( F,istep,xm,ym)
%储存图片
global imax imin jmax jmin
global dt
%F
fig = figure('Visible','off'); % 新建一个figure，并将图像句柄保存到fig
pcolor(xm(imin:imax),ym(jmin:jmax),F(imin:imax,jmin:jmax)') %等高线图
%contourf(xm(imin:imax),ym(jmin:jmax),F(imin:imax,jmin:jmax)',1);
set(gcf,'position',[200,50,600,650]);
c = colorbar('southoutside','Ticks',[0 0.2 0.4 0.6 0.8 1]);
colormap('jet')
c.Label.String = 'Fraction of water';
xlabel('X','fontsize',10);
ylabel('Y','fontsize',10);
legd=['Time = ',num2str(istep*dt,'%02.1f\n')];
legend(legd);
frame = getframe(fig); % 获取frame
img = frame2im(frame); % 将frame变换成imwrite函数可以识别的格式
imwrite(img,['C:\\Users\\dell\\Documents\\MATLAB\\ACFD\\pic_fra\\','时间步数：',num2str(istep,'%06d\n'),'.png']); % 保存到指定目录下，名字为"*.png"
fig2 = figure('Visible','off'); % 新建一个figure，并将图像句柄保存到fig
contourf(xm(imin:imax),ym(jmin:jmax),F(imin:imax,jmin:jmax)',1);
set(gcf,'position',[200,50,600,650]);
c = colorbar('southoutside','Ticks',[0 0.2 0.4 0.6 0.8 1]);
colormap('jet')
c.Label.String = 'Fraction of water';
xlabel('X','fontsize',10);
ylabel('Y','fontsize',10);
legd=['Time = ',num2str(istep*dt,'%02.1f\n')];
legend(legd);
frame2 = getframe(fig2); % 获取frame
img = frame2im(frame2); % 将frame变换成imwrite函数可以识别的格式
imwrite(img,['C:\\Users\\dell\\Documents\\MATLAB\\ACFD\\pic_fra1\\','时间步数：',num2str(istep,'%06d\n'),'.png']); % 保存到指定目录下，名字为"*.png"
fprintf('储存图片:%s\n \n',num2str(istep));
end
8. post_plot
%速度x等高线图
figure(1);
contourf(x(imin:imax),y(jmin:jmax),u(imin:imax,jmin:jmax)',15) %等高线图
xlabel('X','fontsize',10);
ylabel('Y','fontsize',10);
axis equal
c = colorbar('westoutside');
c.Label.String = 'Vel of u';
title('Vel of u');
grid on
%压力图
figure(2);
contourf(xm(2:1+nx),ym(2:1+ny),p(2:1+nx,2:1+ny)',5) %等高线图
colorbar
title('Pressure');
grid on
%残差图
figure(3)
plot(100:100:istep,Residual(:,1),'r-','LineWidth',2);
hold on
plot(100:100:istep,Residual(:,2),'g-','LineWidth',2);
plot(100:100:istep,Residual(:,3),'b-','LineWidth',2);
set(gca,'yscale','log');
xlabel('istep','fontsize',10);
ylabel('Residual','fontsize',10);
grid on
legend('u_R','v_R','p_R');
title('Residual-istep');
%F
figure(4)
contourf(x(imin:imax),y(jmin:jmax),Fn(imin:imax,jmin:jmax)',1) %等高线图
9. Laplace_operator
function [L] = Laplace_operator(nx,ny,dxi,dyi)
%计算Laplace算子矩阵
L=zeros(nx*ny,nx*ny);
for j =1:1:ny
for i =1:1:nx
L(i+(j-1)*nx,i+(j-1)*nx)=2*dxi^2+2*dyi^2;
for ii=i-1:2:i+1
if (ii>0 && ii<=nx) % Interior point
L(i+(j-1)*nx,ii+(j-1)*nx)=-dxi^2;
else % Neuman conditions on boundary
L(i+(j-1)*nx,i+(j-1)*nx)= ...
L(i+(j-1)*nx,i+(j-1)*nx)-dxi^2;
end
end
for jj=j-1:2: j+1
if (jj>0 && jj<=ny) % Interior point
L(i+(j-1)*nx,i+(jj-1)*nx)=-dyi^2;
else % Neuman conditions on boundary
L(i+(j-1)*nx,i+(j-1)*nx)= ...
L(i+(j-1)*nx,i+(j-1)*nx)-dyi^2;
end
end
end
end
%Set pressure in first cell (all other pressures w.r.t to this one)
L(1 ,:)=0; L(1 ,1)=1;
end
10. cal_mu_rho
function [ rho,mu ] = cal_mu_rho( F,rho,mu,imin,imax,jmin,jmax)
%计算密度和mu
rho_water=1;
rho_air=0.5;
nu_water=1e-3;%运动粘性系数
nu_air=0.5e-3;
for j=jmin-1:1:jmax+1
for i=imin-1:1:imax+1
rho(i,j)=rho_air*(1-var(0,1,F(i,j)))+rho_water*var(0,1,F(i,j));
mu(i,j)=(nu_water*(1-var(0,1,F(i,j)))+nu_air*var(0,1,F(i,j)))/ rho(i,j);
end
end
end
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
