文章目錄

• • 繞柱平流問題
  • • 問題描述
    • 觀察量
    • 工具箱實現
  • 周期圓柱繞流問題
  • 固定時間間隔圓柱擾流

NS 方程是偏微分方程領域最重要的一個方程，沒有之一。可惜的是，MATLAB 的自帶的工具箱無法對它進行很好地求解，主要原因在于非線性項（對流項）無法處理，如果把它和右端載荷項放在一塊處理，顯然 MATLAB 工具箱是沒辦法支持這件事情的。

那么，對于不想寫太多代碼，又想用 MATLAB 這個工具求解 NS 方程，有什么好辦法呢？我給大家推薦 FEATool Multiphysics。這個工具夠傻瓜夠強大，能求解很多 CFD 的問題，能和 FEniCS 和 Openfoam 聯動，求足夠完整的文檔，以及完備的社區。我是推薦。下面舉兩個簡單的例子來說明怎么使用它求解 NS 方程。安裝就不用說，一搜就有了。

繞柱平流問題

問題描述

這是一個典型的 Benchmark 的問題了，也可以參考 DFG flow around cylinder benchmark 2D-1, laminar case Re=20。
考慮這樣一個區域（圖畫得很清楚了，不用多解釋）：

我們要求解的問題是：
$$\{\begin{array}{r}\rho \left(\frac{?\mathbf{u}}{?t}+(\mathbf{u}??)\mathbf{u}\right)???\left(\mu \left(?\mathbf{u}+?{\mathbf{u}}^T\right)\right)+p=\mathbf{F}\\ ??\mathbf{u}=0\end{array}$$
其中：
$\rho =1,\mu =0.001$。

左邊界條件：
$${\mathbf{u}}_{\text{inflow?}}=\mathbf{u}(0,y)=4U_{\max }{H}^{?2}(y(H?y),0)$$
這里的 $U_{\max }=0.3$。$H=0.41$ 表示上圖區域的高度。

在這條件下，$U_{\text{mean?}}=0.2$，$Re=\rho U_{\text{mean?}}D/\mu =20$，這就導致了層流。此時，變成了一個穩態問題。

右邊界條件：
$$\nu {?}_{\eta }u?p\eta =0$$
這里的 $\eta$ 是外法向。

其他的邊界是無滑移邊界條件（名詞看起來高大上，其實就是速度為 0）。

觀察量

我們定義阻力系數和升力系數如下：
$$c_d=\frac{2F_d}{\rho U_{\text{mean?}}^{2}D},c_l=\frac{2F_l}{\rho U_{\text{mean?}}^{2}D}$$
其中的阻力和升力定義為：
$$F_d={\int }_S\left(\mu \frac{?{\mathbf{u}}_{\tau }(t)}{?n}{n}_y?p{n}_x\right)dS,F_l=?{\int }_S\left(\mu \frac{?{\mathbf{u}}_{\tau }(t)}{?n}{n}_x+p{n}_y\right)dS$$
其中的 ${\mathbf{u}}_{\tau }$ 是切方向。

工具箱實現

Benchmark 1。
這個問題是 FEATool Multiphysics 的一個算例。打開工具箱，選擇圓柱繞流算例，就是這個例子，點擊 Run，它會自動地給你演示如何做。演示完也就做完了，你可以把變量導入到工作空間，也可以把整個過程生成代碼。看一下代碼：

%% FEATool Multiphysics Version 1.16, Build 22.08.241, Model M-Script file. %% Created on 18-Oct-2022 11:22:20 with MATLAB 9.9.0.1467703 (R2020b) PCWIN64.clc clear close all %% Starting new model. fea.sdim = { 'x', 'y' }; fea.geom = struct; fea = addphys( fea, @navierstokes, { 'u', 'v', 'p' } );%% Geometry operations. gobj = gobj_rectangle( 0, 1, 0, 1, 'R1'); fea = geom_add_gobj( fea, gobj ); gobj = gobj_rectangle( [0], [2.2], [0], [0.41], 'R1' ); fea.geom.objects{1} = gobj; gobj = gobj_ellipse( [ 0.4;0.2 ], 0.2, 0.1, 'E1'); fea = geom_add_gobj( fea, gobj ); gobj = gobj_ellipse( [0.2 0.2], [0.05], [0.05], 'E1' ); fea.geom.objects{2} = gobj; fea.geom = geom_apply_formula( fea.geom, 'R1-E1' );%% Grid operations. fea.grid = gridgen( fea, 'gridgen', 'default', 'hmax', 0.13, 'grading', 0.3 ); fea.grid = gridgen( fea, 'gridgen', 'default', 'hmax', 0.02, 'grading', 0.3 );%% Equation settings. fea.phys.ns.dvar = { 'u', 'v', 'p' }; fea.phys.ns.sfun = { 'sflag1', 'sflag1', 'sflag1' }; fea.phys.ns.eqn.coef = { 'rho_ns', 'rho', 'Density', { 'rho' };'miu_ns', 'mu', 'Viscosity', { 'miu' };'Fx_ns', 'F_x', 'Volume force in x-direction', { '0' };'Fy_ns', 'F_y', 'Volume force in y-direction', { '0' };'u0_ns', 'u_0', 'Initial condition for u', { '0' };'v0_ns', 'v_0', 'Initial condition for v', { '0' };'p0_ns', 'p_0', 'Initial condition for p', { '0' };'miuT_ns', [], 'Turbulent viscosity', { 0 } }; fea.phys.ns.eqn.sdiff = { ' - (miu_ns+miuT_ns)*(2*ux_x + uy_y + vx_y)', ' - (miu_ns+miuT_ns)*(vx_x + uy_x + 2*vy_y)', [] }; fea.phys.ns.eqn.sconv = { 'rho_ns*(u*ux_t + v*uy_t)', 'rho_ns*(u*vx_t + v*vy_t)', [] }; fea.phys.ns.eqn.seqn = { 'rho_ns*u'' - (miu_ns+miuT_ns)*(2*ux_x + uy_y + vx_y) + rho_ns*(u*ux_t + v*uy_t) + p_x = Fx_ns', 'rho_ns*v'' - (miu_ns+miuT_ns)*(vx_x + uy_x + 2*vy_y) + rho_ns*(u*vx_t + v*vy_t) + p_y = Fy_ns', 'ux_t + vy_t = 0' }; fea.phys.ns.eqn.vars = { 'Velocity field', 'sqrt(u^2+v^2)';'x-velocity', 'u';'y-velocity', 'v';'Pressure', 'p';'Vorticity', 'vx-uy';'Velocity field', { 'u', 'v' } }; fea.phys.ns.prop.isaxi = 0; fea.phys.ns.prop.artstab.id = 0; fea.phys.ns.prop.artstab.id_coef = 0.5; fea.phys.ns.prop.artstab.sd = 0; fea.phys.ns.prop.artstab.sd_coef = 0.25; fea.phys.ns.prop.artstab.ps = 1; fea.phys.ns.prop.artstab.ps_coef = 1; fea.phys.ns.prop.artstab.iupw = 0; fea.phys.ns.prop.turb.model = 'laminar'; fea.phys.ns.prop.turb.wallf = 1; fea.phys.ns.prop.turb.inlet = []; fea.phys.ns.prop.active = [ 1;1;1 ]; fea.phys.ns.prop.intb = 0;%% Constants and expressions. fea.expr = { 'h', '0.41';'diam', '0.1';'rho', '1';'miu', '0.001';'umax', '0.3';'umean', '2/3*umax';'fx', 'nx*p+miu*(-2*nx*ux-ny*(uy+vx))';'cd', '2*fx/(rho*umean^2*diam)' };%% Boundary settings. fea.phys.ns.bdr.sel = [ 1, 4, 1, 2, 1, 1, 1, 1 ]; fea.phys.ns.bdr.coef = { 'bw_ns', 'u = 0', 'Wall/no-slip', [], { 1, 1, 1, 1, 1, 1, 1, 1;1, 1, 1, 1, 1, 1, 1, 1;0, 0, 0, 0, 0, 0, 0, 0 }, [], { 0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0 };'bv_ns', 'u = u_0', 'Inlet/velocity', { 'u_0';'v_0';[] }, { 1, 1, 1, 1, 1, 1, 1, 1;1, 1, 1, 1, 1, 1, 1, 1;0, 0, 0, 0, 0, 0, 0, 0 }, [], { 0, 0, 0, '4*umax*y*(h-y)/h^2', 0, 0, 0, 0;0, 0, 0, '0', 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0 };'bn_ns', '-n.(-pI + mu(grad u+grad uT)) = 0', 'Neutral outflow/stress boundary', [], { 0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0 }, [], { 0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0 };'bp_ns', 'p = p_0', 'Outflow/pressure', { [];[];'p_0' }, { 0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0;1, 1, 1, 1, 1, 1, 1, 1 }, { { '-p*nx';'-p*ny';'0' };[] }, { 0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0;0, '0', 0, 0, 0, 0, 0, 0 };'bs_ns', 'n&.u = 0, -n.(-pI + mu(grad u+grad uT)).t = 0', 'Symmetry/slip', [], { 0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0 }, [], { 'solve_hook_bdrslip', 'solve_hook_bdrslip', 'solve_hook_bdrslip', 'solve_hook_bdrslip', 'solve_hook_bdrslip', 'solve_hook_bdrslip', 'solve_hook_bdrslip', 'solve_hook_bdrslip';0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0 } }; fea.phys.ns.bdr.coefi = { 'bcic_ns', 'n.(F_1-F_2) = 0, F_i=-p_iI + mu_i(grad u_i+grad u_iT)', 'Continuity', [], { 0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0 }, [], { 0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0 };'bcij_ns', 'u = u_0', 'Prescribed velocity', { 'u_0';'v_0';[] }, { 1, 1, 1, 1, 1, 1, 1, 1;1, 1, 1, 1, 1, 1, 1, 1;0, 0, 0, 0, 0, 0, 0, 0 }, [], { 0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0 };'bcir_ns', 'p = p_0', 'Prescribed pressure', { [];[];'p_0' }, { 0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0;1, 1, 1, 1, 1, 1, 1, 1 }, [], { 0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0;0, 0, 0, 0, 0, 0, 0, 0 } }; fea.phys.ns.bdr.vars = { 'Viscous force, x-component', '(miu_ns+miuT_ns)*(2*nx*ux+ny*(uy+vx))';'Viscous force, y-component', '(miu_ns+miuT_ns)*(nx*(vx+uy)+2*ny*vy)';'Total force, x-component', '-nx*p+(miu_ns+miuT_ns)*(2*nx*ux+ny*(uy+vx))';'Total force, y-component', '-ny*p+(miu_ns+miuT_ns)*(nx*(vx+uy)+2*ny*vy)' }; fea.phys.ns.prop.intb = 0;%% Solver call. fea = parsephys( fea ); fea = parseprob( fea );fea.sol.u = solvestat( fea, ...'iupw', [ 0, 0, 0 ], ...'linsolv', 'backslash', ...'icub', 'auto', ...'nlrlx', 1, ...'toldef', 1e-06, ...'tolchg', 1e-06, ...'reldef', 0, ...'relchg', 1, ...'maxnit', 20, ...'hook', [], ...'nproc', 2, ...'init', { 'u0_ns', 'v0_ns', 'p0_ns' }, ...'solcomp', [ 1; 1; 1 ] );%% Postprocessing. postplot( fea, ...'surfexpr', 'sqrt(u^2+v^2)', ...'title', 'surface: Velocity field', ...'solnum', 1 ); postplot( fea, ...'surfexpr', 'sqrt(u^2+v^2)', ...'isoexpr', 'p', ...'isolev', 10, ...'arrowexpr', { 'u', 'v' }, ...'arrowspacing', { 10, 10 }, ...'arrowscale', 1, ...'title', 'surface: Velocity field, contour: Pressure, arrow: Velocity field', ...'solnum', 1 ); postplot( fea, ...'surfexpr', 'sqrt(u^2+v^2)', ...'isoexpr', 'p', ...'isolev', 10, ...'arrowexpr', { 'u', 'v' }, ...'arrowspacing', { 10, 10 }, ...'arrowscale', 1, ...'title', 'surface: Velocity field, contour: Pressure, arrow: Velocity field', ...'solnum', 1 ); integral_value = intbdr( 'cd', fea, [5 6 7 8], 2, 1 );

需要指明的是，這份代碼的運行需要在工具箱打開，變量載入空間的時候才能 work。

周期圓柱繞流問題

這也是個 Benchmark2，可以參考 DFG flow around cylinder benchmark 2D-2, time-periodic case Re=100，此時它不再是工具箱的一個算例，我們需要做簡單的修改。幾何剖分沿用上面的兩個例子。

別的不需要變，我們唯一要改的只是 $U_{\max }$，從 0.3 改為 $U_{\max }=1.5$。

在這個情況下，$U_{\text{mean?}}=\frac{2}{3}?1.5=1.0$, $L=2?0.05=0.1$，$\mathrm{Re}=\frac{U_{\text{mean?}}L}{\nu }=\frac{1.0?0.1}{0.001}=100$。這時候算出來雷諾數是 100。這就形成了時間周期行為。

直接修改上述的代碼中的 $U_{\max }$，可以得到。

固定時間間隔圓柱擾流

考慮 Benchmark 3，DFG flow around cylinder benchmark 2D-3, fixed time interval (Re=100)。

事實上，依然還是修改 $U_{\max }=1.5\sin (\pi ?t/8)$。修改為

$$u_{\text{inflow?}}=6?\sin (\pi ?t/8)?(y?(0.41?y))/0.41^2$$

事實上，關于這個，工具箱封裝了一個內置函數 function [ fea, out ] = ex_navierstokes6( varargin ) 直接可以用。其中，可調參數可以參考：ex_navierstokes6 參數表。
通過閱讀代碼可以知道，通過修改 instatbc=false，就可以把 Benchmark 3 改成了 Benchmark 2，即周期圓柱繞流問題。即寫成：

clc clear close all tic [ fea, out ] = ex_navierstokes6( 'instatbc',false ); toc

最后 show 一個結果如下：

這里的總時間 $t=8$ ，時間步長 $dt=0.02$，網格尺寸使得空間自由度 9456。所消耗的時間 8 分鐘。

總結

以上是生活随笔為你收集整理的MATLAB 工具箱傻瓜式求解 NS（Navier Stoke）方程的全部內容，希望文章能夠幫你解決所遇到的問題。

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