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

Abstract

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

Description

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

技术领域

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

背景技术

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

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

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

发明内容

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

为此，本发明提供了一种基于界带有限元和拉格朗日坐标的流体仿真方法，将二维不可压缩流体的计算域Ω按照传统有限元网格剖分成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{\·\·}{\mathrm{\ψ}}+K\left(\mathrm{\ψ}\right)\mathrm{\ψ}=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{\∂\ψ}}{\∂y}=N_y^T{\mathrm{\ψ}}_i,v(x,y)=-\frac{\∂\mathrm{\ψ}}{\∂x}=-N_x{\mathrm{\ψ}}_i$

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

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

$M_i=\underset{n=1}{\overset{n}{\mathrm{\Σ}}}{w}_n\mathrm{\ρ}(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={\∫}_{{\mathrm{\Γ}}_f}\mathrm{\ρ}(x,y)g{N}_x{N}_x^T\mathrm{dl},K_n\left(\mathrm{\ψ}\right)={\∫}_{{\mathrm{\Γ}}_f}\mathrm{\ρ}(x,y){\mathrm{gN}}_x{N}_x^T\left[N_{\mathrm{xy}}^T\mathrm{\ψ}\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{\∂\ψ}}{\∂y}=N_y^T{\mathrm{\ψ}}_i,v(x,y)=-\frac{\∂\mathrm{\ψ}}{\∂x}=N_x{\mathrm{\ψ}}_i$

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

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

$M_i={\∫}_{{\mathrm{\Ω}}_i}\mathrm{\ρ}(x,y)(N_x{N}_x^T+N_y{N}_y^T)\mathrm{dxdy}$

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

$M_i=\underset{n=1}{\overset{n}{\mathrm{\Σ}}}{w}_n\mathrm{\ρ}(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={\∫}_{{\mathrm{\Γ}}_f}\mathrm{\ρ}(x,y)g{N}_x{N}_x^T\mathrm{dl},K_n\left(\mathrm{\ψ}\right)={\∫}_{{\mathrm{\Γ}}_f}\mathrm{\ρ}(x,y){\mathrm{gN}}_x{N}_x^T\left[N_{\mathrm{xy}}^T\mathrm{\ψ}\right]\mathrm{dl}$

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

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

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

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

$M\stackrel{\·\·}{\mathrm{\ψ}}+K\left(\mathrm{\ψ}\right)\mathrm{\ψ}=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秒时的水波形状，由图可见，本发明方法采用了拉格朗日坐标，可以精确跟踪流体中各个粒子的位移，对于自由水面无需做特殊处理，因此可很方便地描述自由水面的不断变化，这是传统采用欧拉坐标难以实现的。因为采用界带有限元，可以精确地处理不可压缩性，计算精度高，没有体积闭锁，这是传统有限单元难以实现的，且在边界处只需要对刚度矩阵和质量矩阵化行化列处理，实际操作方便。仿真试验表明，本发明方法可以很方便和准确地模拟水中孤立波的传播过程。

Claims (9)

1.一种基于界带有限元和拉格朗日坐标的流体仿真方法，将二维不可压缩流体的计算域Ω按照传统有限元网格剖分成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{\·\·}{\mathrm{\ψ}}+K\left(\mathrm{\ψ}\right)\mathrm{\ψ}=0$

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

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

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

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

K_lψ＝ω²Mψ

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

2.根据权利要求1所述一种基于界带有限元和拉格朗日坐标的流体仿真方法，其特征在于步骤(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相同。

3.根据权利要求1所述一种基于界带有限元和拉格朗日坐标的流体仿真方法，其特征在于步骤(c)中，流体中在坐标(x,y)处质点的位移为：

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

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

4.根据权利要求1所述一种基于界带有限元和拉格朗日坐标的流体仿真方法，其特征在于步骤(d)中，对单元进行积分以计算质量矩阵，积分方案选择数值积分，也即

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

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

5.根据权利要求1所述一种基于界带有限元和拉格朗日坐标的流体仿真方法，其特征在于步骤(e)中，根据步骤(c)的位移场表达式，得到流体自由面上Γ_f的位移u和v，在重力作用下，其刚度矩阵的计算表达式为：

K(ψ)＝K_l+K_n(ψ)

$K_1={\∫}_{{\mathrm{\Γ}}_f}\mathrm{\ρ}(x,y)g{N}_x{N}_x^T\mathrm{dl},K_n\left(\mathrm{\ψ}\right)={\∫}_{{\mathrm{\Γ}}_f}\mathrm{\ρ}(x,y){\mathrm{gN}}_x{N}_x^T\left[N_{\mathrm{xy}}^T\mathrm{\ψ}\right]\mathrm{dl}$

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

6.根据权利要求1所述一种基于界带有限元和拉格朗日坐标的流体仿真方法，其特征在于在步骤(a)中，对流体进行有限元网格剖分，可以采用现有有限元软件，如Ansys软件。

7.根据权利要求1所述一种基于界带有限元和拉格朗日坐标的流体仿真方法，其特征在于在步骤(g)和(i)中，求解微分方程的软件可以采用现有的软件包，如大连理工大学开发的SIPESC软件，SIPESC软件是大连理工大学工程力学系开发的工程计算分析软件平台；其功能包括集成开发环境、面向系统集成的活动流程图定制、工程数据库管理系统、开放式结构有限元分析系统、集成优化计算系统等，其中有限元分析系统中集成了方程求解模块、有限元后处理模块等，可以对本发明中微分方程进行求解；但不仅限于SIPESC软件，只要采用步骤(g)和(i)中的微分方程进行分析，均在本发明权利要求范围之内。

8.根据权利要求1所述一种基于界带有限元和拉格朗日坐标的流体仿真方法，其特征在于步骤(h)中，后处理显示流体的动态运动的软件，可以采用现有的后处理软件，如大连理工大学开发的SIPESC.POST、SIPESC.Jifex软件，SIPESC.POST和SIPESC.Jifex均属于SIPESC软件中有限元系统模块，其功能为有限元结果后处理，可视化动画显示等。

9.根据权利要求1所述一种基于界带有限元和拉格朗日坐标的流体仿真方法，其特征在于在步骤(j)中，求解特征值方程，可以采用现有软件，如大连理工大学开发的SIPESC软件，但不仅限于SIPESC软件，只要采用步骤(j)中的特征值方程进行分析，均在本发明权利要求范围之内。