Step4：ナビエストークス方程式を解く
いよいよキャビティ流れについてナビエストークス方程式の数値計算を行います。
まずは初期状態と境界条件を設定します。
【初期状態】
- \(u, v, p = 0\)
【境界条件】
- \(u=1\) at \(y=2\)
- \(u, v=0\) その他
- \(\frac{\partial p}{\partial y}=0\) at \(y=0\)
- \(p=0\) at \(y=2\)
- \(\frac{\partial p}{\partial x}=0\) at \(x=0,2\)
境界条件はFields.hで定義し、Fielsds.cppで処理内容を書きます。
Fields.h（抜粋）
1 2 3 4 | //境界条件 void fixedValueBoundary(Fields::vectorField2d&, std::string&, double&); void specialPBoundary(Fields::vectorField2d&, std::string&); void zeroGradient(Fields::vectorField2d&, std::string&); |
Fielsds.cpp（抜粋）
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 | void Fields::fixedValueBoundary(Fields::vectorField2d& vec, std::string& name, double& pvalue) { if (name == "INLET") { for (int i = 0; i < vec.size(); i++) { int j = vec.size() - 1; vec[i][j].value = pvalue; } } else if (name == "OUTLET") { for (int i = 0; i < vec.size(); i++) { int j = 0; vec[i][j].value = pvalue; } } else if (name == "RIGHT") { for (int j = 1; j < vec.size() - 1; j++) { int i = vec.size() - 1; vec[i][j].value = pvalue; } } else if (name == "LEFT") { for (int j = 1; j < vec.size() - 1; j++) { int i = 0; vec[i][j].value = pvalue; } } else { std::cout << "The boundary name " << name << "is not found in the geometry" << std::endl; } } void Fields::zeroGradient(Fields::vectorField2d& vec, std::string& name) { if (name == "INLET") { for (int i = 0; i < vec.size(); i++) { int j = vec.size() - 1; vec[i][j].value = vec[i][j - 1].value; } } else if (name == "OUTLET") { for (int i = 0; i < vec.size(); i++) { int j = 0; vec[i][j].value = vec[i][j + 1].value; } } else if (name == "RIGHT") { for (int j = 1; j < vec.size() - 1; j++) { int i = vec.size() - 1; vec[i][j].value = vec[i - 1][j].value; } } else if (name == "LEFT") { for (int j = 1; j < vec.size() - 1; j++) { int i = 0; vec[i][j].value = vec[i + 1][j].value; } } else { std::cout << "The boundary name " << name << "is not found in the geometry" << std::endl; } } void Fields::specialPBoundary(Fields::vectorField2d& vec, std::string& name) { if (name == "RIGHT") { for (int j = 1; j < vec.size() - 1; j++) { int i = vec.size() - 1; vec[i][j].value = vec[i - 1][j].y; } } } |
境界条件を関数でまとめておきました。
Step5：残差の計算
ポアソン方程式を求める際に解の収束として、
\begin{align*}
\left|\frac{P_{n^{\prime}+1}-P_{n^{\prime}}}{P_{n^{\prime}}}\right|\tag{19}
\end{align*}
\left|\frac{P_{n^{\prime}+1}-P_{n^{\prime}}}{P_{n^{\prime}}}\right|\tag{19}
\end{align*}
(19)が収束することを確認して次のタイムステップに進むようにします。
収束判定は(19)が$10^{-3}$より小さくなるか50回ループとしています。
では、main.cppにポアソン方程式の収束させるループを追加します。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 | for (int nit = 0; nit < 50; nit++) { P = Pnew; b = -1.0 * (fd::gradXCDS(U) && fd::gradXCDS(U)) - 2.0 * (fd::gradYCDS(U) && fd::gradXCDS(V)) - (fd::gradYCDS(V) && fd::gradYCDS(V)); Pnew = (fd::poissonX(P)) + (fd::poissonY(P)) - (bSpacingRho * b); // P boundary condition FieldsOperations.fixedValueBoundary(Pnew, inlet, leftboundaryValue); FieldsOperations.zeroGradient(Pnew, outlet); FieldsOperations.zeroGradient(Pnew, right); FieldsOperations.zeroGradient(Pnew, left); double l1normal = FieldsOperations.L1norm(Pnew, P); //cout << "L1 norm Pressure Eqn Residual = "<< nit << " == " << l1normal << endl; if (l1normal < targetresidual) { std::cout << "pressure eqn is converged : " << nit << std::endl; std::cout << "L1 norm Pressure Eqn Residual = " << l1normal << std::endl; break; } } // Pressure Eqn end loop double l1normal = FieldsOperations.L1norm(Pnew, P); std::cout << "L1 norm Pressure Eqn Residual = " << l1normal << std::endl; |
ループを抜けたらナビエストークス方程式を解きます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | // du/dt+udu/dx+vdu/dy= -1/ρdp/dx + nu ∇^2 u // dv/dt+udv/dx+vdv/dy= -1/ρdp/dy + nu ∇^2 v Unew = U - dt * (U && fd::gradX(U)) - dt * (V && fd::gradY(U)) - dt * invrho * fd::gradXCDS(P) + dt * visc * (fd::laplacian2d(U)); Vnew = V - dt * (U && fd::gradX(V)) - dt * (V && fd::gradY(V)) - dt * invrho * fd::gradYCDS(P) + dt * visc * (fd::laplacian2d(V)); // UV boundary condition FieldsOperations.fixedValueBoundary(Unew, inlet, ulid); FieldsOperations.fixedValueBoundary(Unew, outlet, zerov); FieldsOperations.fixedValueBoundary(Unew, right, zerov); FieldsOperations.fixedValueBoundary(Unew, left, zerov); FieldsOperations.fixedValueBoundary(Vnew, inlet, zerov); FieldsOperations.fixedValueBoundary(Vnew, outlet, zerov); FieldsOperations.fixedValueBoundary(Vnew, right, zerov); FieldsOperations.fixedValueBoundary(Vnew, left, zerov); |
これで大まかな流れはできましたので、最後に全体のファイルをお見せして終了とします。
Follow @t_kun_kamakiri 