计算流体力学续_全隐式格式_追赶法求解
计算流体力学—张涵信[(P395-P423)+(5.1)]
![计算流体力学—张涵信[(P395-P423)+(5.1)]](https://img.taocdn.com/s3/m/cda3861210a6f524ccbf8551.png)
5.1 研究高阶精度差分格式的必要性在任意贴体曲线坐标系中，支配可压缩黏性非定常流动的无量纲形式的NaVier-Stokes 方程为)ˆˆˆ(1ˆˆˆˆςηξςηξτυυυ∂∂+∂∂+∂∂=∂∂+∂∂+∂∂+∂∂G F E G F E Q L Re （5.1.1） 式中⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=τρωρυρρE u J Q 1ˆ；⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-++++=p U p E p U p U p uU U z y x ττξξρωξρυξρρ)(1ˆJ E ； ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-++++=p V p E p V p V p uV V z y x ττηηρωηρυηρρ)(1ˆJ F ；⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-++++=p W p E p W p W p uW W z y x ττςςρωςρυςρρ)(1ˆJ G ； ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡++++++++=)(Re )(Re )(Re )(Re 01ˆz z y y x x L zz z yz y xz x L yz z yy y xy x L xz z xy y xx x L K K K ξξξτξτξτξτξτξτξτξτξτξυJ E ； ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡++++++++=)(Re )(Re )(Re )(Re 01ˆz z y y x x L zz z yz y xz x L yz z yy y xy x L xz z xy y xx x L K K K ηηητητητητητητητητητηυJ F ； ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡++++++++=)(Re )(Re )(Re )(Re 01ˆz z y y x x L zz z yz y xz x L yz z yy y xy x L xz z xy y xx x L K K K ςςςτςτςτςτςτςτςτςτςτςυJ G ；x xz xy xx x q u K -++=ωτυττ ；y yz yy xy y q u K -++=ωτυττ ；z zz yz xz z q u K -++=ωτυττ ；)](21)1([)](21[222222ωυργρωυρτ+++-=+++=u p u e E ； W V U ，，分别为沿ςηξ，，坐标的逆变速度，它们的表达式为ωςυςςςωηυηηηωξυξξξτττz y x z y x z y x u W u V u U +++=+++=+++=坐标转换关系为)()()(z y x t z y x t z y x t t，，，，，，，，，ςςηηξξτ==== J 为坐标转换Jacobian 行列式⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=z y x z y x z y x ςςςηηηξξξJ应该指出，υυυG F E ˆˆˆ和，中包含黏性应力张量诸分量和导热向量的诸分量，它们所含有的诸导数也应作相应的坐标变换。
计算程序_计算流体力学_对流方程_有限差分法_Lax格式_迎风格式_FTCS格式
% 一维对流方程迎风格式、Lax格式、FTCS格式差分法计算% 潭花林清华大学航天航空学院% FTCS格式对于一维对流方程不稳定，最好不用clcclear all% 1.参数定义dx=1;x1=-18;x2=18;x=x1:dx:x2;L1=length(x);% dt=0.5*dx; % 收敛dt=2*dx; % 不收敛t1=0;t2=t1+80*dt;t=t1:dt:t2;L2=length(t);alpha=1;lambda=alpha*dt/dx;geshi=1; % 迎风格式% geshi=2; % Lax格式% geshi=3; % FTCS格式% 2.显式求解zeta=zeros(L1,L2);for kk=1:3geshi=kk;for ii=1:L1if x(ii)>0zeta(ii,1)=1;elseif x(ii)==0zeta(ii,1)=1/2;elseif x(ii)<0zeta(ii,1)=0;endendendendif geshi==1for ii=2:L1for jj=1:(L2-1)zeta(ii,jj+1)=zeta(ii,jj)-lambda*(zeta(ii,jj)-z eta(ii-1,jj));endzeta(1,jj+1)=zeta(2,jj+1);endzeta1=zeta;else if geshi==2for ii=2:(L1-1)for jj=1:(L2-1)zeta(ii,jj+1)=(zeta(ii+1,jj)+zeta(ii-1,jj))/2-. ..lambda/2*(zeta(ii+1,jj)-zeta(ii-1,jj));endzeta(1,jj+1)=zeta(2,jj+1);zeta(L1,jj+1)=zeta(L1,jj)-lambda*(zeta(L1,jj)-z eta(L1-1,jj));endzeta2=zeta;else if geshi==3for ii=2:(L1-1)for jj=1:(L2-1)zeta(ii,jj+1)=zeta(ii,jj)-lambda/2*(zeta(ii+1,j j)-zeta(ii-1,jj));endzeta(1,jj+1)=zeta(2,jj+1);zeta(L1,jj+1)=zeta(L1,jj)-lambda*(zeta(L1,jj)-z eta(L1-1,jj));endzeta3=zeta;endendendend% 3.绘图% 3.1 t=0figure(1)n=1;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=0时刻的计算结果') %%% 标题axis([-18,18,-0.2,1.2])% 3.2 t=10figure(2)n=(10-t(1))/dt;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=10s时刻的计算结果') %%% 标题% 3.3 t=20figure(3)n=(20-t(1))/dt;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') %%% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=20s时刻的计算结果') %%% 标题% 3.4 t=40figure(4)n=(40-t(1))/dt;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=40s时刻的计算结果') %%% 标题。
计算流体力学数值方法
3-43
计算流体力学
高阶精度可通过采用更多的节点值近似来获 得。一个节点允许的最高精度为1阶，两个节点允 许的最高精度为2阶，依此类推。 理论上讲，某种数值方法的精度越高，随着 网格的加密，误差减小的就越大 。也就是说，采 用高精度的数值方法，只需较少的网格数即可获 得要求的精度。 但是，高阶精度的方法常常需要更多的计算 时间，而且常常会导致解的有界性问题。
计算流体力学
解析解：
d dT (kA ) 0 dx dx d 2T 0 2 dx T c1 x c0 T ( x 0) 100 T ( x 1) 500 T 400x 100
3-21
计算流体力学
控制方程扩散项的离散 梯度扩散项的离散几乎 总是采用中心差分格式：
3-44
计算流体力学
3-5
计算流体力学
d （V） (C A ) SV dt n faces
C un A
1，有限体积法直接对上式进行离散 2，本章只考虑稳态问题，即上式左边第一项为 零
3-6
计算流体力学
有限体积法(FVM) (1) 定义流场求解域几何形 状 (2) 将求解域划分为计算网 格，即一组互不重叠的有限 体或单元。 (3) 基于上述划分的单元对 积分方程进行离散，即用节 点值来近似。 (4) 对得到的离散方程进行 数值求解。
3-11
计算流体力学
对于如图所示的一维控制体，物理量的守恒 可表述为如下关系式： [通量]e (fluxe)- [通量]w (fluxw) ＝源(source) 这里的通量是指穿过单 元表面的输运率。 如果 表示单位质量 的输运量，则总的通量为 对流通量和扩散通量之和， 其中： 对流通量＝ 扩散通量＝
计算流体力学CFD课件
随流体运动的有限控制体模型
连续性方程
质量守恒定律
有限控制体的总质量为：
m dV V
随流体运动的有限控制 体模型
随流体运动的有限控制体模型
连续性方程：
D Dt
V
dV
0
随流体运动的有限控制 体模型
空间位置固定的无穷小微团模型
空间位置固定的无穷小微团模型
连续性方程
质量守恒定律
流出微团的质量流量 ＝微团内质量的减少
动量方程
表面力的两个 来源： 1）压力 2）粘性力
动量方程
粘性力的两个 来源：
1）正应力 2）切应力
动量方程
切应力：与流体剪切变形的时间变化率有关， 如下图中的xy
动量方程
正应力：与流体微团体积的时间变化率有关， 如下图中的xx
动量方程
作用在单位质量流体微团 上的体积力记做 f ，其X
方向的分量为 f x
随流体运动的有限控制 体，同一批流体质点始 终位于同一控制体内
速度散度及其物理意义
速度散度的物理意义：
是每单位体积运动着
的流体微团，体积相对变化的时间变化率。
连续性方程
空间位置固定的有限控制体模型
空间位置固定的有限控制体模型
连续性方程
质量守恒定律
通过控制面S流出控制体的净质量流量 ＝控制体内质量减少的时间变化率
流体微团在流场中的 运动－物质导数的示 意图
物质导数（运动流体微团的时间变化率）
物质导数D/Dt与偏导数/t不同 ，/t是在固定点1时观 察密度变化的时间变化率，该变化由流场瞬间的起伏所引起。
流体微团在流场中的 运动－物质导数的示 意图
物质导数（运动流体微团的时间变化率）
计算流体力学常用的五大类数值方法简介
计算流体力学常用的五大类数值方法简介流体力学数值方法有很多种，其数学原理各不相同，但有二点是所有方法都具备的，即离散化和代数化。
总的来说其基本思想是：将原来连续的求解区域划分成网格或单元子区域，在其中设置有限个离散点(称为节点)，将求解区域中的连续函数离散为这些节点上的函数值；通过某种数学原理，将作为控制方程的偏微分方程转化为联系节点上待求函数值之间关系的代数方程(离散方程)，求解所建立起来的代表方程以获得求解函数的节点值。
不同的数值方法，其主要区别在于求解区域的离散方式和控制方程的离散方式上。
在流体力学数值方法中，应用比较广泛的是有限差分法、有限元法、边界元法、有限体积法和有限分析法，现简述如下。
一、有限差分法这是最早采用的数值方法，它是将求解区域划分为矩形或正交曲线网格，在网格线交点(即节点)上，将控制方程中的每一个微商用差商来代替，从而将连续函数的微分方程离散为网格节点上定义的差分方程，每个方程中包含了本节点及其附近一些节点上的待求函数值，通过求解这些代数方程就可获得所需的数值解。
有限差分法的优点是它建立在经典的数学逼近理论的基础上，容易为人们理解和接受；有限差分法的主要缺点是对于复杂流体区域的边界形状处理不方便，处理得不好将影响计算精度。
二、有限元法有限元法的基本原理是把适定的微分问题的解域进行离散化，将其剖分成相连结又互不重叠的具有一定规则几何形状的有限个子区域(如：在二维问题中可以划分为三角形或四边形；在三维问题中可以划分为四面体或六面体等)，这些子区域称之为单元，单元之间以节点相联结。
函数值被定义在节点上，在单元中选择基函数(又称插值函数)，以节点函数值与基函数的乘积的线性组合成单元的近似解来逼近单元中的真解。
利用古典变分方法(里兹法或伽辽金法)由单元分析建立单元的有限元方程，然后组合成总体有限元方程，考虑边界条件后进而求解。
由于单元的几何形状是规则的，因此在单元上构造基函数可以遵循相同的法则，每个单元的有限元方程都具有相同的形式，可以用标准化的格式表示，其求解步骤也就变得很规范，即使是求解域剖分各单元的尺寸大小不一样，其求解步骤也不用改变，这就为利用计算机编制通用程序进行求解带来了方便。
飞行器流动仿真讲稿第6章-计算流体力学的基本方法资料
x
uin, j
n
n
i1, j
i, j
x
n i, j
vin, j1 vin, j y
vin, j
n
n
i, j 1
i, j
y
用同样的方法可以计算u、v和e的一阶时间偏导数。
第6.2节 MacCormack显式两步格式
取泰勒级数的前两项计算密度的估计值（即时间层的一阶向 前差分）
综上所述，（6-10）式给出的密度ρ的二阶时间偏导 数可以计算出来，然后代入泰勒级数展开（6-5）式 就可显式地计算出t+Δt时刻密度ρ的值。
按照同样的方法，可以显式地计算出t+Δt时刻速度u、 v和e的值，然后用状态方程计算出压强p的值；
Lax-Wendroff格式思路清晰直观，但由于二阶时间 偏导数的推导和计算，需要繁琐的代数运算。
n i, j
t
t av
n i, j
t 2
t
n
i, j
t
n1
i, j
其它流动参数计算方法完全相同。
第6.2节 MacCormack显式两步格式
四、关于MacCormack格式的几点说明
在预估步和校正步中，空间导数的向前差分和向后差分可以 交换，即预估步使用向后差分，而校正步使用向前差分；
1
n i, j
pn i, j 1
pn i, j 1
2y
e n t i, j
uin, j
e e n i1, j
n i1, j
2x
vin, j
e e n i, j1
n i, j 1
2y
pn i, j
n i, j
un i 1,
北航计算流体力学第9课
对于迁移方程: 0t x u au += ()0a > (5-1) 1) Euler 显式格式110n n n n i i i i u u u u a t x+---+=∆∆ (5-2) 其等价微分方程为: (设tc ax∆=∆) ()()221126t x xx xxx a x a x u au c u c u ∆∆+=---+2) MacCormack 格式()()111112n n n n i i i i n ni i i i iu u c u u u u u c u u +-++⎧=--⎪⎨⎡⎤=+--⎪⎣⎦⎩ (5-6) 其等价微分方程为:()()232211624t x xxx xxxx a x a xu au c u c c u ∆∆+=--+-+3) Euler 隐式格式1111102n n n n i i i i u u u u a t x++++---+=∆∆ (5-7) 其等价微分方程为:2323236t x xx xxx a ta t a x u au u u ⎛⎫∆∆∆+=-++ ⎪⎝⎭4) Crank-Nicolson 格式11111110222n n n n n n i i i i i i u u a u u u u t x x ++++-+-⎛⎫---++= ⎪∆∆∆⎝⎭(5-8) 其等价微分方程为:32243224451261202480t x xxx x a t a x a x a t x a t u au u u ⎛⎫⎛⎫∆∆∆∆∆∆+=-+-+++⎪ ⎪⎝⎭⎝⎭一. 对于迁移方程0t x u au += (5-1)设通解为:n In t Imkxm u C e e ω= (5-10)式中I =k ——波数， ω——频率 将通解(5-10)代入方程(5-1),即求得t u In u ω=, x u Imku =代入(5-1), 整理得:0In u aImku ω+=即: n amk ω=- 代入通解式(5-10),得n aImkt Imkx m u C e e -=整理得 ()Imk x at nmu C e -= (5-11)式(5-11)表示各种不同波数(k 不同)的波都以同样的波速a 向下游迁移,如图1所示.所以把方程(5-1)称为迁移方程.tut∆()F x ()F x a t -∆xx at C-=图1. 0t x u au +=的解t x xxx u au u ε+=- ()0,0a ε>> (5-12)它的通解仍为n In t Imkxm u C e e ω= (5-10)求出()3, , t x xxx u In u u Imku u I mk u ω===-将它们代入方程(5-12), 得()3In u aImku I mk u ωε+=整理得 ()2n mk a mk ωε⎡⎤=--⎣⎦将上式代入通解(5-10), 得到(){}2Imk x a mk tn mu C eε⎡⎤--⎣⎦= (5-13)与迁移方程的通解()Imk x at nm u C e-=相比较,可见方程(5-12)的解表示不同波数的波以不同的波速()2a mk ε-向下游迁移, 波数越大, 迁移速度就越小, 如图2所示这种现象称为频散.utx321mk mk mk >>图2. t x xxx u au u ε+=-的解t x xx u au u α+= ()0,0a α>> (5-14)它的通解仍为n In t Imkxm u C e e ω= (5-10)求出 t u , x u 和()2xx u mk u =-代入方程(5-14), 整理得()2In aImk mk ωα=--代入通解(5-10), 整理得()()2mk t Imk x at n mu C ee α--= (5-15)上式表明, 不同波数波都以同一波速a 向下游迁移, 但在迁移的过程中, 由于()()20 mk tet α-→→∞, 初始扰动会越来越小. 这种现象称为耗散. 如图3所示.需要指出的是对于方程 t x xx u au u α+=- , 其通解为()()2mk t Imk x at n mu C ee α-=当t →∞时, ()2mk te α→∞, 即u →∞, 这种现象称为负耗散.u u xxttt∆t∆t∆t∆()()222mk teF x a t α-∆-∆()()2mk teF x a t α-∆-∆()F x ()F x ()()2mk teF x a t α∆-∆()()222mk teF x a t α∆-∆图3. t x xx u au u α+=的解图4. t x xx u au u α+=-的解由Taylor 公式234512624120i i x xx xxx xxxx xxxxx x x x x u u xu u u u u +∆∆∆∆=+∆+++++（1） 234512624120i i x xx xxx xxxx xxxxx x x x x u u xu u u u u -∆∆∆∆=-∆+-+-+（2）式（1）-（2）得中心差: 241126120i i x xxx xxxxx u u x x u u u x +--∆∆=---∆ (3)由式(1)得前差: 23412624120i i x xx xxx xxxx xxxxx u u x x x x u u u u u x +-∆∆∆∆=-----∆(4)由式(2)得:后差: 23412624120i i x xx xxx xxxx xxxxx u u x x x x u u u u u x --∆∆∆∆=+-+-+∆(5)由以上式子可见:1) 中心差分的截断误差项中只包含奇次导数项, 所以由中心差分构成的差分格式具有频散特性;2) 前差和后差的截断误差项中包含有所有导数项, 且第一项为xx u , 因此由前差和后差构成的格式具有耗散特性,且a) 当扰动从上游传向下游时, 即0a >, 后差格式为正耗散; b) 当扰动从下游传向上游时, 即0a <, 前差格式为正耗散.下面分析几个常见格式频散、耗散特性(以0u ua t x∂∂+=∂∂为例) 1. Euler 显式格式xu u at u u n i n i n i n i ∆--=∆--+11 等价微分方程为:()()221126t x xx xxx a x a x u au c u c u ∆∆+=---+显然, 它具有耗散特性(当1c ≤时) 2. MacCormack 格式()()()1111111112n n ni i i i n n n n i i i i n n n i i i u u c u u u u c u u u u u -+++++++⎧=--⎪⎪⎪=--⎨⎪⎪=+⎪⎩ 等价微分方程为:()()232211624t x xxx xxxx a x a xu au c u c c u ∆∆+=--+-+显然, 它是频散格式. 3. Crank-Nicolson 格式⎪⎪⎭⎫ ⎝⎛∆-+∆--=∆-+-++-++x u u x u u a t u u n i n i n i n i n i n i 2221111111 等价微分方程为:3224322441261202480t x xxx xxxxx a t a x a x a t x a t u au u u ⎛⎫⎛⎫∆∆∆∆∆∆+=-+-+++ ⎪ ⎪⎝⎭⎝⎭它是频散格式.对于Euler 中差格式1112n n n n i i i i u u u u at x++---=-∆∆ 其等价微分方程为:()222126t x xx xxx a t a x u au u c u ∆∆+=---+显然, 这是负耗散格式，解发散. 作业55.1 试写出下列微分方程的通解并分析其频散, 耗散特性1. 4 t x x u au u α+=± () 0,0a α>> 2. 6 t x x u au u α+=± () 0,0a α>> 3. 5 t x x u au u ε+=± () 0,0a ε>> 4. 7 t x x u au u ε+=± () 0,0a ε>>5.2 对于迁移方程(0 0t x u au a +=>) 试推导跳点格式的等价微分方程, 并指出它属于哪种格式.。
第3章流体力学连续性方程微分形式
第三节 流体动力学基本方程式
2.质量力 单位质量力在各坐标轴上分量为X,Y,Z，∴质量力为Xdxdydz x方向（牛顿第二运动定律
Fma ）：
d p p dx dx x ( p ) dydz ( p ) dydz X dxdydz dxd x 2 x 2 d
( u ) x x dx dydz
第三节 流体动力学基本方程式
X方向
( u x) dxdydz x
同理可得：
在dt时间内因密度变化而减少的 质量为：
y方向：
z方向：
( u y ) y dxdydz ( u z ) dxdydz z
dxdydz ( )dxdydz t dxdydz t
( u ) ( u ) ( u ) y x z 0 x y z
适用范围：理想、实际、可压缩、不可压缩的恒定流。
（2）不可压缩流体的连续性微分方程
当为不可压缩流时
0 t
Const
u u u y x z 0 x y z
三、粘性流体的运动微分方程
第四节 欧拉运动微分方程的积分
一、在势流条件下的积分
二、沿流线的积分
第三节 流体动力学基本方程式
一、连续性微分方程
在流场内取一微元六面体（如图），边长为dx,dy,dz，中心点O流速为 （ ux,uy,uz ） D' z C' ux dx ux dx A' dz u B' u ux x 2 z x x 2 o’ M uy ux N 以x轴方向为例： C D u 1 x dx dy u u dx 左表面流速 M A x x 2 B o u x x 1 右表面流速 u u dx N x 2 x y ∴ 单位时间内x方向流出流进的质量流量差： ( u ) ( u ) x x 1 1 M M [ u dx ] dydz [ u dx ] dydz x 右 左 x 2 x 2 x
第六章 计算流体力学的基本方法
Nanjing University of Technology
守恒形式
还用欧拉方程进行讨论。
二维流动的适用于CFD计算的守恒型方程：
U F G J t x y
（6-23）
显然，用麦考马克方法和拉克斯-温德罗夫方法，都
可以计算U的分量 、u 、v 、(e V 2 / 2) 在各时间步的
方程（6－24）中的向量F，它在网格点（i＋1，j） 处的值可从下式求出：
F i1 j
Fji
F x
av
x
（6-24）
19
Nanjing University of Technology
空间推进
通过预估－校正法得到这个平均值
预估步
➢用向前差分替代对y的导数:
F x
i
j
J
i j
Gi j 1
u x
p
v y
（6-4）
4
Nanjing University of Technology
拉克斯-温德罗夫方法
➢拉克斯-温德罗夫方法的基础是时间导数的 泰勒展开式。
➢任意选择一个流动参量，为明确起见，选
择密度 。
➢ t 时t刻，同一网格点(i,j)处的密度
可由tt i, j
泰勒级数给出：
tt i, j
(x)2 (y)2 2(y)2 2(x)2
麦考马克方法与拉克斯-温德罗夫方法对比：
➢麦考马克方法在预估步中用向前差分在校正步中用向后差 分，具有二阶精度，与拉克斯-温德罗夫方法具有同样的精 度。 ➢但是麦考马克方法不像拉克斯-温德罗夫方法那样需要计 算二阶时间导数，所以麦考马克方法更容易应用。
15
Nanjing University of Technology
解决工程流体力学问题的三种研究方法
解决工程流体力学问题的三种研究方法下载提示：该文档是本店铺精心编制而成的，希望大家下载后，能够帮助大家解决实际问题。
文档下载后可定制修改，请根据实际需要进行调整和使用，谢谢！本店铺为大家提供各种类型的实用资料，如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等，想了解不同资料格式和写法，敬请关注！Download tips: This document is carefully compiled by this editor. I hope that after you download it, it can help you solve practical problems. The document can be customized and modified after downloading, please adjust and use it according to actual needs, thank you! In addition, this shop provides you with various types of practical materials, such as educational essays, diary appreciation, sentence excerpts, ancient poems, classic articles, topic composition, work summary, word parsing, copy excerpts, other materials and so on, want to know different data formats and writing methods, please pay attention!解决工程流体力学问题的三种研究方法工程流体力学作为研究流体在工程领域中运动和作用的学科，涵盖了广泛的应用和复杂的理论。
第13章_计算流体力学CFD(5)
第二步：
时间步长为 ，空间导数采用中心差分，只对y的 导数采用隐式处理。
交替方向隐式(ADI)方法
第二步：
简化为三对角形式
交替方向隐式(ADI)方法
第二步：
对每一个固定的i，对所有 的j联立形成方程组。 对不同的i，重复上述过程。
交替方向隐式(ADI)方法
考虑二维热传导方程：
V 0
t 在可压缩流动中，右图 速度的棋盘分布经过一 个时间步就会被抹平。
右上角是u的值， 左下角是v的值
交错网格的应用
二维不可压流体压力梯 度采用中心差分：
压力会出现右图的 棋盘式分布
棋盘式的离散压力分 布
交错网格的应用
在交错网格上使用中心 差分就不会出现速度和 压力的棋盘式分布问题。
数值耗散、色散及人工粘性
偏微分方程（修正方程）：
尽管人工粘性降低了解的精度，但通常有助于提高解的稳 定性。
6.6 交替方向隐式(ADI)方法
交替方向隐式(ADI)方法
考虑二维热传导方程： 采用Crank-Nicolson方法（隐式）：
等号右端有五个未知量，不能得到三对角方程组，不能采 用托马斯算法（追赶法）求解。
MacCormack方法
在MacCormack方法中，预估步用向前差分， 校正步用向后差分；也可以预估步用向后差分， 校正步用向前差分。或者在时间推进解法的相 继两个时间步中轮流使用这两种办法。
6.3 粘性流动、守恒形式和空间推进
6.3.1 粘性流动
粘性流动
粘性流动的控制方程是N-S方程。
对定常流动，N-S方程的数学性质更多地表现为 椭圆型的，不能采用Lax-Wendroff方法和 MacCormack方法求解。
计算流体力学—周正贵[(P111-P125) P176]概要
![计算流体力学—周正贵[(P111-P125) P176]概要](https://img.taocdn.com/s3/m/ddce135e5acfa1c7aa00cca7.png)
7.3 网格生成图7.3（a ）为直角坐标系下求解域局部网格，按照节点序号给出网格线对应的贴体坐标数值1i -ξ，i ξ，1i +ξ和1j -η，j η，1j +η。
通常可简单地将节点序号作为在曲线坐标中的坐标值，比如j i j i ==ηξ，。
这样在计算域内即得到正方形网格，1=∆=∆ηξ，如图7.3（b ）所示。
如果在求解域内网格已分布完成，对应每一节点的x y 坐标值 21,,,,,，）（）（=j i y x j i j i 已知，则坐标变换关系可用节点坐标值的差分表达式表示，比如）（j i ,点雅可比可表示为ηξξηηξ∆-⋅∆-=-=-+-+-22)(1,1,,1,1,1,j i j i j i j i j i j i y y x x y x y x J ξη∆-⋅∆---+-+22,1,11,1,j i j i j i j i y y x x图 7.3（a ）求解域内局部网格 图7.3（b ）计算域内局部网格 因此要得到坐标变换关系必须在求解域内分布网格线，即进行网格生成，而后由这些网格线构成曲线贴体坐标。
合理的网格线分布不仅对计算精度有直接影响，甚至会影响计算过程的收敛性。
因此，网格生成是流场数值计算中一个比较关键的问题。
通常要求生成的网格线要平滑，避免有局部太大的扭曲，不同簇网格线要尽可能正交。
网格分布要与物理问题本身相匹配，也就是说疏密分布应与物理量变化率相适应。
比如在求解黏性绕流时，在壁面附面层区流动参数沿壁面法向变化很剧烈，因此沿此方向网格线要密集分布。
对于超音速流动问题在激波附近网格线要加密。
网格生成主要有三种方法：代数方法、微分方程方法和保角变换方法。
微分方程方法可以处理各种类型的不规则边界，具有较强的通用性，因而应用最为广泛；代数方法在工程实际中也有一定的应用；保角变换方法由于适用范围的局限性很大，目前实际应用较少。
在此仅介绍代数方法和微分方程方法。
7.3.1 代数生成方法代数方法实际上是一种插值方法，下面以叶栅通道内H型网格生成来说明这一方法的基本思想。
第13章_计算流体力学CFD(5)总结
空间推进
定常守恒型二维欧拉方程：
对于亚声速流动，上述 方程是椭圆型的，所有 空间推进方法都不适用， MacCormack方法也不 适用。
空间推进
定常守恒型二维欧拉方程：
对于超声速流动，上述方 程是双曲型的，空间推进 方法适用，MacCormack 方法也适用。
空间推进
定常守恒型二维欧拉方程：
MacCormack方法：
偏微分方程（修正方程）：
修正方程等号右端的项是截断误差，如果截断误差的主项 是偶数阶导数，数值解将主要表现出耗散行为；如果主项 是奇数阶导数，数值解将主要表现出色散行为。
数值耗散、色散及人工粘性
偏微分方程（修正方程）：
等号右端的偶数阶导数项起数值耗散的作用，奇数阶导数 项起数值色散的作用。
数值耗散、色散及人工粘性
速度修正量
可以从
得到。
压力修正法的基本原理
压力修正法本质上是一种迭代法，思路如下：
4） 用步骤3)中修正后的压力做为新的p*，回到步 骤2)。重复这个过程，直到速度场满足连续性方程 为止。
这样就得到修正好了的流场。
6.7.4 压力修正公式
压力修正公式
压力修正公式为：
压力修正公式
压力修正公式为：
SIMPLE算法的步 骤如下： 1）在右图所示的交 错网格上分别给出
p
* n
,
u
* n
,
v
* n
数值方法：SIMPLE方法
SIMPLE算法的步 骤如下： 2）求出 u
* n 1
,
v
* n 1
采用动量方程求解。
数值方法：SIMPLE方法
2）
u
计算程序_计算流体力学_对流方程_有限差分法_Lax格式_迎风格式_FTCS格式
% 一维对流方程迎风格式、Lax格式、FTCS格式差分法计算% 潭花林清华大学航天航空学院% FTCS格式对于一维对流方程不稳定，最好不用clcclear all% 1.参数定义dx=1;x1=-18;x2=18;x=x1:dx:x2;L1=length(x);% dt=0.5*dx; % 收敛dt=2*dx; % 不收敛t1=0;t2=t1+80*dt;t=t1:dt:t2;L2=length(t);alpha=1;lambda=alpha*dt/dx;geshi=1; % 迎风格式% geshi=2; % Lax格式% geshi=3; % FTCS格式% 2.显式求解zeta=zeros(L1,L2);for kk=1:3geshi=kk;for ii=1:L1if x(ii)>0zeta(ii,1)=1;elseif x(ii)==0zeta(ii,1)=1/2;elseif x(ii)<0zeta(ii,1)=0;endendendendif geshi==1for ii=2:L1for jj=1:(L2-1)zeta(ii,jj+1)=zeta(ii,jj)-lambda*(zeta(ii,jj)-z eta(ii-1,jj));endzeta(1,jj+1)=zeta(2,jj+1);endzeta1=zeta;else if geshi==2for ii=2:(L1-1)for jj=1:(L2-1)zeta(ii,jj+1)=(zeta(ii+1,jj)+zeta(ii-1,jj))/2-. ..lambda/2*(zeta(ii+1,jj)-zeta(ii-1,jj));endzeta(1,jj+1)=zeta(2,jj+1);zeta(L1,jj+1)=zeta(L1,jj)-lambda*(zeta(L1,jj)-z eta(L1-1,jj));endzeta2=zeta;else if geshi==3for ii=2:(L1-1)for jj=1:(L2-1)zeta(ii,jj+1)=zeta(ii,jj)-lambda/2*(zeta(ii+1,j j)-zeta(ii-1,jj));endzeta(1,jj+1)=zeta(2,jj+1);zeta(L1,jj+1)=zeta(L1,jj)-lambda*(zeta(L1,jj)-z eta(L1-1,jj));endzeta3=zeta;endendendend% 3.绘图% 3.1 t=0figure(1)n=1;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=0时刻的计算结果') %%% 标题axis([-18,18,-0.2,1.2])% 3.2 t=10figure(2)n=(10-t(1))/dt;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=10s时刻的计算结果') %%% 标题% 3.3 t=20figure(3)n=(20-t(1))/dt;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') %%% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=20s时刻的计算结果') %%% 标题% 3.4 t=40figure(4)n=(40-t(1))/dt;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=40s时刻的计算结果') %%% 标题聚乙烯（PE）简介1.1聚乙烯化学名称：聚乙烯英文名称：polyethylene，简称PE结构式：聚乙烯是乙烯经聚合制得的一种热塑性树脂，也包括乙烯与少量α-烯烃的共聚物。
流体力学公式总结
工程流体力学公式总结第二章 流体的主要物理性质流体的可压缩性计算、牛顿内摩擦定律的计算、粘度的三种表示方法。
1．密度 ρ = m /V2．重度 γ = G /V3．流体的密度和重度有以下的关系：γ = ρ g 或 ρ = γ/ g4．密度的倒数称为比体积，以υ表示υ = 1/ ρ = V/m5．流体的相对密度：d = γ流 /γ水 = ρ流 /ρ水6．热膨胀性7．压缩性. 体积压缩率κ8．体积模量9．流体层接触面上的内摩擦力10．单位面积上的内摩擦力（切应力）（牛顿内摩擦定律）11．.动力粘度μ：12．运动粘度ν ：ν = μ/ρ13．恩氏粘度°E ：°E = t 1 / t 2第三章 流体静力学重点：流体静压强特性、欧拉平衡微分方程式、等压面方程及其、流体静力学基本方程意义及其计算、压强关系换算、相对静止状态流体的压强计算、流体静压力的计算（压力体）。
1．常见的质量力：重力ΔW = Δmg 、直线运动惯性力ΔFI = Δm·a离心惯性力ΔFR = Δm·r ω2 .T VV ∆∆=1αp VV ∆∆-=1κV P V K ∆∆-=κ1n A F d d υμ=dnd v μτ±=n v d /d τμ=2．质量力为F 。
：F = m ·am = m (f xi+f yj+f zk)am = F /m = f xi+f yj+f zk 为单位质量力，在数值上就等于加速度实例：重力场中的流体只受到地球引力的作用，取z 轴铅垂向上，xoy 为水平面，则单位质量力在x 、y 、 z 轴上的分量为fx = 0 , fy = 0 , fz = -mg /m = -g式中负号表示重力加速度g 与坐标轴z 方向相反3流体静压强不是矢量，而是标量，仅是坐标的连续函数。
即：p = p (x ,y ,z )，由此得静压强的全微分为:4．欧拉平衡微分方程式单位质量流体的力平衡方程为：5．压强差公式（欧拉平衡微分方程式综合形式）6．质量力的势函数7．重力场中平衡流体的质量力势函数z z p y y p x x p p d d d d ∂∂∂∂∂∂++=d d d d d d 0x p f x y z x y z x∂∂-=ρd d d d d d 0y p f x y z x y z y ∂∂-=ρd d d d d d 0z p f x y z x y z z∂∂-=ρ01=∂∂-x p f x ρ10y p f y ∂∂-=ρ01=∂∂-z p f z ρz z p y y p x x p z f y f x f z y x d d d )d d d (∂∂+∂∂+∂∂=++ρ)d d d (d z f y f x f p z y x ++=ρd (d d d )x y z p f x f y f z dU ρ=++=ρd d d d x y z U U U U x y z =f dx f dy f dz x y z gdz ∂∂∂∂∂∂=++++=-积分得：U = -gz + c*注：旋势判断：有旋无势流函数是否满足拉普拉斯方程：22220x y ψψ∂∂+=∂∂8．等压面微分方程式 .fx d x + fy d y + fz d z = 09．流体静力学基本方程对于不可压缩流体，ρ = 常数。
TJU计算流体力学 计算流体力学
50.00
55.00
7.50 7.00 6.50 6.00 5.50 5.00 4.50 4.00 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 -0.50 -1.00
• 五洼地势
130
120
110
100
90
80
70
60 10 20 30 40 50 60 70 80 90 100 110 120 130 140
6.00 5.50 5.00 4.50 4.00 3.50 3.00 2.50 2.00 1.50 1.00 0.50 0.00 -0.50 -1.00 -1.50 -2.00
25
20
15
Distance /km
10
5
0 0 2 4 6 8 10 12 Distance /km 14 16 18 20
2）静止时不能承受切应力，在微小切应力 作用下，流体会发生不断变形。 3）静止流体内不能抵抗拉力，只能抵抗压 力。
计算流体力学中的基本问题
CFD，Computational Fluid Dynamics
• • • •
基本方程 计算方法 程序设计 前后处理
目的：解释物理现象，研究变化趋势
物理现象和研究重点 基本假定 微积分原理 微体分析 有限体分析 质量守恒、能量守恒和受力平衡 微分方程 数学方法 理论公式
第一章
• • • • • • • • • • • • • 1.2流体力学方程的表达方式 1.2.1坐标系 1.2.2向量表达方式 1.2.2.1连续方程的简化 1.2.2.2运动方程的简化（沿曲线） 1.2.2.3涡流方程 1.2.2.4理想流体方程的积分 1.2.2.5三维流体方程的垂向积分 1.3张量表达方式 1.3.1字母、足标及求和约定 1.3.2张量流体力学方程 1.3.3用多重尺度推倒水波方程 1.4湍流方程
计算流体力学课件完整版
●实验要受测量技术限制，实验周期长、费用高。
☆ 理论研究 ●在研究流体流动规律的基础上，建立了流体流动基 本方程。 ●对于一些简单流动，通过简化求出研究问题的解析 解。
计算流体力学
●对于实际流动问题，通常需运用流体力学基本方程， 借助于计算机求数值解（计算机数值模拟）— 计算流体力学CFD。
Z
skirt.plt X Y
75 50 25
0 -25 -50 -75
-2
Y(M) 0
2
0 2 4 6 10 8 X(M) 12 14
D) 16 Feb 2003 Velocity Vectors
4.5
4 velocity.plt
3.5
3
2.5
2
1.5
Z
Z
(3D) 16 Feb 2003 IJK-Ordered DZ ata
ijkcyl.plt X Y
Z
-0.4 -0.2 Y0 0.2 0.4
1
0.8
0.6
0.4
0.2
0 -0.4 -0.2 0 X 0.2 0.4
Z
jetflow.plXt Y
0.6 0.5 0.4 0.3 0.2 0.1
0 0 Y0.1 0.2
-0.6 -0.4 -0.2 0 X 0.2 0.4 0.6
轴流叶轮计算与实验叶片表面极限流线
计算流体力学
轴流叶轮计算与实验性能比较
计算流体力学
轴流叶轮计算与实验流场结构比较
计算流体力学
第二章 流体力学数值计算数学模型及定解条件
☆本章所涉及的基本方程有两类： ●流体力学基本方程，基本出发点：质量守恒、动量守恒和能
计算流体力学液体表面追踪方法分类-力学论文-物理论文
计算流体力学液体表面追踪方法分类-力学论文-物理论文——文章均为WORD文档，下载后可直接编辑使用亦可打印——1、概述流体模拟一直是CFD领域的一个热点话题，在模拟过程中，流体现象可粗略分为两类：一类为低速流体，包括滴水涟漪、风吹水面、水波叠加等；这类流体最主要的特点是流速较缓，给人以足够的响应时间，流动过程中不会出现大的扭曲和变形，因此抽象出来的物理方程计算量较小，模拟过程更容易控制；另外一类是高速流体，包括液面撞击时产生的破碎、飞溅、翻滚等现象，这类流体的特点是流速较快，要求计算速度比精度更为重要，流动过程中液面产生的扭曲和大变形使液面追踪的物理描述更为复杂，且计算量较大。
2、液体表面追踪方法分类液体表面追踪方法分为两类：网格方法和无网格方法。
这两种方法按采用的坐标系不同又分为拉格朗日法、欧拉法、拉格朗日欧拉混合法；如下所述。
2.1 网格方法2.1.1拉格朗日法在拉格朗日法中，计算网格固定在液体表面，网格的节点即为液体的物质点，因此网格随液体一起运动和变形，液体不会在单元格之间流动，所以该方法能够自然的处理自由表面和物质分界面问题；但对于液体表面的大变形现象，网格会随之发生扭曲和畸变，导致计算结果不稳定。
2.1.2欧拉法在欧拉法中，计算网格以空间坐标为基础划分，网格与液体相互，计算过程中网格固定不变，不存在畸变和相交问题，液体流过单元格空间，空间内的流体状态发生变化。
所以在模拟过程中，液体界面的捕捉相对困难，因此欧拉法通过引入不同的液面追踪模型来离散化控制方程，这种模型的选取决定了数值计算的准确性。
欧拉法的典型代表有MAC（标志网格法）和VOF（流体体积法）。
在MAC法中，流体空间使用欧拉坐标系划分为矩形网格空间，初始时刻，在每个含有流体的网格空间内设置若干个无质量的标记点，标记点以其所在的流场的速度移动，它本身并不参与流场的计算，标记点的移动描绘了整个流场的流动。
模拟时认为含有标记点的网格即为含有流体的网格，因此含有标记点和不含有标记点相邻的网格即为液体的自由面网格，这些网格构成了液体自由面的形状。
计算流体力学4
中心节点P处的压力值没有出现在式(3.4)或(3.5)中。将图3.2 所示的压力场分布代入式(3.4)和(3.5)，离散后的压力梯度在任 何节点处均为0，尽管实际上在空间两个方向上均存在明显的压 力振荡。这样，该压力场将导致在离散后的动量方程中，由压 力产生的源项为0，与均匀压力场所产生的结果完全一样。这显 然是不符合实际的。 同样，若在流场迭代求解过程的某一层次上，在压力场的当 前值中加上一个锯齿状的压力波，则动量方程的离散形式无法 把这一不合理的分量检测出来。它会一直保留到造迭代过程收 敛而且被作为正确的压力场输出(图3.3中的虚线)。
当计算中流体的密度、能量、动量等参数存在相互依赖 关系时，采用耦合式解法具有很大优势，如在3.1.1节中提 到的求解关于u,v,w,p共4个变量的方程组。其主要应用包 括高速可压缩流动、有限速率反应模型等。 耦合式解法中，所有变量整场联立求解应用较普遍，求 解速度较快，而在局部对所有变量联立求解仅用于声变量 动态性极强的场合，如激波捕捉。 设计算区域内的节点数为N，则每一时间步内须求解4N 个方程构成的代数方程组（三个速度方程及一个压力或密 度方程)。总体而言，耦合式解法计算效率较低、内存消耗 大。
为了解决因压力所带来的流场求解难题，人们提出了若 干从控制方程中消去压力的方法。例如，在二维问题中，通 过交叉微分，从两个动量方程中可消去压力，然后可取涡量 和流函数作为变量来求解流场。 涡量—流函数方法成功地解决了直接求解压力所带来的 问题，且在某些边界上，可较容易地给定边界条件，但它也 存在一些明显的弱点，如壁面上的涡量值很难给定，计算量 及存储空间都很大，对于三维问题，自变量为6个，其复杂 性可能越过上述直接求解(u,v,p)的方程组。因此，这类方法 在目前工程中使用并不普遍，而使用最广泛的是求解原始变 量(u,v,p)的分离式解法。 基于原始变量的分离式(segregated)解法的主要思路是：顺 序地、逐个地求解各变量代数方程组，这是相对于联立求解 方程组的藕合式（coupled method）解法而言的。目前使用最 为广泛的是1972年由Patanker和Splding提出的SIMPLE算法。 这种方法将是本章重点介绍的方法。