CN104850689A - 一种基于固定网格技术的流固耦合计算方法 - Google Patents

Abstract

本发明涉及一种基于固定网格技术的流固耦合计算方法，属于流固耦合模拟技术领域。本发明包括步骤：A、调用网格划分模块，采用两套网格：流体和固体组成的整个物理区域，固体区域；B、调用流场计算CFD模块，更新流场状态变量，同时获得流体作用于固体边界节点上的力；C、调用固体计算CSD模块，更新固体状态变量，同时获取固体边界节点上的位移和速度；D、调用结果输出模块，将固体状态变量和流场信息输出到文件，供后处理使用；E、判断是否结束计算。本发明成功避免使用动网格模块，大大节省了网格更新所需要的计算资源；克服了传统的基于移动网格技术时复杂的内部网格更新与边界网格控制运算，不容易准确计算流固界面作用力的缺点。

Description

一种基于固定网格技术的流固耦合计算方法

技术领域

本发明涉及一种基于固定网格技术的流固耦合计算方法，属于流固耦合模拟技术领域。

背景技术

流固耦合问题的难点之一是流体和固体使用不同的数学描述框架。通常，流体运动使用欧拉描述，而固体运动使用拉格朗日描述：传统的流固耦合方法，流体区域的网格随固体边界的运动而动态更新，网格运算和网格质量控制过程异常复杂，得到的边界层附近的速度分布误差较大。为此，本发明提出基于固定网格技术的流固耦合计算方法，避免了使用动网格计算模块，计算过程简单，易于程序设计。

发明内容

本发明提供了一种基于固定网格技术的流固耦合计算方法，以用于克服传统基于移动网格技术在计算固体与流体耦合作用上的不足。

本发明的技术方案是：一种基于固定网格技术的流固耦合计算方法，所述方法的具体步骤如下：

Step1、调用网格划分模块，采用两套网格：流场区域，固体区域；其中流场区域在欧拉描述下采用笛卡尔网格离散，固体区域在拉格朗日描述下使用适体有限元网格离散；

Step2、调用流场CFD(计算流体动力学)计算模块，采用基于分步投影的浸入边界方法，求解不可压缩粘性牛顿流体的流动控制方程，更新流场状态变量，同时获得流体作用于固体边界节点上的力；

Step3、调用固体CSD(计算固体动力学)计算模块，采用有限单元方法，求解固体动力学控制方程，更新固体状态变量，同时获取固体边界节点上的位移和速度；

Step4、将固体状态变量和流场信息输出到文件，供后处理使用；

Step5、判断是否结束计算：

如果△t·n＜T，则进入下一时间步，继续执行步骤Step2、Step3和Step4；

如果△t·n≥T，则结束整个计算；

其中△t为时间步长，T为要求计算的总物理时间，n为时间步数。

所述流场区域包括流体和固体所占据的空间区域，固体区域可包括刚体区域和弹性体区域。

所述步骤Step1中，流场区域在欧拉描述下采用笛卡尔网格划分，其网格单元中心或节点上的流场变量称为欧拉变量，并将网格节点坐标信息x_j输出到文件fcor.txt；固体区域在拉格朗日描述下使用适体有限元网格划分，其网格节点上的变量称为拉格朗日变量，相应网格节点坐标信息X输出到文件scor.txt，其中固体边界节点坐标信息可用表示。

所述步骤Step1中，流场区域在欧拉描述下采用笛卡尔网格划分，其网格单元中心或节点上的流场变量称为欧拉变量，并将网格节点坐标信息x_j输出到文件fcor.txt；固体区域在拉格朗日描述下使用适体有限元网格划分，其网格节点上的变量称为拉格朗日变量，相应网格节点坐标信息X输出到文件scor.txt，其中固体边界节点坐标信息可用表示。

所述步骤Step2中，通过基于分步投影浸入边界方法的CFD计算，实现两大功能：

一是获得流体作用于固体边界节点上的力：

通过流场区域预测速度u′(x_j,t)和近似光滑函数得到的固体边界拉格朗日节点上的速度应该等于给定的固体边界节点的自然速度来实现固体边界节点上力密度的求解，并将结果输出到文件sfor.txt，供步骤Step3中CSD计算模块使用；其中，为固体边界节点坐标的解释，t为时间；

二是更新流场状态变量：

流场区域速度校正值Δu(x_j,t)为：

$\mathrm{\Δu}(x_j,t)=\underset{i}{\mathrm{\Σ}}{C}_{\mathrm{ji}}F(X_B^i,t)\mathrm{h\Δ}{s}_i\mathrm{\Δt}(i=\mathrm{1,2},...M;j=\mathrm{1,2},...,N)---\left(1\right)$

式中，Δs_i为第i段固体边界的面积，h为流体网格间距，M为固体边界节点数，N为流体欧拉网格节点数，C_ji为信息转换矩阵，定义如下：

$C_{\mathrm{ji}}=C(X_B^i-x_j)=\frac{1}{h^3}\mathrm{\φ}\left(\frac{X_B^i-x_j}{h}\right)\mathrm{\φ}\left(\frac{Y_B^i-y_j}{h}\right)\mathrm{\φ}\left(\frac{Z_B^i-z_j}{h}\right)---\left(2\right)$

式中，分别为固体边界节点坐标在x,y,z方向的分量，x_j,y_j,z_j分别为欧拉网格节点坐标x_j在x,y,z方向的分量；函数φ可表示为：

$\mathrm{\&phi;}\left(r\right)=\left\{\begin{array}{cc}1/8(3-2|r|+\sqrt{1+4\left|r\right|-4r^2}),& 0\&le;\left|r\right|<1\\ 1/8(5-2|r|+\sqrt{-7+12\left|r\right|-4r^2}),& 1\&le;\left|r\right|<2\\ 0,& 2\&le;\left|r\right|\end{array}\right.---\left(3\right)$

式中，r为函数φ的自变量；

流场区域速度可由下式更新，

u(x_j,t)＝u′(x_j,t)+Δu(x_j,t)   (4)

并将结果输出到文件fvel.txt，供步骤Step4中的后处理使用。

所述步骤Step3中，通过有限单元方法，求解离散后的固体动力学控制方程，实现两大功能，更新固体状态变量，同时获取固体边界节点上的位移和速度；

一是获取固体边界节点上的位移和速度：

采用有限单元数值离散方法，可获得固体的一组离散非线性方程组，写成矩阵形式有：

$M^s{\stackrel{\&CenterDot;\&CenterDot;}{X}}_t^s+C^s{\stackrel{\&CenterDot;}{X}}_t^s+K^s{X}_t^s=Q_t^s---\left(5\right)$

式中，M^s为质量矩阵，C^s阻尼矩阵，K^s为刚度矩阵，Q^s为外力矢量，包括流体作用于固体边界的力密度和其他外力矢量，和分别代表固体t时刻的加速度矢量，速度矢量和位移矢量；

式(5)求解在时间推进上采用纽马克法，可获得t+Δt时刻有限元节点新的位移矢量速度矢量和加速度矢量提取t+Δt时刻边界节点上的位移矢量速度矢量令并将结果输出到文件sdisvel.txt文件，供步骤Step2中CFD计算模块使用；

二是更新固体状态变量：

通过对式(5)的求解，获得下一时刻有限元节点新的位移矢量通过固体本构方程，可进一步获得固体各节点的应力应变等固体其他状态变量的值，并将结果输出到文件sstress.txt，供步骤Step4中的后处理使用。

本发明的有益效果是：

1、避免使用动网格技术，大量节省计算资源：传统的基于移动网格技术的流固耦合方法，需要借助动网格技术，而对于复杂的流固耦合问题，大幅度的固体运动往往导致流体区域网格更新的失败，而本发明正是弥补这一重要缺陷，在固体与流体耦合作用过程中成功避免使用动网格技术。

2、即使近年来也有少量成功应用动网格技术来执行流固耦合计算，但也很难保证耦合界面信息交换时能量守恒和边界层网格质量有效控制等关键问题，大大降低了模拟精度。本发明采用适宜的近似光滑函数，利用流固界面一致条件，确保流体与固体界面作用力平衡，保证耦合计算的有效性。

3、本发明使用的计算流程物理概念清晰，计算过程简单，易于程序设计，便于处理复杂的流固耦合问题。

附图说明

图1为本发明中的流程图；

图2为本发明中算例几何尺寸和整个系统计算区域的示意图；

图3为本发明中计算得到的悬臂梁自由端的位移随时间的演化过程。

具体实施方式

实施例1：如图1-3所示，一种基于固定网格技术的流固耦合计算方法，所述方法的具体步骤如下：

Step1、调用网格划分模块，采用两套网格：流场区域，固体区域；其中流场区域在欧拉描述下采用笛卡尔网格离散，固体区域在拉格朗日描述下使用适体有限元网格离散；

Step2、调用流场CFD(计算流体动力学)计算模块，采用基于分步投影的浸入边界方法，求解不可压缩粘性牛顿流体的流动控制方程，更新流场状态变量，同时获得流体作用于固体边界节点上的力；

Step3、调用固体CSD(计算固体动力学)计算模块，采用有限单元方法，求解固体动力学控制方程，更新固体状态变量，同时获取固体边界节点上的位移和速度；

Step4、将固体状态变量和流场信息输出到文件，供后处理使用；

Step5、判断是否结束计算：

如果△t·n＜T，则进入下一时间步，继续执行步骤Step2、Step3和Step4；

如果△t·n≥T，则结束整个计算；

其中△t为时间步长，T为要求计算的总物理时间，n为时间步数。

所述流场区域包括流体和固体所占据的空间区域，固体区域可包括刚体区域和弹性体区域。

所述步骤Step1中，流场区域在欧拉描述下采用笛卡尔网格划分，其网格单元中心或节点上的流场变量称为欧拉变量，并将网格节点坐标信息x_j输出到文件fcor.txt；固体区域在拉格朗日描述下使用适体有限元网格划分，其网格节点上的变量称为拉格朗日变量，相应网格节点坐标信息X输出到文件scor.txt，其中固体边界节点坐标信息可用表示。

所述步骤Step1中，流场区域在欧拉描述下采用笛卡尔网格划分，其网格单元中心或节点上的流场变量称为欧拉变量，并将网格节点坐标信息x_j输出到文件fcor.txt；固体区域在拉格朗日描述下使用适体有限元网格划分，其网格节点上的变量称为拉格朗日变量，相应网格节点坐标信息X输出到文件scor.txt，其中固体边界节点坐标信息可用表示。

所述步骤Step2中，通过基于分步投影浸入边界方法的CFD计算，实现两大功能：

一是获得流体作用于固体边界节点上的力：

通过流场区域预测速度u′(x_j,t)和近似光滑函数得到的固体边界拉格朗日节点上的速度应该等于给定的固体边界节点的自然速度来实现固体边界节点上力密度的求解，并将结果输出到文件sfor.txt，供步骤Step3中CSD计算模块使用；其中，为固体边界节点坐标的解释，t为时间；

二是更新流场状态变量：

流场区域速度校正值Δu(x_j,t)为：

$\mathrm{\&Delta;u}(x_j,t)=\underset{i}{\mathrm{\&Sigma;}}{C}_{\mathrm{ji}}F(X_B^i,t)\mathrm{h\&Delta;}{s}_i\mathrm{\&Delta;t}(i=\mathrm{1,2},...M;j=\mathrm{1,2},...,N)---\left(1\right)$

式中，Δs_i为第i段固体边界的面积，h为流体网格间距，M为固体边界节点数，N为流体欧拉网格节点数，C_ji为信息转换矩阵，定义如下：

$C_{\mathrm{ji}}=C(X_B^i-x_j)=\frac{1}{h^3}\mathrm{\&phi;}\left(\frac{X_B^i-x_j}{h}\right)\mathrm{\&phi;}\left(\frac{Y_B^i-y_j}{h}\right)\mathrm{\&phi;}\left(\frac{Z_B^i-z_j}{h}\right)---\left(2\right)$

式中，分别为固体边界节点坐标在x,y,z方向的分量，x_j,y_j,z_j分别为欧拉网格节点坐标x_j在x,y,z方向的分量；函数φ可表示为：

$\mathrm{\&phi;}\left(r\right)=\left\{\begin{array}{cc}1/8(3-2|r|+\sqrt{1+4\left|r\right|-4r^2}),& 0\&le;\left|r\right|<1\\ 1/8(5-2|r|+\sqrt{-7+12\left|r\right|-4r^2}),& 1\&le;\left|r\right|<2\\ 0,& 2\&le;\left|r\right|\end{array}\right.---\left(3\right)$

式中，r为函数φ的自变量；

流场区域速度可由下式更新，

u(x_j,t)＝u′(x_j,t)+Δu(x_j,t)   (4)

并将结果输出到文件fvel.txt，供步骤Step4中的后处理使用。

所述步骤Step3中，通过有限单元方法，求解离散后的固体动力学控制方程，实现两大功能，更新固体状态变量，同时获取固体边界节点上的位移和速度；

一是获取固体边界节点上的位移和速度：

采用有限单元数值离散方法，可获得固体的一组离散非线性方程组，写成矩阵形式有：

$M^s{\stackrel{\&CenterDot;\&CenterDot;}{X}}_t^s+C^s{\stackrel{\&CenterDot;}{X}}_t^s+K^s{X}_t^s=Q_t^s---\left(5\right)$

式中，M^s为质量矩阵，C^s阻尼矩阵，K^s为刚度矩阵，Q^s为外力矢量，包括流体作用于固体边界的力密度和其他外力矢量，和分别代表固体t时刻的加速度矢量，速度矢量和位移矢量；

式(5)求解在时间推进上采用纽马克法，可获得t+Δt时刻有限元节点新的位移矢量速度矢量和加速度矢量提取t+Δt时刻边界节点上的位移矢量速度矢量令并将结果输出到文件sdisvel.txt文件，供步骤Step2中CFD计算模块使用；

二是更新固体状态变量：

通过对式(5)的求解，获得下一时刻有限元节点新的位移矢量通过固体本构方程，可进一步获得固体各节点的应力应变等固体其他状态变量的值，并将结果输出到文件sstress.txt，供步骤Step4中的后处理使用。

实施例2：如图1-3所示，一种基于固定网格技术的流固耦合计算方法，所述方法的具体步骤如下：

一静止刚性方柱浸没在不可压缩流场中，左边进口流动假设为均匀流场，速度大小为v_x＝v₀,v_y＝0。方柱正后方中心附着一个弹性悬臂梁(固体密度ρ_s＝100Kg/m³，弹性模量E＝25Gpa,泊松比ν＝0.35)。在静止时梁轴线平行于远场流动方向，进口速度v₀＝51.3cm/s，雷诺数Re＝ρ^fDv₀/μ^f＝333，式中D＝1cm为刚性方柱的边长，柔性悬臂梁高度为0.06D。

S1：网格划分

流场区域(包括流体和固体区域)在欧拉描述下采用笛卡尔网格划分，网格为均匀四边形网格，网格间距为h＝0.01D，并将单元坐标信息输出到文件fcor.txt，固体区域(此算例包括刚体区域和弹性体区域，即刚性方柱和弹性悬臂梁)，在拉格朗日描述下使用有限元四边形网格划分，间距为0.01D，相应网格节点信息文件输出到scor.txt文件。

S2：调用流场CFD计算模块，更新流场状态变量，同时获得流体作用于固体边界节点上的力；

通过流场区域预测速度u′(x_j,t)和近似光滑函数得到的固体边界拉格朗日节点上的速度应该等于给定的固体边界节点的自然速度。对刚性方柱，其固体边界节点的自然速度对弹性悬臂梁边界节点的自然速度即等于固体CSD计算得到的固体边界的速度，在初始时刻，可取从而实现固体边界节点上力密度的求解，并将结果输出到文件sfor.txt，供下一步骤S3中CSD计算模块使用。同时，更新流场速度，并将结果输出到文件fvel.txt，供步骤S4中的后处理使用。

S3：调用固体CSD计算模块，更新固体状态变量，同时获取固体边界节点上的位移和速度

采用有限单元数值离散方法，可获得固体的一组离散非线性方程组，写成矩阵形式有：式中，M^s为质量矩阵，C^s阻尼矩阵，K^s为刚度矩阵，Q^s为外力矢量，包括流体作用于固体边界的力密度和其他外力矢量，和分别代表固体t时刻的加速度矢量，速度矢量和位移矢量。

在时间推进上采用纽马克法，可获得t+Δt时刻有限元节点新的位移矢量速度矢量和加速度矢量提取t+Δt时刻边界节点上的位移矢量速度矢量令并将结果输出到文件sdisvel.txt文件，供步骤S2中CFD计算模块使用。

S4：调用结果输出模块，将作用在固体边界上的力以及流场信息输出到文件，供后处理软件读取显示。

S5：时间推进

在一个时间步内计算完成后，转入下一个时间步，重复上述步骤S2-S4，直到满足计算时间要求停止计算。本例取时间步长Δt＝0.001s，计算总时间T＝12s，共计时间步n＝12000。

图3为悬臂梁自由端的位移随时间的演化过程。从图可见，弹性悬臂梁在流体动力作用下经历了一个长期的不稳定振动过程，约8s后振动趋于稳定的周期振动，梁自由端的振幅约1.4cm左右。

上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

Claims (5)

1.一种基于固定网格技术的流固耦合计算方法，其特征在于：所述方法的具体步骤如下：

Step1、调用网格划分模块，采用两套网格：流场区域，固体区域；其中流场区域在欧拉描述下采用笛卡尔网格离散，固体区域在拉格朗日描述下使用适体有限元网格离散；

Step2、调用流场CFD计算模块，采用基于分步投影的浸入边界方法，求解不可压缩粘性牛顿流体的流动控制方程，更新流场状态变量，同时获得流体作用于固体边界节点上的力；

Step3、调用固体CSD计算模块，采用有限单元方法，求解固体动力学控制方程，更新固体状态变量，同时获取固体边界节点上的位移和速度；

Step4、将固体状态变量和流场信息输出到文件，供后处理使用；

Step5、判断是否结束计算：

如果△t·n＜T，则进入下一时间步，继续执行步骤Step2、Step3和Step4；

如果△t·n≥T，则结束整个计算；

其中△t为时间步长，T为要求计算的总物理时间，n为时间步数。

2.根据权利要求1所述的基于固定网格技术的流固耦合计算方法，其特征在于：所述流场区域包括流体和固体所占据的空间区域，固体区域可包括刚体区域和弹性体区域。

3.根据权利要求1所述的基于固定网格技术的流固耦合计算方法，其特征在于：所述步骤Step1中，流场区域在欧拉描述下采用笛卡尔网格划分，其网格单元中心或节点上的流场变量称为欧拉变量，并将网格节点坐标信息x_j输出到文件fcor.txt；固体区域在拉格朗日描述下使用适体有限元网格划分，其网格节点上的变量称为拉格朗日变量，相应网格节点坐标信息X输出到文件scor.txt，其中固体边界节点坐标信息可用表示。

4.根据权利要求1所述的基于固定网格技术的流固耦合计算方法，其特征在于：所述步骤Step2中，通过基于分步投影浸入边界方法的CFD计算，实现两大功能：

一是获得流体作用于固体边界节点上的力：

通过流场区域预测速度u′(x_j,t)和近似光滑函数得到的固体边界拉格朗日节点上的速度应该等于给定的固体边界节点的自然速度来实现固体边界节点上力密度的求解，并将结果输出到文件sfor.txt，供步骤Step3中CSD计算模块使用；其中，为固体边界节点坐标的解释，t为时间；

二是更新流场状态变量：

流场区域速度校正值Δu(x_j,t)为：

$\mathrm{\&Delta;u}(x_j,t)=\underset{i}{\mathrm{\&Sigma;}}{C}_{\mathrm{ji}}F(X_B^i,t)\mathrm{h\&Delta;}{s}_i\mathrm{\&Delta;t}(i=\mathrm{1,2},...M;j=\mathrm{1,2},...,N)---\left(1\right)$

式中，Δs_i为第i段固体边界的面积，h为流体网格间距，M为固体边界节点数，N为流体欧拉网格节点数，C_ji为信息转换矩阵，定义如下：

$C_{\mathrm{ji}}=C(X_B^i-x_j)=\frac{1}{h^3}\mathrm{\&phi;}\left(\frac{X_B^i-x_j}{h}\right)\mathrm{\&phi;}\left(\frac{Y_B^i-y_i}{h}\right)\mathrm{\&phi;}\left(\frac{Z_B^i-z_j}{h}\right)---\left(2\right)$

式中，分别为固体边界节点坐标在x,y,z方向的分量，x_j,y_j,z_j分别为欧拉网格节点坐标x_j在x,y,z方向的分量；函数φ可表示为：

$\mathrm{\&phi;}\left(r\right)=\left\{\begin{array}{cc}1/8(3-2|r|+\sqrt{1+4\left|r\right|-4r^2}),& 0\&le;\left|r\right|<1\\ 1/8(5-2|r|+\sqrt{-7+12\left|r\right|-4r^2}),& 1\&le;\left|r\right|<2\\ 0,& 2\&le;\left|r\right|\end{array}\right.---\left(3\right)$

式中，r为函数φ的自变量；

流场区域速度可由下式更新，

u(x_j,t)＝u′(x_j,t)+Δu(x_j,t)   (4)

并将结果输出到文件fvel.txt，供步骤Step4中的后处理使用。

5.根据权利要求1所述的基于固定网格技术的流固耦合计算方法，其特征在于：所述步骤Step3中，通过有限单元方法，求解离散后的固体动力学控制方程，实现两大功能，更新固体状态变量，同时获取固体边界节点上的位移和速度；

一是获取固体边界节点上的位移和速度：

采用有限单元数值离散方法，可获得固体的一组离散非线性方程组，写成矩阵形式有：

$M^s{\stackrel{\&CenterDot;\&CenterDot;}{X}}_t^s+C^s{\stackrel{\&CenterDot;}{X}}_t^s+K^s{X}_t^s=Q_t^s---\left(5\right)$

式中，M^s为质量矩阵，C^s阻尼矩阵，K^s为刚度矩阵，Q^s为外力矢量，包括流体作用于固体边界的力密度和其他外力矢量，和分别代表固体t时刻的加速度矢量，速度矢量和位移矢量；

式(5)求解在时间推进上采用纽马克法，可获得t+Δt时刻有限元节点新的位移矢量速度矢量和加速度矢量提取t+Δt时刻边界节点上的位移矢量速度矢量令并将结果输出到文件sdisvel.txt文件，供步骤Step2中CFD计算模块使用；

二是更新固体状态变量：

通过对式(5)的求解，获得下一时刻有限元节点新的位移矢量通过固体本构方程，可进一步获得固体各节点的应力应变等固体其他状态变量的值，并将结果输出到文件sstress.txt，供步骤Step4中的后处理使用。