CN103970989A - 一种基于流固界面一致条件的浸入边界流场计算方法 - Google Patents

Abstract

本发明涉及一种基于流固界面一致条件的浸入边界流场计算方法，属于计算流体力学及其流固耦合模拟技术领域。本发明包括步骤：A、调用网格划分模块，采用两套网格：流场区域，固体边界区域；B、调用流场计算模块，求得流场区域的预测值；C、调用固体边界力密度计算模块，获取作用在固体边界上的作用力密度；D、调用流场速度校正模块，获得流场区域的校正值，并更新流场速度；E、调用结果输出模块，将作用在固体边界上的力以及流场信息输出到文件，供后台读取显示；F、判断是否结束计算。本发明避免使用动网格技术，大量节省计算资源；克服传统浸入边界法近壁区速度求解时复杂的插值运算，不易计算固体边界作用力密度的缺点。

Description

一种基于流固界面一致条件的浸入边界流场计算方法

技术领域

本发明涉及一种基于流固界面一致条件的浸入边界流场计算方法，属于计算流体力学及其流固耦合模拟技术领域。

背景技术

流固耦合问题的难点之一是流体和固体使用不同的数学描述框架。通常，流体运动使用欧拉描述，而固体运动使用拉格朗日描述：传统的浸入边界法，固体边界区域的速度大多采用线性（双线性）插值获得，往往插值过程异常复杂，得到的边界层附近的速度分布误差较大。为此，本发明提出基于流固界面一致条件的浸入边界流场计算方法，避免了近壁区速度求解时复杂的插值过程，计算过程简单，易于程序设计。

浸入边界法提供了一种方式，使这两种框架相互联系在一起。而浸入边界法将浸入流体的固体模化成Navier-Stokes动量方程中的体力，这样整个物理区域（流体区域和固体边界区域）可以看成一个流场进行求解。浸入边界法采用两套网格，整个物理区域在欧拉描述下采用笛卡尔网格离散，固体边界区域在拉格朗日描述下使用适体曲线网格离散。拉格朗日变量和欧拉变量的信息交换通过近似光滑函数实现。两套网格互不关联，流场求解不再使用动网格技术，从而避免了网格畸形以及离散网格几何不守恒等问题，能够有效地处理刚体与流体的相互作用问题。

发明内容

本发明提供了一种基于流固界面一致条件的浸入边界流场计算方法，以用于克服传统浸入边界法在计算刚体与流体耦合作用上的不足。

本发明的技术方案是：一种基于流固界面一致条件的浸入边界流场计算方法，所述方法的具体步骤如下：

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

B、调用初、边值设定模块，对流场区域的边界设定流场边界条件，对流场区域原始变量设定初始值；

C、调用流场速度预测模块，采用分步投影方法，求解不可压缩粘性牛顿流体的流动控制方程，求得流场区域原始变量的预测值；

D、调用固体边界力密度计算模块，通过流固界面速度一致条件，获得固体边界离散的拉格朗日点上的作用力密度；

E、调用流场速度校正模块，计算流场区域的速度校正值，并更新流场速度；

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

G、判断是否结束计算：

如果△t·n＜T，则进入下一时间步，继续执行步骤C、D、E和F；

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

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

所述流场区域包括流体和固体所占据的空间区域。

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

所述步骤B中，所述流场边界条件主要包括（1）与固体接触面的边界速度都设为零，即无滑移无渗透壁面边界条件；（2）已知流速作为计算域进口边界条件；（3）出口采用自由出口边界条件。

所述步骤B中，对流场区域原始变量设定初始值中的原始变量包括速度变量和压力变量；

所述步骤C中，流场速度预测模块通过求解压力泊松方程，获取流场区域的压力，进而获得流场区域的预测速度 ，并通过程序接口，提取笛卡尔网格单元中心上的预测速度，输出到文件fvel.txt。

所述步骤D中，

为满足流固界面无滑移无渗透的一致边界条件，通过流场区域速度和近似光滑函数得到的固体边界拉格朗日点上的速度应该等于给定的固体边界的自然速度，来实现固体边界力密度的求解，其中为固体边界离散曲线网格点的初始构型坐标，为时间，变量下标和分别表示固体边界离散曲线网格第个节点和流场区域欧拉网格第个单元，固体边界力密度可写成矩阵的形式，即求解代数方程组

（1）

式中

（2）

（3）

（4）

（2）式和（3）式中，为流场区域的网格间距，为第段固体边界的面积， 为固体边界第个拉格朗日离散点上的自然速度， 为步骤C中求得的流场区域第个单元中心的预测速度，式（4）中表示固体边界拉格朗日点上的固体边界力密度，假如固体边界一系列的拉格朗日点离散，整个物理域由固定的等距的欧拉网格点离散，（2）式和（3）式中的和为信息转换矩阵，定义如下：

（5）

（6）

式中，函数可表示为：

（7）

其中，通过转换矩阵将拉格朗日点上物理量转换到欧拉网格点上，通过转换矩阵将欧拉网格点上物理量转换到拉格朗日点上；

通过式（1）求得固体边界离散拉格朗日点上的固体边界力密度，并将结果输出到文件sfor.txt。

所述步骤E中，流场区域速度校正值为：

（8）

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

（9）

式中为第个欧拉网格点上下一步的速度，可表示步骤C中的预测速度和式（8）计算得到的校正速度之和。

本发明的有益效果是：

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

2、本发明提出的耦合方法克服传统的浸入边界法求解近壁区速度复杂的插值过程，适应性更强，应用更简便，使之能更有效地预测固体和流体的耦合作用，在流体力学领域和流固耦合领域应用广泛。

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

附图说明

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

图2为本发明中整个物体系统的计算区域的示意图；

图3为本发明中计算得到的不同流向断面的平均速度，图中表示时间平均的方向的速度值。

具体实施方式

实施例1：如图1-3所示，一种基于流固界面一致条件的浸入边界流场计算方法，所述方法的具体步骤如下：

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

B、调用初、边值设定模块，对流场区域的边界设定流场边界条件，对流场区域原始变量设定初始值；

C、调用流场速度预测模块，采用分步投影方法，求解不可压缩粘性牛顿流体的流动控制方程，求得流场区域原始变量的预测值；

D、调用固体边界力密度计算模块，通过流固界面速度一致条件，获得固体边界离散的拉格朗日点上的作用力密度；

E、调用流场速度校正模块，计算流场区域的速度校正值，并更新流场速度；

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

G、判断是否结束计算：

如果△t·n＜T，则进入下一时间步，继续执行步骤C、D、E和F；

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

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

所述流场区域包括流体和固体所占据的空间区域。

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

所述步骤B中，所述流场边界条件主要包括（1）与固体接触面的边界速度都设为零，即无滑移无渗透壁面边界条件；（2）已知流速作为计算域进口边界条件；（3）出口采用自由出口边界条件。

所述步骤B中，对流场区域原始变量设定初始值中的原始变量包括速度变量和压力变量；

所述步骤C中，流场速度预测模块通过求解压力泊松方程，获取流场区域的压力，进而获得流场区域的预测速度，并通过程序接口，提取笛卡尔网格单元中心上的预测速度，输出到文件fvel.txt。

所述步骤D中，

通过流场区域速度和近似光滑函数得到的固体边界拉格朗日点上的速度应该等于给定的固体边界的自然速度，来实现固体边界力密度的求解，其中为固体边界离散曲线网格点的初始构型坐标，为时间，变量下标和分别表示固体边界离散曲线网格第个节点和流场区域欧拉网格第个单元，固体边界力密度可写成矩阵的形式，即求解代数方程组

（1）

式中

（2）

（3）

（4）

（2）式和（3）式中，为流场区域的网格间距，为第段固体边界的面积， 为固体边界第个拉格朗日离散点上的自然速度，为步骤C中求得的流场区域第个单元中心的预测速度，式（4）中表示固体边界拉格朗日点上的固体边界力密度，（2）式和（3）式中的和为信息转换矩阵，定义如下：

（5）

（6）

式中，函数可表示为：

（7）

其中，通过转换矩阵将拉格朗日点上物理量转换到欧拉网格点上，通过转换矩阵将欧拉网格点上物理量转换到拉格朗日点上；

通过式（1）求得固体边界离散拉格朗日点上的固体边界力密度，并将结果输出到文件sfor.txt。

所述步骤E中，流场区域速度校正值为：

（8）

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

（9）

式中为第个欧拉网格点上下一步的速度，可表示步骤C中的预测速度和式（8）计算得到的校正速度之和。

实施例2：如图1-3所示，一种基于流固界面一致条件的浸入边界流场计算方法，所述方法的具体步骤如下：

求解区域内方柱高度为H，计算区域为，流场中浸没一刚性方柱，坐标原点位于方柱左下侧顶点。流动雷诺数定义为，时间步长，方柱边界离散采用间距为0.01H。

S1：网格划分

流场区域（包括流体和刚性方柱区域）在欧拉描述下采用笛卡尔网格划分，网格为均匀四边形网格，网格间距为h=0.01H，并将单元坐标信息输出到文件fcor.txt，固体边界区域在拉格朗日描述下使用适体曲线网格划分，间距为0.01H，相应网格节点信息文件输出到scor.txt文件。

S2：初、边值计算

边值条件：

计算区域左边为进口边界条件，为，，上下为无滑移无穿透边界条件，即，，右边为出口边界条件，。

初值条件：

设定流场区域速度初始值，，流场压力。

S3：流场区域速度预测

采用分步投影方法，数值求解不可压缩粘性牛顿流体的流动控制方程，通过求解压力泊松方程，获取流场区域的压力，进而获得流场的预测速度，并通过程序接口，提取笛卡尔网格单元中心上的预测速度，输出到文件fvel.txt。

S4：固体边界力密度计算

为满足流固界面无滑移无渗透的一致边界条件，通过流场区域速度和近似光滑函数得到的固体边界拉格朗日点上的速度应该等于给定的固体边界的自然速度，来实现固体边界力密度的求解，其中为固体边界离散曲线网格点的初始构型坐标，为时间，变量下标和分别表示固体边界离散曲线网格第个节点和流场区域欧拉网格第个单元。

S4.1：求转换矩阵和

已知固体边界一系列的拉格朗日点离散，整个物理域由固定的等距的欧拉网格点离散，由步骤S1中的单元坐标信息文件fcor.txt和固体边界区域网格节点信息文件scor.txt，由下式（1）和（2）可求得转换矩阵和：

（1）

（2）

式中，函数可表示为：

（3）

S4.2：求系数矩阵A和B

已知流场区域的网格间距h=0.01H，时间步长，方柱边界离散采用间距为0.01H，第段固体边界的面积，固体边界第个拉格朗日离散点上的自然速度，以及步骤S3中求得的流场区域第个单元中心的预测速度，可由下面式（3）和式（4）求得系数矩阵A和B；

（3）

（4）

S4.3：求固体边界力密度

固体边界力密度可写成矩阵的形式，通过高斯-约当方法求解下面的代数方程组（5），即可获得固体边界拉格朗日离散点上的边界力密度，即

（5）

S5：流场区域速度更新

流场区域速度校正值可由下面（6）式获得：

（6）

流场区域速度可由下面（7）式更新

（7）

式中为第个欧拉网格点上下一步的速度，可表示步骤C中的预测速度和式（6）计算得到的校正速度之和。

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

图3显示了不同流向断面时间平均的速度分布。从图可见，在刚性方柱的上方区域，流道束窄，流向速度大大增加，最大值达到进口速度的2.9倍，同时受刚性方柱前缘影响，流动出现分离现象，在和断面都出现负的流向速度，在刚性方柱上表面附近形成明显的漩涡结构。从断面上速度分布看，在刚性方柱后方速度呈现正负交替变化，多处出现速度驻点，预示贴近后壁面局部呈现多个小漩涡结构，体现附壁流动的复杂性和不规则性。在断面时，在左右出现速度驻点，在下方速度为负值，说明流体在刚性方柱后下方形成大的漩涡结构。在和断面时，断面流向速度已趋于平缓，速度分布类似于平板间的泊潇叶流动。

S7：时间推进

在一个时间步内计算完成后，转入下一个时间步，重复上述步骤S3-S6，直到满足计算时间要求停止计算。本例取时间步长=0.005s，计算总时间T=20s，共计时间步。

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

Claims (8)

1.一种基于流固界面一致条件的浸入边界流场计算方法，其特征在于：所述方法的具体步骤如下：

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

B、调用初、边值设定模块，对流场区域的边界设定流场边界条件，对流场区域原始变量设定初始值；

C、调用流场速度预测模块，采用分步投影方法，求解不可压缩粘性牛顿流体的流动控制方程，求得流场区域原始变量的预测值；

D、调用固体边界力密度计算模块，通过流固界面速度一致条件，获得固体边界离散的拉格朗日点上的作用力密度；

E、调用流场速度校正模块，计算流场区域的速度校正值，并更新流场速度；

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

G、判断是否结束计算：

如果△t·n＜T，则进入下一时间步，继续执行步骤C、D、E和F；

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

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

2.根据权利要求1所述的基于流固界面一致条件的浸入边界流场计算方法，其特征在于：所述流场区域包括流体和固体所占据的空间区域。

3.根据权利要求1所述的基于流固界面一致条件的浸入边界流场计算方法，其特征在于：所述步骤A中，流场区域在欧拉描述下采用笛卡尔网格划分，其网格单元中心或节点上的流场变量称为欧拉变量，并将单元坐标信息输出到文件fcor.txt；固体边界区域在拉格朗日描述下使用适体曲线网格划分，其网格节点上的变量称为拉格朗日变量，相应网格节点坐标信息输出到文件scor.txt。

4.根据权利要求1所述的基于流固界面一致条件的浸入边界流场计算方法，其特征在于：所述步骤B中，所述流场边界条件主要包括（1）与固体接触面的边界速度都设为零，即无滑移无渗透壁面边界条件；（2）已知流速作为计算域进口边界条件；（3）出口采用自由出口边界条件。

5.根据权利要求1所述的基于流固界面一致条件的浸入边界流场计算方法，其特征在于：所述步骤B中，对流场区域原始变量设定初始值中的原始变量包括速度变量和压力变量。

6.根据权利要求1所述的基于流固界面一致条件的浸入边界流场计算方法，其特征在于：所述步骤C中，流场速度预测模块通过求解压力泊松方程，获取流场区域的压力，进而获得流场区域的预测速度 ，并通过程序接口，提取笛卡尔网格单元中心上的预测速度，输出到文件fvel.txt。

7.根据权利要求1所述的基于流固界面一致条件的浸入边界流场计算方法，其特征在于：所述步骤D中，

通过流场区域速度和近似光滑函数得到的固体边界拉格朗日点上的速度应该等于给定的固体边界的自然速度，来实现固体边界力密度的求解，其中为固体边界离散曲线网格点的初始构型坐标，为时间，变量下标和分别表示固体边界离散曲线网格第个节点和流场区域欧拉网格第个单元，固体边界力密度可写成矩阵的形式，即求解代数方程组

（1）

式中

（2）

（3）

（4）

（2）式和（3）式中，为流场区域的网格间距，为第段固体边界的面积， 为固体边界第个拉格朗日离散点上的自然速度，为步骤C中求得的流场区域第个单元中心的预测速度，式（4）中表示固体边界拉格朗日点上的固体边界力密度，（2）式和（3）式中的和为信息转换矩阵，定义如下：

（5）

（6）

式中，函数可表示为：

（7）

其中，通过转换矩阵将拉格朗日点上物理量转换到欧拉网格点上，通过转换矩阵将欧拉网格点上物理量转换到拉格朗日点上；

通过式（1）求得固体边界离散拉格朗日点上的固体边界力密度，并将结果输出到文件sfor.txt。

8.根据权利要求1所述的一种基于流固界面一致条件的浸入边界流场计算方法，其特征在于：所述步骤E中，流场区域速度校正值为：

（8）

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

（9）

式中为第个欧拉网格点上下一步的速度，可表示步骤C中的预测速度和式（8）计算得到的校正速度之和。