计算流体力学

1、计算流体力学作业第一部分用微元法推导出连续性方程，N-S方程，以及牛顿流体的化简形式解：由质量守恒定律，单位时间内流体微元体中质量的增加，等于同一时间间隔内流入微元体的净质量，得质量守恒方程：t+(u)zu)x+(u)y=0即为连续性方程：t+u=0动量守恒定律：微元体中流体的动量对时间的变化率等于外界作用在该微元体上各种力之和。Fx=max以Fx为例，Fx分为体积力（Body forces），与表面力（Surface forces），表面力又分为压力（Pressure）以及粘性应力（Viscous），其中粘性应力分为普通的（如xx）与剪切力（如xy，xy表示在y方向上作用于与x所垂直的平面的

2、力）x方向的体积力为：fxdzdydz在abcd平面上，x方向的表面力为yxdxdz与abcd相平行，在efgh平面上，x方向的表面力为yx+yxydydxdz注意二者方向相反，综合为yx+yxydy-yxdxdz，为x轴正方向同理可得另外两个平行平面的的合力。另：x方向的压力为：pdydzx负方向的压力为p+pxdxdydz压力合力为x方向：p-p+pxdxdydz综上，Fx=-px+xxx+yxy+zxzdxdydz+fxdxdydz微元体的质量：m=dxdydz微元体加速度：ax=DuDt则DuDt=-px+xxx+yxy+zxz+fx将式子左侧变形为全微分形式DuDt=ut+Vu其中，

3、ut=(u)t-ut且uV=uV+Vu则DuDt=(u)t+uV如此可得x方向的N-S方程：(u)t+uV=-px+xxx+yxy+zxz+fx同理得y方向：(v)t+vV=-py+xyx+yyy+zyz+fyZ方向：(w)t+wV=-pz+xzx+yzy+zzz+fz对于不可压缩流体：u=t=0对于牛顿流体：DuDt=f-p+2u将上述式子带入N-S方程，联立可得：DuDt=-p+2u+F其中，u=u,v,wT;F=Fx,Fy,FzT2 用有限体积法离散二维对流扩散方程xu+yv=xx+yy+S,将得出代数方程：app=anbnb+Su并确定式中系数。解：对扩散方程在上图控制容积内积分，有：

4、sVxudV+VyvdV=VxxdV+VyydV+VSdV由图可知，控制容积的边界面积（长度）Aw=Ae=y,An=As=x则上述方程式可以写为uAe-uAw+uAn-uAs=Axe-Axw+Axn-Axs+Sxy计算通过控制容积截面的场变量值或者导数值：穿过东侧边界的对流量：Ce=ueAee=ueAe(E+p)/2穿过东侧边界的扩散量：Ie=eAexe=eAeE-PxPE穿过西侧边界的对流量：Cw=uwAww=uwAw(P+w)/2穿过西侧边界的扩散量：Iw=wAwxw=wAwE-PxWP穿过北侧边界的对流量：Cn=vnAnn=vnAn(N+P)/2穿过北侧边界的扩散量：In=nAnyn=n

5、AnN-PyPN穿过南侧边界的对流量： Cs=vsAss=vsAs(P+S)/2穿过南侧边界的扩散量：Is=sAsys=sAsN-PySP源项：SV=Su+SPP此时，另F=uA或F=vA，D=Ay或D=Ax则：Fe=(u)eAe, Fw=(u)wAw, Fn=(u)nAn, Fn=(u)nAn,De=eAexPE, Dw=wAwxWP,Dn=nAnyPN,Ds=sAsySP带入通过各界面的扩散量及对流量，然后积分方程变为：Fe2E+P-Fw2P+W+Fn2N+P-FS2P+S=DeE-w-DwP-w+Dn N-P-Ds P-S+Su+SPp按照节点变量整理，最终可得：Dw+Fw2+De-Fe

6、2+Ds+Fs2+Dn-Fn2+Fe-Fw+Fn-Fs-SPP=Dw+Fw2W+De-Fe2E+Ds+Fs2S+Dn-Fn2N即app=aWW+aEE+aSS+aNN+Su式中，aW=Dw+Fw2, aE=De-Fe2, aS=Ds+Fs2, aN=Dn-Fn2,aP=aW+aE+aS+aN+F-SP ,F=Fe-Fw+Fn-Fs2设某场变量经过对流扩散过程从一维区域的x=0点输运到x=L点，输运方程为ddxu=ddxddx其中：流体密度=1.0kg/m3，L=2.0m,扩散系数=0.2kg/(m.s)边界条件x=0=0=1,x=L=L=0求（1）当流速u=0.1m/s，对对流项采用中心差分格

7、式将区域离散成十个节点时，用有限体积法求在区域内的分布并作图。（2）当流速u=2.5m/s，对对流项采用迎风格式将区域离散成二十个节点时，用有限体积法求在区域内的分布并作图。解：求解区域的离散网格系统如图所示：每一个控制容积的长度x=0.2m,此时F=u,D=/x,Fe=Fw=F,De=Dw=D,对所有控制容积都成立，中间节点29的离散方程为app=aWW+aEE，边界点1与边界点10需要特殊处理。对节点1所在的控制容积积分，有VxudV=VxxdV（此时点1为P）则(uA)e-uAw=Addxe-(Addx)w在西侧界面时，有w=A=1,(ddx)w=P-Ax/2则Fe2P+E-FAA=De

8、E-P-2DA(P-A)按照节点场变量整理，有De+Fe2+2DAP=0*W+De-Fe2E+(2DA+FA)A即 aPP=aww+aee+Su式中aw=0，aE=De-Fe2, aP=aW+aE+Fe-Fw-SP,SP=-2DA+Fw, Su=2DA+FAA当u=常数时，aP表达式中（Fe-Fw)=0。同理，对节点5所在的控制容积积分，得：FBB-Fw2P+E=2DBB-P-Dw(P-W)按节点场变量整理，有Dw-Fw2+2DBP=Dw+Fw2W+0*E+(2DB-FB)B即aPP=aww+aee+Su式中aW=Dw+Fw2, aE=0, aP=aw+aE+(Fe-Fw)-SPSP=-(2D

9、B-Fe), Su=2DB-FBB综上，我们得到下表所列的离散方程系数和等效源项计算公式节点aWaESPSu10D-F/2-(2D+F)(2D+F)A2到9D+F/2D-F/20010D+F/20-(2D-F)(2D-F) B第一种计算情况时，u=0.1ms，F=u=0.1,D=x=0.20.2=1,A=1,B=0,带入上表，可得离散方程系数节点aWaESuSPaP100.951.1-2.13.052到91.050.95002.0101.0500-2.13.15将方程写成增广矩阵形式（未知数为各节点处的值）3.05-1.05000000002.1-0.952-1.05000000000-0.9

10、52-1.05000000000-0.952-1.05000000000-0.952-1.05000000000-0.952-1.05000000000-0.952-1.05000000000-0.952-1.05000000000-0.952-1.05000000000-0.952.950用追赶法解此方程组，得：0.970940.8203480.6840990.5608260.4492930.3483820.2570820.1744770.0997390.03211912345678910=%图有问题（2）此时x=220=0.1m,当流动为正方向（从左往右）时，采用迎风格式即为：w=W；e=

11、P此时离散方程为FeP-FwW=DeE-P-Dw(P-W) 按照节点场变量整理为：Dw+De+FeP=Dw+FwW+Dee在边界点1处，对流项采用迎风格式：FeP-FAA=DeE-P-DA(P-A)同理，在边界点20处：FBP-FwW=DBB-P-Dw(P-W)按照节点场变量整理为aPP=aww+aee+Su的形式，系数值见下表节点aWaESuS_P10D（2D+F）A-（2D+F）219D+FD0020D+F02DB-2D表中u=0.1ms，F=u=2.5,D=x=0.20.1=2,A=1,B=0带入上表，得离散方程系数：节点aWaESuSP1026.5-6.52194.5200204.50

12、0-4写成方程组的形式，解得：各节点处的值分别为1=0.999999952 2=0.9999997953=0.999999443 4=0.9999986515=0.999996868 6=0.9999928567=0.99998383 8=0.9999635219=0.999917826 10=0.99981501311=0.999583683 12=0.99906319213=0.997892085 14=0.99525709415=0.989328366 16=0.97598872617=0.945974537 18=0.87844261319=0.726495782 20=0.38461

13、5414附：1: function x = zhuigan(b,a,c,f)%追赶法求解三对角方程组%b=diag(B);%a=diag(B,-1);%c=diag(B,1);%a,b,c，f,x都是列向量a=0,a'n,=size(b);alpha=zeros(n,1);y=zeros(n,1);x=zeros(n,1);for i=1 alpha(i)=c(i)/b(i);endfor i=2:n-1 alpha(i)=c(i)/(b(i)-a(i)*alpha(i-1);endfor i=1 y(i)=f(i)/b(i);endfor i=2:n y(i)=(f(i)-a(i)*

14、y(i-1)/(b(i)-a(i)*alpha(i-1);endfor i=n x(i)=y(i);endfor i=n-1:-1:1 x(i)=y(i)-alpha(i)*x(i+1);end2:function T = zhongxinchafen(L,rho,a,b,G,u,n)%一维情形下离散方程的形成%差分格式%A的第一列为aw，第二列为ae，第三列为Su，第四列为Sp，第五列为ap；%先输n值，n值为网格数量%看情况改变F与D的值deltax=L/n;x=deltax/2:deltax:L-deltax/2;D=G/deltax; F=rho*u;A(1,1)=0;A(1,2)=D

15、-F/2;A(1,3)=(2*D+F)*a;A(1,4)=-(2*D+F);A(1,5)=A(1,1)+A(1,2)-A(1,4);for i=2:n-1 A(i,1)=D+F/2; A(i,2)=D-F/2; A(i,3)=0; A(i,4)=0; A(i,5)=A(i,1)+A(i,2)-A(i,4);endA(n,1)=D+F/2;A(n,2)=0;A(n,3)=(2*D-F)*b;A(n,4)=-(2*D-F);A(n,5)=A(n,1)+A(n,2)-A(n,4);p=A(:,5);w=A(:,1);w=-w(w=0);e=A(:,2);e=-e(e=0);f=A(:,3);T=a,

16、zhuigan(p,w,e,f)',b;plot(0,x,L,T,'-',0,x,L,T,'*');xlabel('距离')ylabel('变量分布')end3:function T1 = yingfeng1(L,rho,a,b,G,u,n)%一维情形下离散方程的形成%迎风%u为流速%A的第一列为aW，第二列为aE，第三列为Su，第四列为Sp，第五列为apdeltax=L/n;x=deltax/2:deltax:L-deltax/2;D=G/deltax; F=rho*u;A(1,1)=0;%节点1的aWA(1,2)=D;%节点1的aEA(1,3)=(2*D+F)*a;%节点1的SuA(1,4)=-(2*D+F);%节点1的SpA(1,5)=A(1,1)+A(1,2)-A(1,4);for i=2:n-1 A(i,1)=D+F; A(i,2)=D; A(i,3)=0; A(i,4)=0; A(i,5)=A(i,1)+A(i,2)-A(i,4);endA(n,1)=D+F;%节点5的