一种基于界带有限元和拉格朗日坐标的流体仿真方法

摘要

本发明提出了一种基于界带有限元和拉格朗日坐标的不可压缩流体仿真分析方法，是将二维不可压缩流体的计算域Ω按照传统有限元网格剖分成Ne个单元，每个单元为Ωi；构造单元Ωi的位移插值场，根据位移插值场构造流体的动力微分方程，求解该动力微分方程得到流体的各种物理参数，从而进行流体的运动分析；其特征在于用界带有限单元法构造位移插值场；并基于拉格朗日坐标描述法得到流体的动力微分方程。本发明将拉格朗日坐标方法与界带有限元方法结合来解决不可压缩流体的运动仿真问题，目的是利用界带有限元精度高和拉格朗日坐标下边界处理方便，通用性好的优势，提高分析的计算效率和精度。

说明书

技术领域

本发明涉及流体仿真分析技术，构建了一种基于界带有限元和拉 格朗日坐标的流体仿真方法。

背景技术

流体的仿真分析广泛应用于各个不同工程领域，如水利工程、航 空飞行、高速列车等。流体的仿真分析，从理论上看主要有欧拉坐标 和拉格朗日坐标两种，从数值分析手段来说，主要有限差分法和有限 单元法两种。其中有限差分方法使用规则的矩形网格，因此对于在不 规则体中的流体仿真分析需要加密网格，这导致计算量大幅增加。有 限单元法在固体结构仿真分析中得到广泛运用，该方法的特点是适用 于不规则问题，但是需要变分原理。

目前，流体的仿真分析主流是基于欧拉坐标，在欧拉坐标系可以 导出不可压缩流体的运动微分方程，即著名的纳维-斯托克斯方程，对 于该方程的计算最常用的方法是有限差分法。但是基于欧拉坐标的流 体仿真分析在处理自由面流体问题时，分析困难，计算格式复杂。由 于自由面随时间不断变化，其计算区域也不再规则，因此用有限差分 分析困难，而且基于欧拉坐标的流体方程，很难建立变分原理，从而 导致有限元分析也很困难。目前仅能对某些特殊的流体运动问题，建 立变分，使用有限元分析。在分析不可压缩流体时，采用普通有限元， 还存在体积闭锁等问题，影响计算精度。

如果运用拉格朗日坐标，则可以很容易得到流体的动能表达式和 势能表达式，并可以利用哈密尔顿变分原理导出拉格朗日坐标系的流 体动力微分方程。根据动能表达式和势能表达式，还可以利用有限元 进行分析。但是在拉格朗日坐标下，以流体中质点的位移为基本未知 量，以位移为基本未知量建立的传统有限单元不能完全满足不可压缩 条件，计算精度不好。实际上在固体不可压缩材料的仿真分析中，有 学者提出界带有限元方法，通过引入流函数来代替位移作为未知量， 并利用界带单元进行分析，所构造的单元完全满足不可压缩条件。界 带单元可以在传统有限元网格上进行分析，因此可以利用现有有限元 网格剖分技术，而且具有计算精度高、所需的自由度小等优点。目前 还没有关于流体仿真分析中采用界带有限元方法的应用。

发明内容

本发明针对现有技术不足，提出一种分析二维不可压缩流体的仿真 分析方法，该方法将拉格朗日坐标方法与界带有限元方法结合来解决 不可压缩流体的运动仿真问题，目的是利用界带有限元精度高和拉格 朗日坐标下边界处理方便，通用性好的优势，提高分析的计算效率和 精度。

为此，本发明提供了一种基于界带有限元和拉格朗日坐标的流体 仿真方法，将二维不可压缩流体的计算域Ω按照传统有限元网格剖分成 N_e个单元，每个单元为Ω_i；构造单元Ω_i的位移插值场，根据位移 插值场构造流体的动力微分方程，求解该动力微分方程得到流体的各 种物理参数，从而进行流体的运动分析；本发明用界带有限单元法构 造位移插值场；并基于拉格朗日坐标描述法得到流体的动力微分方程； 具体方法如下：

(a)将二维不可压缩流体的计算域Ω采用传统有限元网格剖分成 N_e个单元，每个单元为Ω_i；

(b)在单元Ω_i上建立流函数ψ(x,y)的插值场时，以单元Ω_i为本体， 将Ω_i周边的单元视为Ω_i的界带，将这些单元的节点合作进行插值，构造 单元Ω_i上的插值函数ψ(x,y)；

(c)将步骤(b)中所述的流函数表达式求偏导，得到流体中在坐标 (x,y)处质点的位移表达式；

(d)根据上一步的位移场表达式，得到每个单元上流体的质量矩阵 M_i，并把所有单元的质量矩阵计算出来后，通过累加得到总体质量矩 阵M；

(e)根据步骤(c)的位移场表达式，得到流体的刚度矩阵：

K(ψ)＝K_l+K_n(ψ)

其中K是刚度矩阵，ψ是所有节点的流函数值所组成的向量，K_l是线性 刚度矩阵，K_n是非线性刚度矩阵；

(f)流体的边界包括两种，一种为自由面Γ_f，一种为不可穿过边界Γ_n， 根据不可穿过边界Γ_n上所有有限元单元节点的编号，把刚度矩阵和质量 矩阵中相应编号的行和列划去，得到描述流体运动的非线性微分方程：

$M\stackrel{··}{\psi }+K\left(\psi \right)\psi =0$

(g)利用非线性微分方程求解软件求解上述非线性微分方程，得到 不同时间点上的流函数向量ψ；

(h)根据流函数向量ψ，得到流体中各个质点在不同时间上的位移； 根据流体中各质点在不同时间点的位移，利用有限元后处理程序，可 以得到流体的动态仿真图；

(i)在步骤(g)中，刚度矩阵包含非线性和线性两个部分，如果非线 性因素较小而无需考虑时，可以得到线性的动力微分方程此时利用微分方程求解软件求解该线性微分方程，可以得到线性情况 下的流函数；

(j)如果要分析流体的振动模态和频率，利用特征值求解软件求解 下面的特征值方程

K_lψ＝ω²Mψ

其中，ω即为水的振动频率。

在步骤(b)中，界带单元上流函数的插值函数具体格式为：

ψ(x,y)＝N^Tψ_i,N^T＝p^T(x,y)P^-1

其中

p^T(x,y)＝(1,x,y,x²,xy,y²,…)

P^T＝[p(x₁,y₁),p(x₂,y₂),…,p(x_N,y_N)]

ψ^T＝(ψ(x₁,y₁),ψ(x₂,y₂),…,ψ(x_N,y_N))

N^T是形函数向量。(x_j,y_j)是界带单元的节点坐标，共N个，包含单元本 体节点和本体单元的周边单元节点坐标，N是界带单元的节点总数。 p^T(x,y)是插值多项式，其项数与N相同。

在步骤(c)中，流体中在坐标(x,y)处质点的位移为：

$u(x,y)=\frac{\text{∂ψ}}{\partial y}=N_y^T{\psi }_i,v(x,y)=-\frac{\partial \psi }{\partial x}=-N_x{\psi }_i$

其中N_x和N_y分别表示形函数向量N对坐标求偏导。

在步骤(d)中，对单元进行积分以计算质量矩阵，积分方案选择数 值积分，也即

$M_i=\underset{n=1}{\overset{n}{\Sigma }}{w}_n\rho (x_n,y_n)(N_x{N}_x^T+N_y{N}_y^T)$

其中，(x_n,y_n)为单元内部的数值积分点，w_n为积分点的权函数，n是积 分点数目。

在步骤(e)中，根据步骤(c)的位移场表达式，得到流体自由面上Γ_f的 位移u和v，在重力作用下，其刚度矩阵的计算表达式为：

K(ψ)＝K_l+K_n(ψ)

$K_1={\int }_{{\Gamma }_f}\rho (x,y)g{N}_x{N}_x^T\mathrm{dl},K_n\left(\psi \right)={\int }_{{\Gamma }_f}\rho (x,y){\mathrm{gN}}_x{N}_x^T\left[N_{\mathrm{xy}}^T\psi \right]\mathrm{dl}$

其中K是刚度矩阵，ψ是所有节点的流函数值所组成的向量，K_l是线性 刚度矩阵，K_n是非线性刚度矩阵，ρ是流体的密度，g是重力加速度， Γ_f是流体的自由面。

在步骤(a)中，对流体进行有限元网格剖分，可以采用现有有限元 软件，如Ansys软件。

在步骤(g)和(i)中，求解微分方程的软件可以采用现有的软件包， 如大连理工大学开发的SIPESC软件，SIPESC软件是大连理工大学工 程力学系开发的工程计算分析软件平台。其功能包括集成开发环境、 面向系统集成的活动流程图定制、工程数据库管理系统、开放式结构 有限元分析系统、集成优化计算系统等，其中有限元分析系统中集成 了方程求解模块、有限元后处理模块等，可以对本发明中微分方程进 行求解。

在步骤(h)中，后处理显示流体的动态运动的软件，可以采用现有 的后处理软件，如大连理工大学开发的SIPESC.POST、SIPESC.Jifex 软件，SIPESC.POST和SIPESC.Jifex均属于SIPESC软件中有限元系 统模块，其功能为有限元结果后处理，可视化动画显示等。

在步骤(j)中，求解特征值方程，可以采用现有软件，如大连理工大 学开发的SIPESC软件。

本发明的有益积极效果：

1.本发明基于拉格朗日坐标下的流体动力方程，较传统欧拉坐标 分析而言，物理意义更为明显，处理自由边界条件更为方便和精确。 由于利用了有限元方法，与有限差分方法相比，在处理复杂流体边界 时，更加方便和精确。

2.本发明基于界带有限元，建立流函数形成的位移场，该位移场 满足不可压缩条件，避免了传统有限单元的体积闭锁等缺陷，网格可 以利用传统有限单元网格，与传统单元不同在于，界带单元把传统单 元作为本体单元，比利用本体单元周边的单元节点，从而构造出高阶 位移场，因此所需的计算自由度要少于传统有限单元，计算效率和计 算精度更高。

3.本发明可以使用于不同工况分析中，包括线性流体模态分析、 线性流体运动仿真分析和非线性流体运动仿真分析。目前已用于模拟 非线性物理学中的孤立子传播、水池中水的晃动等问题的分析。

附图说明

附图1是本发明实施流程图。

附图2是三角形界带单元。图中，i,j,k,l,m,n分别代表节点的编 号。

附图3是四边形界带单元，图中1～12分别代表节点编号。

附图4是某水池初始时刻二维垂直面形状。

附图5是某水池20秒时二维垂直面形状。

附图6是某水池40秒时二维垂直面形状。

具体实施方式

下面结合图1具体阐述本发明的具体实施方式:

具体实施方式一：本实施方式用于分析非线性不可压缩流体的运 动仿真问题。

(a)将流体的计算域Ω采用传统有限元网格剖分成N_e个单元，每个 单元为Ω_i；

(b)在单元Ω_i上建立流函数ψ(x,y)的插值场。以单元Ω_i为本体，将 单元Ω_i周边的单元视为该单元的界带，将这些单元的节点合作进行插值， 构造单元Ω_i上的插值函数，插值函数的具体格式为：

ψ(x,y)＝N^Tψ_i,N^T＝p^T(x,y)P^-1

其中

p^T(x,y)＝(1,x,y,x²,xy,y²,…)

P^T＝[p(x₁,y₁),p(x₂,y₂),…,p(x_N,y_N)]

ψ^T＝(ψ(x₁,y₁),ψ(x₂,y₂),…,ψ(x_N,y_N))

N^T是形函数向量。(x_j,y_j)是界带单元的节点坐标，共N个，包含单元本 体节点和本体单元的周边单元节点坐标，N是界带单元的节点总数。 p^T(x,y)是插值多项式，其项数与N相同。

以图2为例，进一步阐述界带单元插值函数的构造，图2给出三 角形界带单元，该界带单元由四个普通的线性三节点三角形单元组成， 其中阴影部分为界带单元的单元本体，该单元本体附近还有3个周边 单元。当需要计算单元本体的流函数时，取单元本体[i,j,k]的三个节点， 并联合三个周边单元，共六个节点，共同插值，此时界带单元的插值 函数为：

ψ(x,y)＝N^Tψ_i,N^T＝p^T(x,y)P^-1

p^T(x,y)＝(1,x,y,x²,xy,y²)

P^T＝[p(x_i,y_i),p(x_j,y_j),p(x_k,y_k),p(x_l,y_l),p(x_m,y_m),p(x_n,y_n)]

ψ^T＝(ψ(x_i,y_i),ψ(x_j,y_j),ψ(x_k,y_k),ψ(x_l,y_l),ψ(x_m,y_m),ψ(x_n,y_n))

根据图2所示三角形界带单元可以看到，界带单元可以利用传统单元 的网格剖分，构造出高价插值场，计算精度要好于传统单元的精度。 图2所示三角形界带单元仅为本发明的实施例，并非限制本发明的保 护范围。对于任何传统有限单元，以其中某个单元为本体，并利用单 元本体的周边单元与该单元本体共同构造界带单元的插值函数，均属 于界带单元函数，故凡运用界带单元思想构造的其他单元，均包含在 本发明的保护范围内。

(c)将步骤(b)中所述的流函数表达式求偏导，得到流体中在坐标 (x,y)处质点的位移：

$u(x,y)=\frac{\text{∂ψ}}{\partial y}=N_y^T{\psi }_i,v(x,y)=-\frac{\partial \psi }{\partial x}=N_x{\psi }_i$

其中N_x和N_y分别表示形函数向量N对坐标求偏导。

(d)根据上一步的位移场表达式，得到每个单元上流体的质量矩阵：

$M_i={\int }_{{\Omega }_i}\rho (x,y)(N_x{N}_x^T+N_y{N}_y^T)\mathrm{dxdy}$

其中M_i是单元质量矩阵，ρ是流体的密度。把所有单元的质量矩阵计算 出来后，通过叠加得到总体质量矩阵M。在该步骤中，对单元进行积 分以计算质量矩阵，积分方案选择数值积分，也即

$M_i=\underset{n=1}{\overset{n}{\Sigma }}{w}_n\rho (x_n,y_n)(N_x{N}_x^T+N_y{N}_y^T)$

其中，(x_n,y_n)为单元内部的数值积分点，w_n为积分点的权函数，n是积分 点数目。

(e)根据步骤(c)的位移场表达式，得到流体自由面上的位移表达式 Γ_f，在重力作用下，其刚度矩阵的计算表达式为：

K(ψ)＝K_l+K_n(ψ)

$K_1={\int }_{{\Gamma }_f}\rho (x,y)g{N}_x{N}_x^T\mathrm{dl},K_n\left(\psi \right)={\int }_{{\Gamma }_f}\rho (x,y){\mathrm{gN}}_x{N}_x^T\left[N_{\mathrm{xy}}^T\psi \right]\mathrm{dl}$

其中K是刚度矩阵，ψ是所有节点的流函数值所组成的向量，K_l是线性 刚度矩阵，K_n是非线性刚度矩阵，ρ是流体的密度，g是重力加速度。 Γ_f是流体的自由面。在该步骤中，刚度矩阵的积分只在流体的自由面上 进行，计算时仍然采用数值积分，即

$\left(\begin{array}{cc}{K}_1=\underset{j=1}{\overset{J}{\Sigma }}{h}_j\rho (x_j,y_j)g{N}_x{N}_x^T& K_n\left(\psi \right)=\underset{j=1}{\overset{J}{\Sigma }}{h}_j\rho (x_j,y_j)g{N}_x{N}_x^T\left[N_{\mathrm{xy}}^T\psi \right]\end{array}\right)$

其中，(x_j,y_j)为自由面上的数值积分点，h_j为积分点的权函数，J是积分 点数目。

(f)流体的边界包括两种，一种为自由面Γ_f，一种为Γ_n，表示流体 不可穿过该边界。找出在边界Γ_n上所有节点的编号，把刚度矩阵和质量 矩阵中相应的行和列划去，得到描述流体运动的非线性微分方程：

$M\stackrel{··}{\psi }+K\left(\psi \right)\psi =0$

(g)利用非线性微分方程求解软件求解上述非线性微分方程，得到 不同时间点上的流函数向量ψ，求解微分方程的软件可以采用现有的软 件包，如大连理工大学开发的SIPESC软件；

(h)根据流函数向量ψ和步骤(c)中的位移表达式，得到流体中各个 质点在不同时间上的位移。根据流体中各质点在不同时间点的位移， 利用有限元后处理程序(如大连理工大学开发的SIPESC.POST、Jifex 等软件)，可以得到流体的动态仿真图。

具体实施方式二：本实施方案用于线性不可压缩流体的仿真分析， 与具体实施方式一不同点在于步骤(e)，此时不考虑非线性刚度矩阵部 分，可以得到线性的动力微分方程此时利用微分方程求解 软件求解该线性微分方程，可以得到线性情况下的个质点在不同时刻 的流函数，其余步骤相同。

具体实施方式三：本实施方式用于分析不可压缩流体振动模态， 与具体实施方式一和二的不同点在于步骤(e)、(f)、(g)和(h)，此时不考 虑非线性刚度矩阵部分，得到特征值问题：

K_lψ＝ω²Mψ

其中，ω为流体的振动频率。求解特征值方程，可以采用现有软件， 如大连理工大学开发的SIPESC软件。根据流函数模态ψ和具体实施方 式一中步骤(c)给出的位移表达式，得到流体中各个质点的位移。利用 有限元后处理程序(如大连理工大学开发的SIPESC.POST、Jifex等软 件)，可以得到流体的振动模态。

仿真实例：利用本发明方法，对某矩形水池中孤立波传播进行仿真。 该孤立水波从0秒时分离为两个孤立波，分别向两边传播，在60秒时， 孤立波碰到水池壁而反射，反射产生了部分尾波。首先，采用矩形单 元对矩形水池网格剖分，网格剖分如图4所示，图4中的水面是初始 时刻水波的形状。其次，在每个单元上建立流函数的插值函数，界带 单元采用如图3所示的四边形界带单元，插值节点包括单元本体和周 边单元的节点，其中单元本体的节点编号分别为：[3,4,5,6]，周边单元 的节点编分别为：[1,2,7,8,9,10,11,12]，共12个节点，其插值形函数表 达式为：

ψ(x,y)＝N^Tψ_i,N^T＝p^T(x,y)P^-1

p^T(x,y)＝(1,x,y,x²,xy,y²,x³,x²y,xy²,y³,x³y,xy³)

P^T＝[p(x₁,y₁),p(x₂,y₂),…,p(x₁₂,y₁₂)]

ψ^T＝(ψ(x₁,y₁),ψ(x₂,y₂),…,ψ(x₁₂,y₁₂))

根据单元插值函数，形成描述水波运动的控制方程求解该微分方程得到各时刻的流函数向量ψ，根据ψ计算得到各个时刻 水中质点的位移，可以利用后处理软件得到不同时刻水波变形图。

图5和6分别给出在20秒和40秒时的水波形状，由图可见，本 发明方法采用了拉格朗日坐标，可以精确跟踪流体中各个粒子的位移， 对于自由水面无需做特殊处理，因此可很方便地描述自由水面的不断 变化，这是传统采用欧拉坐标难以实现的。因为采用界带有限元，可 以精确地处理不可压缩性，计算精度高，没有体积闭锁，这是传统有 限单元难以实现的，且在边界处只需要对刚度矩阵和质量矩阵化行化 列处理，实际操作方便。仿真试验表明，本发明方法可以很方便和准 确地模拟水中孤立波的传播过程。