有時候在NS方程中我們需要構造有精確解的例子去測試格式的精度，這時候往往都是先給出精確解的表達式，然后帶入原NS方程求出源項。如何求這個源項呢，對于可壓流NS方程來說手動求導數會很麻煩。所以，這里我們用matlab符號計算來求源項。

clear all syms x y u1 u2 u3 u4 pressure gamma txx txy tyx tyy mu Pr kTx kTy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%各變量定義%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %x x方向坐標 %y y方向坐標 %u1 密度 %u2 x方向動量 %u3 y方向動量 %u4 總能 %pressure 壓強 %gamma 氣體常數 %txx txy tyx tyy 粘性應力張量 %mu 粘性系數 %Pr Prandtl數 %kTx $\kappa T_x$ %kTy $\kappa T_y$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% u1 = sin(2*(x+y)) + 4; u2 = sin(2*(x+y)) + 4; u3 = sin(2*(x+y)) + 4; u4 = (sin(2*(x+y)) + 4)^2;pressure = (gamma - 1)*(u4 - 0.5*(u2^2+u3^2)/u1); txx = 2.0/3.0 * mu * (2*diff(u2/u1,x)-diff(u3/u1,y)); txy = mu*(diff(u2/u1,y)+diff(u3/u1,x)); tyx = mu*(diff(u3/u1,x)+diff(u2/u1,y)); tyy = 2.0/3.0*mu*(2*diff(u3/u1,y)-diff(u2/u1,x)); kT = mu * gamma / Pr *(u4/u1-0.5*(u2^2+u3^2)/(u1*u1)); kTx = diff(kT,x); kTy = diff(kT,y);f1x = u2; f1y = u3; g1x = 0; g1y = 0;f2x = u2^2/u1 + pressure; f2y = u2*u3/u1; g2x = txx; g2y = txy;f3x = u2*u3/u1; f3y = u3^2/u1 + pressure; g3x = txy; g3y = tyy;f4x = u2/u1*(u4+pressure); f4y = u3/u1*(u4+pressure); g4x = u2/u1*txx + u3/u1*txy + kTx; g4y = u2/u1*txy + u3/u1*tyy + kTy;d1 = diff(f1x,x) + diff(f1y,y) - diff(g1x,x) - diff(g1y,y); d2 = diff(f2x,x) + diff(f2y,y) - diff(g2x,x) - diff(g2y,y); d3 = diff(f3x,x) + diff(f3y,y) - diff(g3x,x) - diff(g3y,y); d4 = diff(f4x,x) + diff(f4y,y) - diff(g4x,x) - diff(g4y,y);d1 = simplify(d1) d2 = simplify(d2) d3 = simplify(d3) d4 = simplify(d4)

轉載于:https://www.cnblogs.com/yuehq/p/5411957.html

總結

以上是生活随笔為你收集整理的NS方程由精确解求源项matlab代码的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。