三维粘性流的数值模拟

12001.2.8三维粘性流的数值模拟1.控制方程三维非定常可压缩NS方程在直角坐标下的形式为zhygxfzhygxftqvvv（1）这里TE,w,v,u,qTzzzyzxzzzyzxvTyzyyyxyzyyyxvTxzxyxxxzxyxxvT2T2T2zTkwvu,,,,0hyTkwvu,,,,0gxTkwvu,,,,0fw)pE(,pw,vw,uw,whv)pE(,vw,pv,uv,vgu)pE(,uw,uv,pu,uf式中222wvu21p11E（2）为总能量，粘性系数由Sutherland公式给出CTCTTT05.100（3）这里的smkg107894.150，K15.288T0，K4.110C。2而Prckp为热传导系数，Pr是Prandtl数。粘性应力张量ywzv,zw34yv32xu32xwzu,zw32yv34xu32xvyu,zw32yv32xu34zyyzzzzxxzyyyxxyxx另外还要给出状态方程RTp（5）2.无量纲化取来流的密度、速度VV、温度T、黏性系数和飞行器的特征长度L为特征量，定义无量纲量22,,,,,p,u,v,w,E,T,tVxxzptxyzLLLLVuvwETVVVVT利用这些无量纲量，可将方程组（1）无量纲化，得到zhygxfRe1zhygxftqvvv（6）这里，在省略无量纲量上面的“－”后，q，无粘通量f、g、h，总能量E和粘性应力张量等表达式的形式不变，定压比热TMV11Tc11Tp1TpcccRRcc222vpppp从而粘性通量无量纲化后（省略无量纲量上面的“－”）成为3T2zzzyzxzzzyzxvT2yzyyyxyzyyyxvT2xzxyxxxzxyxxvzTMPr1wvu,,,,0hyTMPr1wvu,,,,0gxTMPr1wvu,,,,0fSutherland公式本身就是无量纲形式的，而状态方程的无量纲形式为TM1TVcTVpTVRTVRTVpp2222222省略无量纲量上面的“－”，就是TM1p2（7）以上的推导用以到下列关系式：Mayer公式Rccvp，音速pc，来流马赫数cVM，以及来流的状态方程RTp。现在来改写粘性通量。为简化下面的推导，记zx,yx,xx321，v3v2v1hf,gf,ff定义向量Tp,w,v,u,q~和矩阵Pr1wvu34cPr110000000000034000000M21140u320w000000000003200000000M,00u32v000000000000320000000M3121Pr1wv34ucPr110000003400000000000M2220v32w0000000320000000000000M,00uv32000000000320000000000M3221Pr1w34vucPr110340000000000000000M23350vw320000320000000000000000M,0u0w32000032000000000000000M3231则方程组（6）中的粘性通量可写成3,2,1ixq~Mf31jjiji（8）3.坐标变换通过坐标变换,,zz,,,yy,,,xx（9）可将方程组（6）变换到计算域，形式为vvvHGFRe1HGFtQ（10）式中Te,w,v,u,JJqQTTTTV)pe(,pzwV,pyvV,pxuV,VJhzgyfxJGU)pe(,pzwU,pyvU,pxuU,UJhzgyfxJF6TTW)pe(,pzwW,pyvW,pxuW,WJhzgyfxJH,U,V,WTvvvvTvvvvTvvvvFJfghuvwxyzxyzGJfghuvwxyzxyzHJfghuvxyzxyzw这里yyxxJ1z,zzxxJ1y,zzyyJ1xyyxxJ1z,zzxxJ1y,zzyyJ1xyyxxJ1z,zzxxJ1y,zzyyJ1x（11）而zzzyyyxxx,,Dz,y,xDJ是坐标变换的Jacobi行列式。方程组（10）的推导过程如下。用坐标变换的Jacobi行列式J乘7以方程组（6）的两边，得zhygxfJRe1zhygxfJtqJvvv（12）利用求导的链式规则，对上式左边的无粘通量，有hzJzJzJhzgyfxJgyJyJyJhzgyfxJfxJxJxJhzgyfxJhzJgyJfxJhzgyfxJhzJgyJfxJhzgyfxJhzJgyJfxJhzgyfxJzhygxfJzhygxfJzhygxfJzhzhzhJygygygJxfxfxfJzhygxfJ但是由（11）式80zyzyzyzyzyzyzyzyzyzyzyzyzyzyzyzyzyzyxJxJxJ222222222222同理0xJxJxJ0xJxJxJ就有HGFhzgyfxJhzgyfxJhzgyfxJzhygxfJ同样地，对（12）式右边的粘性通量，也有vvvvvvHGFzhygxfJ注意到坐标变换（9）与时间无关，J可直接进入（12）式左边的时间导数项，就得到方程组（10）。94.粘性通量的简化记321,,，v3v2v1HF,GF,FF则由（8）式，粘性通量又可写成333333111111jiiiimmnmnmmnmnjmmnmjnqqFJfJMJMxxxxx31jj31m31nmnnjmi31jj31m31nmnnjmiqqq~MxxJq~MxxJ3,2,1iQJ1qq~MxxJ31jj31m31nmnnjmi（13）定义矩阵q~qM，则1Mqq~，可用消去法求得。先求得11wvuwvu21000w000v000u00001p,w,v,u,E,w,v,u,q~qM222再求逆矩阵矩阵中的最后一行是由（2）式（已无量纲化）得到的。即11pE,wwE,vvE,uuE,wvu21E22210