驱动方腔流SIMPLE方法

驱动方腔流SIMPLE方法

问题描述:二维驱动方腔流

请在矩形区域上数值地求解不可压NS方程,你的报告将包含但不仅限于以下的讨论:

  • 问题描述和数值格式

  • 格式的时空精度分析

  • 比较非线性项处理的守恒格式和非守恒格式的不同

  • 不同雷诺数的影响(分别取100、1000、3000)

数值格式:SIMPLE格式

问题的数学表达

考虑二维空间的不可压NS方程:
u t + p x = − ( u 2 ) x − ( u v ) y + 1 R e ( u x x + u y y ) v t + p y = − ( u v ) x − ( v 2 ) y + 1 R e ( v x x + v y y ) u x + v y = 0 \begin{aligned} u_{t}+p_{x} &=-\left(u^{2}\right)_{x}-(u v)_{y}+\frac{1}{R e}\left(u_{x x}+u_{y y}\right) \\ v_{t}+p_{y} &=-(u v)_{x}-\left(v^{2}\right)_{y}+\frac{1}{R e}\left(v_{x x}+v_{y y}\right) \\ u_{x}+v_{y} &=0 \end{aligned} ut+pxvt+pyux+vy=(u2)x(uv)y+Re1(uxx+uyy)=(uv)x(v2)y+Re1(vxx+vyy)=0
其中,求解区域为: Ω = [ 0 , l x ] × [ 0 , l y ] \Omega=\left[0, l_{x}\right] \times\left[0, l_{y}\right] Ω=[0,lx]×[0,ly],四条边界按东南西北依次记为 E S W N \mathbf{ESWN} ESWN。边界条件考虑无滑移的边界条件,即:
u ( x , l y ) = u N ( x ) v ( x , l y ) = 0 u ( x , 0 ) = u S ( x ) v ( x , 0 ) = 0 u ( 0 , y ) = 0 v ( 0 , y ) = v W ( y ) u ( l x , y ) = 0 v ( l x , y ) = v E ( y ) \begin{array}{rl}{u\left(x, l_{y}\right)=u_{N}(x)} & {v\left(x, l_{y}\right)=0} \\ {u(x, 0)=u_{S}(x)} & {v(x, 0)=0} \\ {u(0, y)=0} & {v(0, y)=v_{W}(y)} \\ {u\left(l_{x}, y\right)=0} & {v\left(l_{x}, y\right)=v_{E}(y)}\end{array} u(x,ly)=uN(x)u(x,0)=uS(x)u(0,y)=0u(lx,y)=0v(x,ly)=0v(x,0)=0v(0,y)=vW(y)v(lx,y)=vE(y)
结合不可压条件,方程中的非线性项写为非守恒的模样,为:
( u 2 ) x + ( u v ) y = u u x + v u y ( u v ) x + ( v 2 ) y = u v x + v v y \begin{array}{l}{\left(u^{2}\right)_{x}+(u v)_{y}=u u_{x}+v u_{y}} \\ {(u v)_{x}+\left(v^{2}\right)_{y}=u v_{x}+v v_{y}}\end{array} (u2)x+(uv)y=uux+vuy(uv)x+(v2)y=uvx+vvy
这其实就是NS方程中的对流项。

交错网格

在这里插入图片描述

和一般的差分说不同的是,交错网格将不同的速度分量,以及压力标记在不同的位置。如图所示,实心圆处表示速度分量 U U U的位置,空心圆处表示速度分量 V V V的位置,网格的中点用叉叉表示的地方是压力 P P P所在的位置。那么可以看到,在考虑边界条件的时候, U U U在左右两边刚好达到边界,而上下两边则需要在边界处还要多往外扩充半格。对于 V V V也是一样的道理。而压力 P P P的边界,在四个面上都往外半格。这些扩充出去的半个格点处的值要怎么求呢?用最为普通的差分处理方法,结合边界条件即可,十分简单。这就是所谓的"交错网格"。

求解步骤

不再琐碎地叙述非线性项怎么处理,粘性项怎么处理,压力怎么矫正等等,下面给出一个顺序的求解步骤描述,希望读者看着这些步骤,结合交错网格的编号,很容易就写出程序。

非线性项处理

一般问题,方程中的非线性项只能显式处理,差分如是,有限元也如是(非线性项不好做变分),但是显式处理的问题就是有些时候不太稳定,在这里表现的就是对CFL数有要求。
U ∗ − U n Δ t = − ( ( U n ) 2 ) x − ( U n V n ) y V ∗ − V n Δ t = − ( U n V n ) x − ( ( V n ) 2 ) y \begin{array}{l}{\frac{U^{*}-U^{n}}{\Delta t}=-\left(\left(U^{n}\right)^{2}\right)_{x}-\left(U^{n} V^{n}\right)_{y}} \\ {\frac{V^{*}-V^{n}}{\Delta t}=-\left(U^{n} V^{n}\right)_{x}-\left(\left(V^{n}\right)^{2}\right)_{y}}\end{array} ΔtUUn=((Un)2)x(UnVn)yΔtVVn=(UnVn)x((Vn)2)y
通过这个,我们就可以算一步到 U ∗ , V ∗ U^*,V^* U,V。问题在于,等式右边的非线性项要怎么离散(这里举的是守恒的表示,非守恒也一样)。这里的困难在于,我们用的是交错格式,离散变量都没取在相同的位置,况且在做个导数离散,导数值都指不定到哪些点上去。怎么办呢?我们可以通过已有点的值,通过平均的方法,生成更多点值,用这些生成的点的值,就可以构造正确位置的导数离散。
因为中心差分不稳定,所以非线性项的离散我们一般用中心差分和迎风结合的方法,什么意思呢?
我们分别计算 U U U的水平平均(到了每个格子中心),以及 U U U的垂直平均(到了网格交点上),即:
( U ˉ h ) i + 1 2 , j = U i , j + U i + 1 , j 2 ( U ˉ v ) i , j + 1 2 = U i , j + U i , j + 1 2 \begin{aligned}\left(\bar{U}^{h}\right)_{i+\frac{1}{2}, j} &=\frac{U_{i, j}+U_{i+1, j}}{2} \\\left(\bar{U}^{v}\right)_{i, j+\frac{1}{2}} &=\frac{U_{i, j}+U_{i, j+1}}{2} \end{aligned} (Uˉh)i+21,j(Uˉv)i,j+21=2Ui,j+Ui+1,j=2Ui,j+Ui,j+1
类似,定义网格差分量为:
( U ~ h ) i + 1 2 , j = U i + 1 , j − U i , j 2 ( U ~ v ) i , j + 1 2 = U i , j + 1 − U i , j 2 \begin{aligned}\left(\tilde{U}^{h}\right)_{i+\frac{1}{2}, j} &=\frac{U_{i+1, j}-U_{i, j}}{2} \\\left(\tilde{U}^{v}\right)_{i, j+\frac{1}{2}} &=\frac{U_{i, j+1}-U_{i, j}}{2} \end{aligned} (U~h)i+21,j(U~v)i,j+21=2Ui+1,jUi,j=2Ui,j+1Ui,j
V V V也做类似的定义,就不再敲表达式了。不同的是,对 V V V来说,做垂直的平均和差分,值在格子中心,而水平的,到了格点(线交点)上。那么,计算格式为:
U ∗ − U Δ t = − ( ( U ‾ h ) 2 − γ ∣ U ‾ h ∣ U ~ h ) x − ( ( U ‾ v V ‾ h ) − γ ∣ V ‾ h ∣ U ~ v ) y V ∗ − V Δ t = − ( ( U ‾ v V ‾ h ) − γ ∣ U ‾ v ∣ V ~ h ) x − ( ( V ‾ v ) 2 − γ ∣ V ‾ v ∣ V ~ v ) y \begin{aligned} \frac{U^{*}-U}{\Delta t} &=-\left(\left(\overline{U}^{h}\right)^{2}-\gamma\left|\overline{U}^{h}\right| \tilde{U}^{h}\right)_{x}-\left(\left(\overline{U}^{v} \overline{V}^{h}\right)-\gamma\left|\overline{V}^{h}\right| \tilde{U}^{v}\right)_{y} \\ \frac{V^{*}-V}{\Delta t} &=-\left(\left(\overline{U}^{v} \overline{V}^{h}\right)-\gamma\left|\overline{U}^{v}\right| \tilde{V}^{h}\right)_{x}-\left(\left(\overline{V}^{v}\right)^{2}-\gamma\left|\overline{V}^{v}\right| \tilde{V}^{v}\right)_{y} \end{aligned} ΔtUUΔtVV=((Uh)2γUhU~h)x((UvVh)γVhU~v)y=((UvVh)γUvV~h)x((Vv)2γVvV~v)y
这样,不难看到,对右边做中心离散(单格),位置就和等式左边的照上了,还是 U , V U,V U,V原来所在的位置。这里的 γ \gamma γ取值如下:
γ = min ⁡ ( 1.2 ⋅ Δ t ⋅ max ⁡ ( max ⁡ i , j ∣ U i , j ∣ , max ⁡ i , j ∣ V i , j ∣ ) , 1 ) \gamma=\min \left(1.2 \cdot \Delta t \cdot \max \left(\max _{i, j}\left|U_{i, j}\right|, \max _{i, j}\left|V_{i, j}\right|\right), 1\right) γ=min(1.2Δtmax(i,jmaxUi,j,i,jmaxVi,j),1)

粘性项处理

隐式处理粘性项:
U ∗ ∗ − U ∗ Δ t = 1 R e ( U x x ∗ ∗ + U y y ∗ ∗ ) V ∗ ∗ − V ∗ Δ t = 1 R e ( V x x ∗ ∗ + V y y ∗ ∗ ) \begin{aligned} \frac{U^{* *}-U^{*}}{\Delta t} &=\frac{1}{R e}\left(U_{x x}^{* *}+U_{y y}^{* *}\right) \\ \frac{V^{* *}-V^{*}}{\Delta t} &=\frac{1}{R e}\left(V_{x x}^{* *}+V_{y y}^{* *}\right) \end{aligned} ΔtUUΔtVV=Re1(Uxx+Uyy)=Re1(Vxx+Vyy)
二阶导数的离散直接采用最一般的二阶中心离散方法即可,利用已定义的位置上的值。

压力校正

对于顶盖驱动方腔问题,压力用的是齐次的诺依曼边界条件。校正的步骤如下:

  • 计算离散的散度值: F n = ∇ ⋅ U n F^{n}=\nabla \cdot \mathbf{U}^{n} Fn=Un。容易想到,这个值在格中心。

  • 求解泊松方程: − Δ P n + 1 = − 1 Δ t F n -\Delta P^{n+1}=-\frac{1}{\Delta t} F^{n} ΔPn+1=Δt1Fn。因为压力的定义,也是在格子中心,这个求解是自然而然的。

  • 计算压力梯度: G n + 1 = ∇ P n + 1 \mathbf{G}^{n+1}=\nabla P^{n+1} Gn+1=Pn+1。显然,这个梯度包含了两个方面,对 x x x
    的导数(刚好落到 U U U的位置),和对 y y y的导数(刚好落到 V V V的位置)。

  • 更新速度域: U n + 1 = U n − Δ t G n + 1 \mathbf{U}^{n+1}=\mathbf{U}^{n}-\Delta t \mathbf{G}^{n+1} Un+1=UnΔtGn+1

当然,这里的 U U U取的实际是上一步得到的 U ∗ ∗ U^{**} U。用公式来表达这个过程,就是:
求解:
− Δ P n + 1 = − 1 Δ t ∇ ⋅ U n -\Delta P^{n+1}=-\frac{1}{\Delta t} \nabla \cdot \mathbf{U}^{n} ΔPn+1=Δt1Un
得到压力 P P P后,求解:
U n + 1 − U ∗ ∗ Δ t = − ( P n + 1 ) x V n + 1 − V ∗ ∗ Δ t = − ( P n + 1 ) y \begin{aligned} \frac{U^{n+1}-U^{* *}}{\Delta t} &=-\left(P^{n+1}\right)_{x} \\ \frac{V^{n+1}-V^{* *}}{\Delta t} &=-\left(P^{n+1}\right)_{y} \end{aligned} ΔtUn+1UΔtVn+1V=(Pn+1)x=(Pn+1)y

精度分析

往观这些求解格式,空间上离散采用的基本上算是中心差分,时间上采用的欧拉离散(显示和隐式),所以目测这个格式的精度应该是 O ( Δ x 2 + Δ t ) O(\Delta x ^2+\Delta t) O(Δx2+Δt)
另外一点,我们再取 γ \gamma γ的时候,其表达式前面的1.2取得不好,对精度有影响。
时间和空间的精度分析,做泰勒展开可以很容易得到,为了节省时间,报告也马上就要提交了,这里就不再细述了。

非线性项的守恒格式和非守恒格式

以上提到的都是守恒格式非线性项的处理,对于非守恒格式,可以这个样子处理非线性项:
U ∗ − U n Δ t = − U U x − V U y V ∗ − V n Δ t = − U V x − V V y \begin{aligned} \frac{U^{*}-U^{n}}{\Delta t} &=-UU_{x}-VU_{y} \\ \frac{V^{*}-V^{n}}{\Delta t} &=-UV_{x}-VV_{y} \end{aligned} ΔtUUnΔtVVn=UUxVUy=UVxVVy
同样地,我们需要做一些平均,使得等式左右两边对应的点照上。以第一个式子为例:
U x U_x Ux的离散,我们首先得对 U U U x x x方向上的平均,把值挪到网格中心,得到 U ˉ h \bar U^h Uˉh。对于 V V V,我们可以做 U U U周围的四点平均,把值挪到 U U U所在的位置。对于 U y U_y Uy的离散,我们可以做 U U U的竖直平均得到网格交点上的值 U ˉ v \bar U^v Uˉv后,再做中心差分近似,得到的值刚好也在 U U U的位置。
思路就是这样,程序修改也很简单。考虑到非守恒总是不如守恒的,以及时间问题,就不再详细去做。

不同雷诺数的结果比较

可视化

为了画图看效果,我们可以画一下流线,画流线,需要用到流函数,取流函数的等高线,就是流场。所谓的流线,其实就是各个点顺着速度走的方向,串在一块。那么流函数就要在垂直速度的方向上进行递增,这样,它的水平集才能体现流线。
q q q来表示流函数,即有:
( ∇ q ) ⊥ = u ⟺ − q y = u  and  q x = v (\nabla q)^{\perp}=\mathbf{u} \quad \Longleftrightarrow \quad-q_{y}=u \text { and } q_{x}=v (q)=uqy=u and qx=v
两边儿都加个旋度算子,就变成了:
− Δ q = ∇ × ( ∇ q ) ⊥ = ∇ × u ⟺ − Δ q = − q y y − q x x = u y − v x -\Delta q=\nabla \times(\nabla q)^{\perp}=\nabla \times \mathbf{u} \quad \Longleftrightarrow \quad-\Delta q=-q_{y y}-q_{x x}=u_{y}-v_{x} Δq=×(q)=×uΔq=qyyqxx=uyvx
所以说,求流函数很简单,同样是交错网格,我们把流函数的值定义在网格交点上。那么,要得到流函数,只要求解泊松方程:
− Δ Q n = ( U n ) y − ( V n ) x -\Delta Q^n = (U^n)_y-(V^n)_x ΔQn=(Un)y(Vn)x

不同雷诺数下的观察

使用MIT的程序,跑了跑,贴几张图。MIT的程序可以在这里下载到。
下面我修改程序,把雷诺数设为100,1000,3000,来看看情况。对于标量函数 P P P来说,用一张热图(颜色图)就可以刻画趋势,对于矢量,除了大小,还有方向,要怎么呈现呢?方向可以用带箭头的图,以及流线来体现,大小还是可以用色图,对吧。我把流线图和速度矢量图画在一张图上,压力图画一张图。

在这里插入图片描述

从图中可以看出,在雷诺数为100时,在角落没有形成涡。

在这里插入图片描述

从图中可以看出,在雷诺数为1000时,右下角形成了一个涡。

在这里插入图片描述

从图中可以看出,在雷诺数为1000时,左下角和右下角各形成了一个涡。

简单修改后的代码和注释

function  mit18086_navierstokes
%NAVIERSTOKES
%    Solves the incompressible Navier-Stokes equations in a
%    rectangular domain with prescribed velocities along the
%    boundary. The solution method is finite differencing on
%    a staggered grid with implicit diffusion and a Chorin
%    projection method for the pressure.
%    Visualization is done by a colormap-isoline plot for
%    pressure and normalized quiver and streamline plot for
%    the velocity field.
%    The standard setup solves a lid driven cavity problem.
%-----------------------------------------------------------------------
clc
clear
close all
Re = 3000;     % Reynolds number100 1000 3000  %%%%%%%%%%%%%%%%%%%
dt = 1e-2;    % time step
tf = 50;    % final time 10 30 50%%%%%%%%%%%%%%%%%%%%%%%
lx = 1;       % width of box
ly = 1;       % height of box
nx = 100;      % number of x-gridpoints
ny = 100;      % number of y-gridpoints
nsteps = 10;%round(tf/dt/10);  % number of steps with graphic output
%-----------------------------------------------------------------------
nt = ceil(tf/dt); dt = tf/nt;
x = linspace(0,lx,nx+1); hx = lx/nx;
y = linspace(0,ly,ny+1); hy = ly/ny;
[X,Y] = meshgrid(y,x);
%-----------------------------------------------------------------------
% initial conditions
U = zeros(nx-1,ny); V = zeros(nx,ny-1);
% boundary conditions
uN = x*0+1;    vN = avg(x)*0;
uS = x*0;      vS = avg(x)*0;
uW = avg(y)*0; vW = y*0;
uE = avg(y)*0; vE = y*0;
%-----------------------------------------------------------------------
Ubc = dt/Re*([2*uS(2:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hx^2+...
    [uW;zeros(nx-3,ny);uE]/hy^2);
Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hx^2+...
    [2*vW(2:end-1);zeros(nx-2,ny-1);2*vE(2:end-1)]/hy^2);

fprintf('initialization')
%拉普拉斯算子的矩阵逼近
Lp = kron(speye(ny),K1(nx,hx,1))+kron(K1(ny,hy,1),speye(nx));
Lp(1,1) = 3/2*Lp(1,1);
%symamd,重排求解的线性方程组,使高斯消元的时候,即可能少了产生非零元
perp = symamd(Lp); Rp = chol(Lp(perp,perp)); Rpt = Rp';
Lu = speye((nx-1)*ny)+dt/Re*(kron(speye(ny),K1(nx-1,hx,2))+...
    kron(K1(ny,hy,3),speye(nx-1)));
peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut = Ru';
Lv = speye(nx*(ny-1))+dt/Re*(kron(speye(ny-1),K1(nx,hx,3))+...
    kron(K1(ny-1,hy,2),speye(nx)));
perv = symamd(Lv); Rv = chol(Lv(perv,perv)); Rvt = Rv';
Lq = kron(speye(ny-1),K1(nx-1,hx,2))+kron(K1(ny-1,hy,2),speye(nx-1));
perq = symamd(Lq); Rq = chol(Lq(perq,perq)); Rqt = Rq';

fprintf(', time loop\n--20%%--40%%--60%%--80%%-100%%\n')
for k = 1:nt
    % treat nonlinear terms
    gamma = min(1.2*dt*max(max(max(abs(U)))/hx,max(max(abs(V)))/hy),1);
    %首先扩充速度矩阵,把矩阵倒着看,即行标对应x,列标对应y
    Ue = [uW;U;uE]; Ue = [2*uS'-Ue(:,1) Ue 2*uN'-Ue(:,end)];
    Ve = [vS' V vN']; Ve = [2*vW-Ve(1,:);Ve;2*vE-Ve(end,:)];
    Ua = avg(Ue')'; Ud = diff(Ue')'/2;
    Va = avg(Ve);   Vd = diff(Ve)/2;
    %计算非线性项,UVx和UVy
    UVx = diff(Ua.*Va-gamma*abs(Ua).*Vd)/hx;
    UVy = diff((Ua.*Va-gamma*Ud.*abs(Va))')'/hy;
    %计算U方x和V方y
    Ua = avg(Ue(:,2:end-1));   Ud = diff(Ue(:,2:end-1))/2;
    Va = avg(Ve(2:end-1,:)')'; Vd = diff(Ve(2:end-1,:)')'/2;
    U2x = diff(Ua.^2-gamma*abs(Ua).*Ud)/hx;
    V2y = diff((Va.^2-gamma*abs(Va).*Vd)')'/hy;
    %更新内部点的值
    U = U-dt*(UVy(2:end-1,:)+U2x);
    V = V-dt*(UVx(:,2:end-1)+V2y);

    % implicit viscosity
    %把U拉成一列来搞
    rhs = reshape(U+Ubc,[],1);
    u(peru) = Ru\(Rut\rhs(peru));
    U = reshape(u,nx-1,ny);
    rhs = reshape(V+Vbc,[],1);
    v(perv) = Rv\(Rvt\rhs(perv));
    V = reshape(v,nx,ny-1);

    % pressure correction
    %压力校正,计算速度域的散度
    %计算右端项的矩阵形式,reshape,求解线性系统,在reshape回来,更新
    %内部点的值
    rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1);
    p(perp) = -Rp\(Rpt\rhs(perp));
    P = reshape(p,nx,ny);
    %通过-\nabla P更新内部点的值
    U = U-diff(P)/hx;
    V = V-diff(P')'/hy;

    % visualization
    if floor(25*k/nt)>floor(25*(k-1)/nt), fprintf('.'), end
    if k==1|floor(nsteps*k/nt)>floor(nsteps*(k-1)/nt)
        % stream function
        % 计算速度域的旋度
        %流函数,解线性方程组,为了可视化,齐次狄利克雷边界点包含在矩阵Q中
        rhs = reshape(diff(U')'/hy-diff(V)/hx,[],1);
        q(perq) = Rq\(Rqt\rhs(perq));
        Q = zeros(nx+1,ny+1);
        Q(2:end-1,2:end-1) = reshape(q,nx-1,ny-1);
        % clf, contourf(avg(x),avg(y),P',20,'w-'), hold on
        Ue = [uS' avg([uW;U;uE]')' uN'];
        Ve = [vW;avg([vS' V vN']);vE];
        Len = sqrt(Ue.^2+Ve.^2+eps);

        figure(1);
        subplot(1,2,2);
        contourf(avg(x),avg(y),P',20,'w-');%压力图
        axis equal, axis([0 lx 0 ly])
        %    axis('')
        p = sort(p);
        caxis(p([8 end-7]));
        title('pressure P');

        subplot(1,2,1);
        contour(x,y,Q',40,'k-');%流线图 300 100 40%%%%%%%%%%%%%%%%%%%%%%%%%%
        axis equal, axis([0 lx 0 ly])
        hold on;
        %  p = sort(p);
        %  caxis(p([8 end-7]))

        %      subplot(2,2,2);
        stride = 2;
        x_temp = x(1:stride:end);
        y_temp = y(1:stride:end);
        Ue_temp = Ue(1:stride:end,1:stride:end);
        Ve_temp = Ve(1:stride:end,1:stride:end);
        Len_temp = Len(1:stride:end,1:stride:end);

        quiver(x_temp,y_temp,(Ue_temp./Len_temp)',(Ve_temp./Len_temp)',.4,'g.')%速度的箭头矢量图
        axis equal, axis([0 lx 0 ly])
        % p = sort(p);
        % caxis(p([8 end-7]))
        title('velocity vector U & stream function Q');

        %   hold off,
        %        axis equal, axis([0 lx 0 ly])
        %        p = sort(p);
        %        caxis(p([8 end-7]))
        %  caxis('auto')
        suptitle(sprintf('Re = %0.1g   t = %0.2g',Re,k*dt));
        drawnow
    end
end
fprintf('\n')

%=======================================================================

function B = avg(A,k)
if nargin<2, k = 1; end
if size(A,1)==1, A = A'; end
if k<2, B = (A(2:end,:)+A(1:end-1,:))/2; else, B = avg(A,k-1); end
if size(A,2)==1, B = B'; end

function A = K1(n,h,a11)
% a11: Neumann=1, Dirichlet=2, Dirichlet mid=3;
A = spdiags([-1 a11 0;ones(n-2,1)*[-1 2 -1];0 a11 -1],-1:1,n,n)'/h^2;
已标记关键词 清除标记
©️2020 CSDN 皮肤主题: 点我我会动 设计师:上身试试 返回首页