计算水动力学.docx

《计算水动力学.docx》由会员分享，可在线阅读，更多相关《计算水动力学.docx（57页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 质量守恒方程(mass conservation equation)单位时间内微元体中流体质量的增加可同一时间间隔内流入该微元体的净质量bp工。（夕）工。（夕酎）工_n111 Vbp工。（夕）工。（夕酎）工_n111 Vdt dx引入矢量符号动量方程dy dzV -a=div(a) = L+2- dx dy1 + V.(pu) = 0dz峰后继可压流体的质量守恒方程，又称连续方程(P)/ % 6% 6 屯二+ 加 pun) = -+ + + + F dtdx dx dv &等+而,（m）一dv Ox 以， dz + 而，(PWn)= 一XZdp drx.Gt _+- + dz dx d dz

2、对于牛顿流体，粘性应力T与流体的变形率成比例%=29+2加(11)C Ar =2/ + z7r(u) 办t=2/=+ 2d4(u) dz动力粘度(dynamic viscosity)2 第二粘度(Second viscosity), = r = lixy 户卜(du c、十 l Sr；5u 加1股取2 = 2/3=r =XZ 2X/ 一十一(5z dx产 =y另一求取导数近似值的方法是将函数拟合成一插值函 数，然后对其求导。例如，假设用分段线性插值，我们 可以得到FDS或BDS近似(取决于笫2个点位于点片的 左边还是右边)例如，利用点X-、x/IIxn将函数拟合成抛物 线，那么由此插值函数可求

3、得看处的一阶导数%+14州(4占 + 人叫+1)%+14州(4占 + 人叫+1)&+1(乙片)2 。一 (A勺+1)2 + 创(A&+1)2 -(勺)2An = Xi - Xi-1(| =(|)1+a1(x-x1) + a2(x-xi)2儿+i =G + Axi+i ai + Axti a2(解出a和a2 % 二虫-Axt a】 +Ax： a7其它多项式，如样条函数(spline)等，亦可以用作 构造插值函数，然后求得导数近似。一般，一阶导数近 似值具有与所用的多项式插值函数的次数(degree)相同 的阶数(Oder)的截断误差。一次多项式相应的一阶导数 近似值具有一阶的截断误差。二阶导数近

4、似值BDS BDSBDS BDSFDS ZJ1) + 仇一 1(3+1 - Xi)仇(&+1 Xj-1)(g+i Xi)2xi -Xi-i)CDS CDSCDSd+i(3 - &-1) + 血 i(%+i - Xi)-仇(数+i - 3i)紧凑型格紧凑型格，(勺+1 -0-i)(&+i - 0)(勺-勺一1)式(compact schemes)基于多项式拟合，使用邻近节点的值和导数值来确定系数:一般地，紧凑型格式可以写成:()i+1=a1+2a2(x1+1-xi) + 3a3(xl+1-x1)2+.有限差分方程的建立基本思想,上调攵网格点上的变量的差分形式代替微分，结果得到一个包括所有节点上未

5、知量的联立代数方程组计算流体力学中应用最多的离散化数值方法作为计算技术它历史最悠久，理论上相对成熟1）直接用差商逼近代入微分方程中的所有各阶偏导数分别选择适当的差商逼近，并考虑逼近的截断误差精度，从而 将微分方程改写为代数的差分方程;同时得到整个差分方程对微分方程的逼近的精度。 基本步骤-借助Taylor级数展开给出各阶导数的差商表达式-将方程中的各阶导数用相应的差商表达式代替-整理化简d2n 例：dx 2 = 00x 0时,迎风格式。 ii）仍然采用线性拟合，但采用A、C两点进行插值llD = 4 + (AX - 4/V)11 N 2AxABC1+1 =u-_1+|(l-c)(u*1+1-u

6、_1)整理 =;(%+%或j 2“为1旦=0(Lax格式)N2 Axiii)线性拟合，采用B、C两点外插,iv)抛物线拟合，用A、B、C三点的值进行二次曲线拟合 一设：过A、B、C三点uA、uB、uC的值，为二次抛物线PG)=。+A、B、C三点的x坐标可以简单地给为：- Ax , 0, AxiiA 二% - Ax % + Ax = %uc = % + Ax% + Ax2%，+制(一)+为2AX(-r/A/)22AX-n+i n aAt z 1 aAtaAt1 Ax)(ua-2ub+uc)= u；_aAtAxu； -1嗜)-1 aAt2 Ax%2u；+%整理得= u；-u7。哈味-211；+喷乙

7、u； -u：-+aAtii01嗜)a2AtAx2 Ax2i嗜-2u：+%(一Lax-Wendroff格式) 其出发方程可认为：du du a1 Z S、/、 1- ci =. 5 dx 2 dx3)分裂差分算子的离散化方法假设微分方程中的微分算子L(u)可作“和”分裂，即Z(w) = Za(w)+Z2 (w).那么微分方程詈=z（，）的差分方程可由i = ； ： (% %) = 1 c(即Step 2: /i =甲+ $闷 -2下+ 标=(1 +阂卜下 JJJJ J 、J Jdu du d2n1- a = cr-dt dx dx“和”式分裂为wdudta-二0dxdu d2n-Crr d d2

8、L = -aher-dx dx-(2)可按卜列格式计算（例取FTCS格式）Step 1:Step 2:理论依据 由以上两步，可消去中间步结果.即：%+1=(1 +s 电1 4“、.)；=1-c(“y)+s5；csC;；最后一项容易验证cs6/ = F。鸟0(A.3 Ax AV二0(加 2)对于时间前差格式：=七=。(加)上述格式与源方程的误差项cs 5 与TE同阶,可并入巳考虑，而不 影响差分的最后精度。证毕。差分方程的相容性、收敛性及稳定性: 实例：对流现象的方程和边界条件0tT,0x111”。，J?一v正+乂芯）存1+ 小）=At/Ax趋于0才相容，即AtvvAx截断误差E：=vAx2（?

9、占）:12m-uh 7 - 2acos0 +Vl-4a2 sill2 9放大因子 G二l + 2c稳定条件无条件稳定要求： Z/Zf。2）、对流方程的差分格式：招+ 含二FTCS差分格式ii+1 -11? u? - U? .J J I a 川 口广1At2 Ax时间方向用前向差分, 空间方向用中心差分稳定性分析得放大因子G = 1 - i s sin 61稳定条件：无条件不稳定迎风差格式十a + | 瑶二呜 t2 ZkrCourant 数+0(显式格式)稳定性分析得放大因子G = (l- 2| sin2 -)-/5 sin 02时间方向用前向差分, 空间 方向 用 迎风差分稳定条件：s 1(C

10、FL条件)C ourant-Friediichs-Lewy蛙跳(leap-frog)格式喝+ 一 吗 7 ,:+1 - 1稳定性分析得放大因子G = (1-52 sin2 0)1 2 - z 5 sin 0 稳定条件：s4=偶数的网格节点(j,n)称为“偶点”，把n+j=奇数的 网格节点(pi)称为“奇点”。用蛙跳格式计算时，奇点和偶点将各 自形成一组解。其结果在邻点之间有显著的“波动”现象，这当然 不是所期望的。Roache(1972)提出，采用平均n层和n+1层上的解 来克服这个问题。 Lax-Wendroff(LW)格式蛙跳格式是三层格式，在计算1严1时要用到iP和1pL因而 计算存贮量

11、要比二层格式大。其次，在起算时必须给出二层初 始数据u。和u1通常第二层8=1 )的计算选用具有同样精度的二 层格式，如Lax-Wendioff(LW)格式。c$2w；+1 =国(泊一吟-】)+ 另(以加一2H +认为duXdtNt1 A d2U1 2 a f ll-Ar=1a、?2 dt1 t 2 ax2稳定性分析得放大因子G 二(1 -2s2 sin2 -)-issin2 0稳定条件：S0 7/x2/x稳定性分析得放大因子l-z（O.55sin0）、G + i（0.5ssme） 稳定条件：无条件稳定截断误考是。（汽ArD3）、对流扩散方程的差分格式： FTCS差分格式a | Mx；+1 整

12、理得产=-Ar2稳定条件:一 ）；十 1 + （1 2。）；十 + 5）_ 乙乙（显式格式）1. 2夫= V（格式Re数）Zkr2I。应xVon Neumann稳定性分析&产-4-(27 一 $)；十 1 + (1 2)； + s)吗一 1乙乙：二 EA* 7 1 二 YAeJkXj= 用/6-)JJJa | Ai4+1 代入整理并求G二才得:G = l + 2a(coskAx-l)-is sinkAx -1 + 2a(cos9-1)-/5 sin9g2如se-D蓝配篝黑胪G = 1 + 2ct（cos9-1）-/5 sin 9（此式在复平面上表示实半轴 为2o, 虚半轴为s的椭圆）G = 1

13、 + 2ct（cos9-1）-/5 sin 9（此式在复平面上表示实半轴 为2o, 虚半轴为s的椭圆）显然应:2al5 1(5 一2s 1b22?力对椭期二a7+4=1其顶点的曲率半经分别为: b-故还应有：故还应有：5 n+l .一% I9十皿% 一 为 士% (I -3)二x2J稳定条件：1a&（ - ）-1）一声明1 2%+吟7）_ 0e = sgn(a)-稳定条件：“京 2+，c 02 = 0a 0F _ kL仙+ a |2kr)= 037)”（粘性项增加但稳定性改善）蛙跳(leap, frog)格式%+一谤-1U v-j2稳定条件：长=华 5 =譬=想2、一维热传导方程相容性、稳定性

14、的推导Lax等价定理在采用一种差分格式写成差分方程并讨论其相容、收敛与稳定的性质时可以看到， 相容性是比拟容易证明的，只须用Taylor级数展开就行；稳定性的证明以后可以介 绍诸如Von Neumann法等来解决，而收敛性证明一般比拟麻烦，但Lax等价定理很好 地解决了这个难题。Lax定理：对给定的适定的线性初值问题，用有限差分方程去近似时，假设差分方程 满足相容性条件，那么其稳定性是其收敛性的充分必要条件。简单地表达为： 稳定性器热攵敛件初竹，相谷这个定律的作用是明显的。对适定的线性初值问题，只须证明其相容性与稳定性， 收敛性自然满足，从而省却证明收敛性的麻烦。一维热传导方程：比0,0） =

15、/x）=6。） u（Lri）=倏（力 采用FTCS差分格式时间方向用前向差分, 空间方向用中心差分十一；泊 一2+_】/ 、八 1 c r J区 # S20;尸 12，- 1) 得=与tOJ = f(与)G =,J)% 配0/) =, S = 1,2,.N) 吗=我(GJ吗=+ （1 2。） + 噌一（显式格式）=声由时间n层通过简单运算就可以得到时间n+1层上的u。（1）、相容性及精度推导Taylor公式瓦;+i = /Taylor公式瓦;+i = /代入可”=办；+1 + （1 2。） +。吟一【储）:+。3）=（养+。叱）上式说明FTCS格式在网格加密直至间距趋于零时就趋向于原微分方程差

16、分格式与 原微分方程相容的。将上式改写为：+039它对于时间是一阶的，对于空间是二阶的。在时间方向上具有一阶精度，在空 间方向上具有二阶精度。（2）、稳定性推导Von Neumann稳定性分析（Fourier分析）：Fourier级数展开：；=收回）， ;口无X=1An是波数为kx（波长入二2兀依）的分量在时间n层的振幅，与x无关。定义相位角。二kx x,代入差分方程，对于每个Fourier分量有：=/普 +o(/(De +4%心TN -2/6)消去评:An+l =Anl + o(ee + e-70 - 2)放大因子:1cos(x) = etx+e.-lx)JW+1G = = 1 + o(/

17、+ 1汨-2)G = 1 + 2 er (cos 0 - 1) = 1 - 4cr sin 2 因|G|1时，差分格式稳定；所以，上式中OK 0.5满足该要求, 即选取的At、Ax应该满足o = vAt/Ax2 Vo.5判稳小结：1、将任一Fouiiei分量加代入差分方程2、求放大因子G=An+i/A 并令|G|13、求出应该满足|G| 1的条件x第四章小结1、FVM的原理将控制方程在控制容积上积分从而得到离散化方程的离散化方法。具体步骤：将控制方程在控制容积上积分；假定适当的分布函数(distribution function)-阶梯分布-线性分布将分布函数代入并完成积分，整理化简得离散化方

18、程。第六章小结1、有限分析法FAM基本思想：把求解区域划分成许多矩形网格，四个网格组成一个单元；在局部单 元内将微分方程线性化，在单元边界上为一近似函数；局部单元内求解微分方程的 解析解；建立单元中心的和其周围八个结点之间的迭代关系式。格式推导：对象：对流扩散方程标准的椭圆型方程：2A(|)x +2B(|)y Rxx +(|)yy +G步骤：1)、线性化 2)、齐次化2 A耙 +2B% =幄 +（|）” +G2/0 +2B*=a + %G（44 + /）2（4 +炉）瑜（X）一A Jd7W Nf NEWC。（。,0）加。叫一1一EC lxSCSE么（x）3）、边界近似：假设变量近似满足指数多项

19、式分布。w = aw（e2HfJ - i） + w7! + cw ， = -h （j）E = aE （e2BtJ -1） + bErj + cE ,占=h .s =%（/.-1） + 44 + % , r/ = -k（f）N = aN （e2A -1） + bN + cN , r/ = k4）、有限分析解80（O,O） = ZCM+CpG /1=1=Cne 如e+ Cnw如 w+ Cnc 如 c + Cse巾se + Csc,sc + Cec4ec +CpG式中：Cn为有限分析系数A B h k （dx dy）C最=一2 cos1l4/7PBk e0 2 cosh 魔、Ah ,网 e +e4c

20、cb1l4/?cos1lBI（1）Ce-Ah . -Bke +e4coslL4/?cosliBt（14月）=/出 2 跳。C SE-Ah . Bk e +e4cos1l4/7cosIiBA:(1-i)二产%Ah -Bk-Bk e,NC 2ccdiBk2cosIlMR =e-2”翳爵（ooP2 = 4A/zcosh A/zcosh Bk cothm=(1 尸 4,(A2Tz2 + / )2 cosh(t Jr + 炉 + %h = (m - 0.5)CpCp2(A2 +B2)aa(c/ve + CEC + CSE Gwv一 wc + sw )+ Bk(Cne + Cnc + Caw se sc

21、Csw)1hktanh Ah(l 一 己)=-tanh Bk(l 一召) 2A- 2B（3）优、缺点优点：1）、稳定性好。有限分析系数的大小取决于单元雷诺数Ah、Bk及网格长宽比k/h=0； 有限分析系数 CSC、CNC、CWC、CEC、CSW、CSE、CNW、CNE 及 Pl、 P2、E2的值都在。和1之间；八个系数之和等于1这一特性能保证迭代格式 具有较好的稳定性。2）、对称性好。由于椭圆型方程及矩形网格单元都具有对称性，因而有限分析系数 也具有许多对称性。NEECNW、 NCX !、 I、 iWC 7SE当A2=-A1时，两种情况的有限分析系数关于NS轴对称；当B2=B1时，两种 情况的

22、有限分析系数关于WE轴对称；当A2=A1B2=B1时，两种情况的有限分析 系数关于NW-SE轴对称；当A2=B1、B2=A1、。=1/01时，两种情况的有限分析 系数关于SW-NE轴对称；根据这一特性，在分析计算时，我们只需考虑A、B0 的情况，其它情况可以由对称性得到。3）、具有自动迎风效应。网格雷诺数在运动方程中，相当于水流的流速分量，因此, 当水流从西南方向流来，流向（风向）为东北方向，西南方向点对中心点的影响最 大，相应有限分析系数CWS （权系数）最大，其次为CSC、CWC两点，并且随着 Ah、Bk的增大，对流作用逐渐增大，西南方向点对中心点的影响逐渐最大，有限 分析系数CWS逐渐增

23、大。缺点：指数、级数、交错级数，每次计算都要耗费较多的计算时间，在单元雷诺数较大 时（如Ah、Bk35）,直接采用公式计算，会出现错误的负系数或大于1的正系数, 这与上述有限分析系数的特性相矛盾。造成这一错误，并不是有限分析法的理论问题，而是有限分析系数的计算问题。 随着单元雷诺数增大，交错级数收敛所需的项数和所需的有效位数也逐渐增多，一 般计算机的八位有效为数不能满足计算要求，即使双精度变量，也不能解决根本问 题。另外，一般计算机的数值范围为（10-381038）,在指数计算时也应注意防止计 算机溢出。2、混和有限分析法HFAM（1）基本思想：在局部单元线性化；维对流扩散方程的有限分析解；线

24、形算子迭加原理；导出二维对流扩散方程的有限分析解；建立单元中心点和其周围四个结点之 间的迭代关系式。格式推导：线性化一维对流扩散方程：2A6-限二GxXi-i Xi汨+/I11 X图5-4 单元网格图一维分析解：0 = C +会 工/_,%,.+式中系数Cl C2由单元两边界结点的函数值储/M+I确定利用线形算子迭加原理，在单元内将二维对流扩散方程变为xeXi.bx i+iyeyj-i，yj+ixexi.bx i+i,yeyj-i,yj+iLxO=2A(j)x -限=GxLy(|)=2B(|)y -(j)yy =GyL(|)=Ly(|)+Ly(|)=G% =，C力 n+CpGH=1-CncMc

25、 + Cscsc+ CwcOwc + CecOec +CpG式中为混合有限分析系数，反映了周围四点对中心的影响的全系数。C 一妁从 p. _ Bh SINH (Ah) 2 SINH(Bk)P = 2Ak coth Ah + 2B/icoth Bk（3）优、缺点1）、防止了交错级数带来的问题，2）、指数溢出问题依然存在，3）、不能反映单元四 个角点对中心点的影响。4）、对称性、迎风效应：网格雷诺数Ah、Bk在运动方程 中，相当于水流的流速分量。当水流从西南方向流来，流向（风向）为东北方向。随着Ah、Bk的增大，对流作用逐渐增大，西方和南方两点对中心点的影响逐渐增 大，相应有限分析系数CSC、CW

26、C也逐渐最大。自动迎风效应只是在Ah、Bk较 小时有效，当Ah、Bk5时，CSCeCWCM）.5,混合有限分析法又蜕化为一维简单迎 风格式。3、有限近似法FPM基本思想：1）、在网格单元内构造方程的基本解;常数4、一次线性函数为x+4V 双曲型函数。4, 一产）三次函数纬（V 一392 +/-3x）2）、由这些基本解构造出一个函数；在六边形网格区域内满足拉普拉斯方程，式中有六个系数区（i=l,6）。3）、由六边形6个节点值/ = L6）来确定6个系数叼4）、确定中心点的值。二弓格式推导：v - v V 0对于粘性为常数的不可压流体x y z u Reynolds平均方程基本思想湍流山两种流动形

27、式叠加而成，一是时间平均流动，二是瞬时脉 动流动。p = Xp Reynolds平均控制方程组连续方程普+ V.(pu) = 0动量方程l + div(puu) dt=grad u)-+8(/w)dx dydz+ dfv( pvu)= %v(/z gtad v)一年 +dx dy dzIfdzv(xwu)5(/?vwd+ S非守恒型控制方程的推导：守恒与非守恒形式的控制方程可P0dt守恒与非守恒控制方程的形式+ divpu(p) - divV grad e)+ 5守恒型方程(conservative form)对流项写作散度形式物理量都在微分符号内正六边形有限近似解7点格式(系数计算):cp

28、= a、+ a2x + a3y + a4(x2 -y2) + a5xy + a6(x3 -3xy2)23 9|=。+2r+4厂 +6 厂(P) % + 4 生 +a3r 42 +5/20 3 = % 4 2 +a3r 4 厂223(p4 = ax -a2r + aAr -a6r一 a2r 2 a3r - y 为/ +a5r2 +、。6 =% + 2a2ra3r2a4r2 -? a5r2 a(6己解方程组得4 = % C(p，G =卜 =1,2,6z=l=1特点6 在系数矩阵均中，第一列6个元素%均等于1,第一列6个元素对应的代数余子式之和541 正好等于系数矩阵行列式可，由此可以得到 去&=1

29、,说明该有限近似六点格式具有系77 %数之和为1的特性。有限体积法：FVM的原理将控制方程在控制容积上积分从而得到离散化方程的离散化方法 具体步骤将控制方程在控制容积上积分；假定适当的分布函数(distribution function)-阶梯分布_线性分布将分布函数代入并完成积分，整理化简得离散化方程阶梯型分布函数控制容积上均匀分布(为一常数)控制容积代表点(节点)处的值为分布值: 节点间线性分布：小尸+(x _ X/ )(X-Xp)xexw,xpxexp,xE分布函数说明：梯形分布主要用于- 计算控制容积上待求变量的值- 源项，非导数项- 非稳定项线性分布主要用于- 待求变量的梯度值- 控

30、制界面处待求变量值分布函数说明：梯形分布主要用于- 计算控制容积上待求变量的值- 源项，非导数项- 非稳定项线性分布主要用于- 待求变量的梯度值- 控制界面处待求变量值网格参数：名称与定义网格参数：名称与定义 (8x)w = (8x) + (M 节点W - P之间的距离 (&)e = (&)+e + (&)节点P - E之间的距离 (8)+“，控制界面W-节点P之间的距离 (收)1节点P-控制界面e之间的距离 Ax = (3x)=,+(取)7 控制容积 桃，e 左、右控制面 (8x)w = (8x) + (M 节点W - P之间的距离 (&)e = (&)+e + (&)节点P - E之间的距离 (8)+“，控制界面W-节点P之间的距离 (收)1节点P-控制界面e之间的距离 Ax = (3x)=,+(取)7 控制容积 桃，e 左、右控制面将方程在控制容积PXw , Xj上积分dx+ PStZv = OJw1(一)J.(&r)1 f.吗公 w dx