MATLAB 工具箱傻瓜式求解 NS（Navier Stoke）方程

Posted 2023-02-01 陆嵩

文章目录

• • 绕柱平流问题
  • • 问题描述
    • 观察量
    • 工具箱实现
  • 周期圆柱绕流问题
  • 固定时间间隔圆柱扰流

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{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}\right)-\nabla \cdot \left(\mu \left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right)\right)+p=\mathbf{F}\\ \nabla \cdot \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 {\partial }_{\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{\partial {\mathbf{u}}_{\tau }(t)}{\partial n}{n}_y-p{n}_x\right)dS,F_l=-{\int }_S\left(\mu \frac{\partial {\mathbf{u}}_{\tau }(t)}{\partial 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.objects1 = 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.objects2 = 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}\cdot 1.5=1.0$, $L=2\cdot 0.05=$

多目标优化求解基于matlab非支配排序灰狼优化（NS-GWO）算法求解多目标优化问题含Matlab源码 2015期

MATLAB标定工具箱 相机标定

如何用MATLAB计算矩阵的行列式

文件上传工具类——傻瓜式上传文件

opencv级联分类器快速训练工具傻瓜式训练软件教程