驱动方腔流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}
ΔtU∗−Un=−((Un)2)x−(UnVn)yΔtV∗−Vn=−(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,j−Ui,j=2Ui,j+1−Ui,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}
ΔtU∗−UΔtV∗−V=−((Uh)2−γ∣∣∣Uh∣∣∣U~h)x−((UvVh)−γ∣∣∣Vh∣∣∣U~v)y=−((UvVh)−γ∣∣∣Uv∣∣∣V~h)x−((Vv)2−γ∣∣∣Vv∣∣∣V~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⋅Δt⋅max(i,jmax∣Ui,j∣,i,jmax∣Vi,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}
ΔtU∗∗−U∗ΔtV∗∗−V∗=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=−Δt1∇⋅Un
得到压力
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+1−U∗∗ΔtVn+1−V∗∗=−(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}
ΔtU∗−UnΔtV∗−Vn=−UUx−VUy=−UVx−VVy
同样地，我们需要做一些平均，使得等式左右两边对应的点照上。以第一个式子为例：
对
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)⊥=u⟺−qy=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=−qyy−qxx=uy−vx
所以说，求流函数很简单，同样是交错网格，我们把流函数的值定义在网格交点上。那么，要得到流函数，只要求解泊松方程：
−
Δ
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;
扫一扫
3924
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
