一种新的投影浸入边界法的速度校正方法

摘要

本发明涉及一种新的投影浸入边界法的速度校正方法，属于计算流体力学及其流固耦合模拟技术领域。本发明包括A、调用网格划分模块，采用两套网格：流场区域欧拉网格，固体边界区域拉格朗日网格；B、调用流场速度预测模块，求得流场区域的预测值；C、调用欧拉网格点属性判断模块，判断欧拉网格点属于固体内部区域还是固体外部区域；D、调用流场速度校正模块，更新流场速度；E、调用结果输出模块，将作用在流固边界上的力以及流场信息输出到文件，供后台读取显示；F、判断是否结束计算。本发明克服传统投影浸入边界法近壁区速度求解时复杂的插值运算，节约了计算资源，提高了计算效率。

说明书

技术领域

本发明涉及一种新的投影浸入边界法的速度校正方法，属于计算流体力学及其流固耦合 模拟技术领域。

背景技术

流固耦合问题的难点之一是流体和固体使用不同的数学描述框架。通常，流体运动使用 欧拉描述，而固体运动使用拉格朗日描述：传统的浸入边界法，固体边界区域流场速度的校 正大多采用线性(双线性)插值获得，往往插值过程计算异常复杂，计算工作量大，计算 效率低。为此，本发明提出一种新的投影浸入边界法的速度校正方法，将流场区域划分为纯 流场区域和次流场区域，纯流场区域速度不再校正，对次流场区域，先进行欧拉网格点属性 判断，如果属于固体内部，直接让速度等于零，如果属于固体外部，再根据流固边界作用力 密度进行修正。避免了近壁区速度求解时复杂的插值过程，计算工作量大大降低，提高计算 效率。

发明内容

本发明提供了一种新的投影浸入边界法的速度校正方法，以用于克服传统浸入边界法在 近壁区速度求解方法上计算异常复杂，计算工作量大，计算效率低的问题，避免了近壁区 速度求解时复杂的插值过程，计算工作量大大降低，提高计算效率。

本发明新的投影浸入边界法的速度校正方法是这样实现的：所述方法的具体步骤如下：

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

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

C、调用欧拉网格点属性判断模块，判断欧拉网格点属于固体内部区域还是固体外部区 域；

D、调用流场速度校正模块，更新流场速度；

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

F、判断是否结束计算：

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

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

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

所述步骤A中，所述流场区域包括流体和固体所占据的空间区域；在空间上将流场区域 划分为纯流场区域和次流场区域，纯流场区域指不受固壁边界层影响并且不包含固体在内的 空间区域，次流场区域指受固壁边界层影响并完全包含固体在内的空间区域。

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

所述步骤B中，流场速度预测模块是指在初、边值条件下采用分步投影方法求解不可压 缩粘性牛顿流体的流动控制方程，获取流场区域的压力，进而获得流场区域的预测速度 u′(x_j,t)，并通过程序接口，提取笛卡尔网格节点上的预测速度u′(x_j,t)，输出到文件fvel.txt。

所述步骤C中，欧拉网格点属性判断模块具体步骤如下：

C1、固体边界离散节点外法线方向求解模块；对任一固体边界离散节点X(s_i)，其相邻两 个节点分别为X(s_i-1)和X(s_i+1)，通过此三点作一条抛物线近似该段曲线，以s为参数，该抛 物线方程表示为

$\left(\begin{array}{c}X\left(s\right)=a_x^i{s}^2+b_x^{i}s+c_x^i\\ Y\left(s\right)=a_y^i{s}^2+b_y^{i}s+c_y^i\end{array}\right)---\left(1\right)$

式中为二次抛物线方程(1)的系数，由此三点的坐标唯一确定；

该段曲线上任意一点的外法线方向表示为

$\left(\begin{array}{c}{N}_x=\frac{-Y_s}{\sqrt{X_s^2+Y_s^2}}\\ N_y=\frac{X_s}{\sqrt{X_s^2+Y_s^2}}\end{array}\right)---\left(2\right)$

此处，分别表示x,y方向的单位矢量，式中X_s,Y_s分别表示(1)式中X(s)和Y(s)对 参数s求一次导数；

C2、任一流体网格节点方向导数求解模块；对任一流体网格节点x_j，求得与之最近的固 体边界离散节点X(s_min)，从该固体边界离散节点X(s_min)到流体网格节点x_j的方向导数 ${\overrightarrow{n}}^j=(n_x^j\overrightarrow{i},n_y^j\overrightarrow{j})$表示为

$\left(\begin{array}{c}{n}_x^j=\frac{x_j-X\left(s_{\min }\right)}{\sqrt{{(x_j-X\left(s_{\min }\right))}^2+{(y_j-Y\left(s_{\min }\right))}^2}}\\ n_y^j=\frac{y_j-Y\left(s_i^{\min }\right)}{\sqrt{{(x_j-X\left(s_{\min }\right))}^2+{(y_j-Y\left(s_{\min }\right))}^2}}\end{array}\right)---\left(3\right)$

依据式(2)，对固体边界离散节点X(s_min)，其外法线方向表示为

$\left(\begin{array}{c}{N}_x^{\min }=\frac{-Y_s^{\min }}{\sqrt{{\left(X_s^{\min }\right)}^2+{\left(Y_s^{\min }\right)}^2}}\\ N_y^{\min }=\frac{X_s^i}{\sqrt{\sqrt{{\left(X_s^{\min }\right)}^2+{\left(Y_s^{\min }\right)}^2}}}\end{array}\right)---\left(4\right)$

式中分别表示(1)式中X(s)和Y(s)对参数s的一次导数在离散节点X(s_min)的 值；

C3、任一流体网格节点属性判断模块；对任一流体网格节点x_j，定义

$\delta =n_x^j·N_x^{\min }+n_y^j·N_y^{\min }---\left(5\right)$

如果δ≥0，则认为该流体网格节点位于固体外部，如果δ<0，则认为该流体网格节点 位于固体内部。

所述步骤D中，更新流场速度的具体步骤如下：

D1、流固界面力密度求解模块；通过流场次区域速度u′(x_j,t)和δ(x-X(s,t))近似光 滑函数得到的固体边界拉格朗日点上的速度u(X(s_i)等于给定的固体边界的自然速度 来实现流固边界力密度的求解，并将结果输出到文件sfor.txt；其中s为固体边 界离散曲线网格点的初始构型坐标，t为时间，变量下标i和j分别表示固体边界离散曲线网 格第i个节点和流场区域欧拉网格第j个单元；

D2、流场区域速度校正模块；

D2.1、如果流体网格节点x_j位于纯流场区域，该节点速度不需要校正，即

u(x_j,t)＝u′(x_j,t)             (6)

D2.2、如果流体网格节点x_j位于次流场区域，调用欧拉网格点属性判断模块：

D2.2.1、如果流体网格节点x_j位于固体内部，流体网格节点x_j速度式中u(x_j,t)由 下式更新

u(x_j,t)＝0             (7)

D2.2.2、如果流体网格节点x_j位于固体外部，其速度校正值Δu(x_j,t)为：

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

式中Δs_i为第i段固体边界的面积，式中的C_ji为信息转换矩阵，定义如下：

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

式中，h为流场区域的网格间距，函数φ表示为：

$\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|<0\\ 0,& 2\le \left|r\right|\end{array}\right)---\left(10\right)$

流体网格节点x_j速度u(x_j,t)由下式更新

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

本发明的有益效果是：

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

2、本发明提出的投影浸入边界法的速度校正方法克服传统的浸入边界法求解近壁区速度 复杂的插值过程，采用适宜的近似光滑函数，适应性更强，利用流固界面一致条件，确保流 固耦合系统能量守恒，保证耦合计算的有效性。

3、将流场区域划分为纯流场区域和次流场区域，速度校正针对性更强，大大提高计算效 率，使之能更有效地预测固体和流体的耦合作用，在流体力学领域和流固耦合领域将得到广 泛应用。

附图说明

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

图2为本发明中整个流固耦合系统的计算区域的示意图；

图3、图4为本发明中计算得到的某一时刻流固边界的作用力密度；

图5为本发明中计算得到的某一时刻流场的速度分布。

具体实施方式

实施例1：如图1-5所示，一种新的投影浸入边界法的速度校正方法，所述方法的具体步 骤如下：

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

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

C、调用欧拉网格点属性判断模块，判断欧拉网格点属于固体内部区域还是固体外部区 域；

D、调用流场速度校正模块，更新流场速度；

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

F、判断是否结束计算：

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

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

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

所述步骤A中，所述流场区域包括流体和固体所占据的空间区域；在空间上将流场区域 划分为纯流场区域和次流场区域，纯流场区域指不受固壁边界层影响并且不包含固体在内的 空间区域，次流场区域指受固壁边界层影响并完全包含固体在内的空间区域。

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

所述步骤B中，流场速度预测模块是指在初、边值条件下采用分步投影方法求解不可压 缩粘性牛顿流体的流动控制方程，获取流场区域的压力，进而获得流场区域的预测速度 u′(x_j,t)，并通过程序接口，提取笛卡尔网格节点上的预测速度u′(x_j,t)，输出到文件fvel.txt。

所述步骤C中，欧拉网格点属性判断模块具体步骤如下：

C1、固体边界离散节点外法线方向求解模块；对任一固体边界离散节点X(s_i)，其相邻两 个节点分别为X(s_i-1)和X(s_i+1)，通过此三点作一条抛物线近似该段曲线，以s为参数，该抛 物线方程表示为

$\left(\begin{array}{c}X\left(s\right)=a_x^i{s}^2+b_x^{i}s+c_x^i\\ Y\left(s\right)=a_y^i{s}^2+b_y^{i}s+c_y^i\end{array}\right)---\left(1\right)$

式中为二次抛物线方程(1)的系数，由此三点的坐标唯一确定；

该段曲线上任意一点的外法线方向表示为

$\left(\begin{array}{c}{N}_x=\frac{-Y_s}{\sqrt{X_s^2+Y_s^2}}\\ N_y=\frac{X_s}{\sqrt{X_s^2+Y_s^2}}\end{array}\right)---\left(2\right)$

此处，分别表示x,y方向的单位矢量，式中X_s,Y_s分别表示(1)式中X(s)和Y(s)对 参数s求一次导数；

C2、任一流体网格节点方向导数求解模块；对任一流体网格节点x_j，求得与之最近的固 体边界离散节点X(s_min)，从该固体边界离散节点X(s_min)到流体网格节点x_j的方向导数 ${\overrightarrow{n}}^j=(n_x^j\overrightarrow{i},n_y^j\overrightarrow{j})$表示为

$\left(\begin{array}{c}{n}_x^j=\frac{x_j-X\left(s_{\min }\right)}{\sqrt{{(x_j-X\left(s_{\min }\right))}^2+{(y_j-Y\left(s_{\min }\right))}^2}}\\ n_y^j=\frac{y_j-Y\left(s_i^{\min }\right)}{\sqrt{{(x_j-X\left(s_{\min }\right))}^2+{(y_j-Y\left(s_{\min }\right))}^2}}\end{array}\right)---\left(3\right)$

依据式(2)，对固体边界离散节点X(s_min)，其外法线方向表示为

$\left(\begin{array}{c}{N}_x^{\min }=\frac{-Y_s^{\min }}{\sqrt{{\left(X_s^{\min }\right)}^2+{\left(Y_s^{\min }\right)}^2}}\\ N_y^{\min }=\frac{X_s^i}{\sqrt{\sqrt{{\left(X_s^{\min }\right)}^2+{\left(Y_s^{\min }\right)}^2}}}\end{array}\right)---\left(4\right)$

式中分别表示(1)式中X(s)和Y(s)对参数s的一次导数在离散节点X(s_min)的 值；

C3、任一流体网格节点属性判断模块；对任一流体网格节点x_j，定义

$\delta =n_x^j·N_x^{\min }+n_y^j·N_y^{\min }---\left(5\right)$

如果δ≥0，则认为该流体网格节点位于固体外部，如果δ<0，则认为该流体网格节点 位于固体内部。

所述步骤D中，更新流场速度的具体步骤如下：

D1、流固界面力密度求解模块；通过流场次区域速度u′(x_j,t)和δ(x-X(s,t))近似光 滑函数得到的固体边界拉格朗日点上的速度u(X(s_i)等于给定的固体边界的自然速度 来实现流固边界力密度的求解，并将结果输出到文件sfor.txt；其中s为固体边 界离散曲线网格点的初始构型坐标，t为时间，变量下标i和j分别表示固体边界离散曲线网 格第i个节点和流场区域欧拉网格第j个单元；

D2、流场区域速度校正模块；

D2.1、如果流体网格节点x_j位于纯流场区域，该节点速度不需要校正，即

u(x_j,t)＝u′(x_j,t)                     (6)

D2.2、如果流体网格节点x_j位于次流场区域，调用欧拉网格点属性判断模块：

D2.2.1、如果流体网格节点x_j位于固体内部，流体网格节点x_j速度式中u(x_j,t)由 下式更新

u(x_j,t)＝0                (7)

D2.2.2、如果流体网格节点x_j位于固体外部，其速度校正值Δu(x_j,t)为：

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

式中Δs_i为第i段固体边界的面积，式中的C_ji为信息转换矩阵，定义如下：

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

式中，h为流场区域的网格间距，函数φ表示为：

$\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|<0\\ 0,& 2\le \left|r\right|\end{array}\right)---\left(10\right)$

流体网格节点x_j速度u(x_j,t)由下式更新

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

实施例2：如图1-5所示，一种新的投影浸入边界法的速度校正方法，所述方法的具体步 骤如下：

流场区域内包含某一水轮机活动导叶翼型，其翼弦线长度为L，最大厚度为D，计算区域 为15L×10L，坐标原点位于翼弦线的中点。流动雷诺数定义为Re＝ρU_bD/μ＝200，时间步 长Δt＝0.002，方柱边界离散采用间距为0.01L。

S1：网格划分

流场区域(包括流体和水轮机活动导叶翼型区域)在欧拉描述下采用笛卡尔网格划分， 网格为均匀四边形网格，网格间距为h＝0.01L。并将流场区域划分为纯流场区域和次流场区域， 本例中次流场区域为1.4L×1.4L，次流场区域中心为坐标原点。并将次流场区域网格节点坐标 信息输出到文件fcor.txt，固体边界区域在拉格朗日描述下使用适体曲线网格划分，间距为 0.01L，相应网格节点信息文件输出到scor.txt文件。

S2：流场区域速度预测

边值条件：

计算区域左边为进口边界条件，为u＝U_b，v＝0，上下为滑移无穿透边界条件，即u＝U_b， v＝0，右边为出口边界条件，

初值条件：

设定流场区域速度初始值u＝U_b，v＝0，流场压力p＝0。

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

S3：欧拉网格点属性判断模块

S3.1、固体边界离散节点外法线方向求解模块。对任一固体边界离散节点X(s_i)，其相邻 两个节点分别为X(s_i-1)和X(s_i+1)，通过此三点可作一条抛物线近似该段曲线，以s为参数， 该抛物线方程可表示为

$\left(\begin{array}{c}X\left(s\right)=a_x^i{s}^2+b_x^{i}s+c_x^i\\ Y\left(s\right)=a_y^i{s}^2+b_y^{i}s+c_y^i\end{array}\right)---\left(1\right)$

式中为二次抛物线方程(1)的系数，可由此三点的坐标唯一确定。

该段曲线上任意一点的外法线方向可表示为

$\left(\begin{array}{c}{N}_x=\frac{-Y_s}{\sqrt{X_s^2+Y_s^2}}\\ N_y=\frac{X_s}{\sqrt{X_s^2+Y_s^2}}\end{array}\right)---\left(2\right)$

此处，分别表示x,y方向的单位矢量。式中X_s,Y_s分别表示(1)式中X(s)和Y(s)对 参数s求一次导数。

S3.2、任一流体网格节点方向导数求解模块。对任一流体网格节点x_j，求得与之最近的 固体边界离散节点X(s_min)，从该固体边界离散节点X(s_min)到流体网格节点x_j的方向导数 ${\overrightarrow{n}}^j=(n_x^j\overrightarrow{i},n_y^j\overrightarrow{j})$可表示为

$\left(\begin{array}{c}{n}_x^j=\frac{x_j-X\left(s_{\min }\right)}{\sqrt{{(x_j-X\left(s_{\min }\right))}^2+{(y_j-Y\left(s_{\min }\right))}^2}}\\ n_y^j=\frac{y_j-Y\left(s_i^{\min }\right)}{\sqrt{{(x_j-X\left(s_{\min }\right))}^2+{(y_j-Y\left(s_{\min }\right))}^2}}\end{array}\right)---\left(3\right)$

依据式(2)，对固体边界离散节点X(s_min)，其外法线方向可表示为

$\left(\begin{array}{c}{N}_x^{\min }=\frac{-Y_s^{\min }}{\sqrt{{\left(X_s^{\min }\right)}^2+{\left(Y_s^{\min }\right)}^2}}\\ N_y^{\min }=\frac{X_s^i}{\sqrt{\sqrt{{\left(X_s^{\min }\right)}^2+{\left(Y_s^{\min }\right)}^2}}}\end{array}\right)---\left(4\right)$

式中分别表示(1)式中X(s)和Y(s)对参数s的一次导数在离散节点X(s_min)的 值。

S3.3、任一流体网格节点属性判断模块。对任一流体网格节点x_j，定义

$\delta =n_x^j·N_x^{\min }+n_y^j·N_y^{\min }---\left(5\right)$

如果δ≥0，就认为该流体网格节点位于固体外部，如果δ<0，就认为该流体网格节点 位于固体内部。

S4：速度校正计算

S4.1、流固界面作用力密度计算

为满足流固界面无滑移无渗透的一致边界条件，通过流场区域速度u′(x_j,t)和 δ(x-X(s,t))近似光滑函数得到的固体边界拉格朗日点上的速度u(X(s_i)应该等于给定的固体 边界的自然速度来实现流固界面力密度的求解，并将结果输出到文件sfor.txt。 其中s为固体边界离散曲线网格点的初始构型坐标，t为时间，变量下标i和j分别表示固体 边界离散曲线网格第i个节点和流场区域欧拉网格第j个单元。

S4.2：流场区域速度校正。

S4.2.1、如果流体网格节点x_j位于纯流场区域，该节点速度不需要校正，即

u(x_j,t)＝u′(x_j,t)                   (6)

S4.2.2、如果流体网格节点x_j位于次流场区域，调用欧拉网格点属性判断模块。

S4.2.2.1、如果流体网格节点x_j位于固体内部，流体网格节点x_j速度式中u(x_j,t)可 由下式更新

u(x_j,t)＝0                 (7)

S4.2.2.2、如果流体网格节点x_j位于固体外部，其速度校正值Δu(x_j,t)为：

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

式中Δs_i为第i段固体边界的面积，式中的C_ji为信息转换矩阵，定义如下：

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

式中，h为流场区域的网格间距，函数φ可表示为：

$\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|<0\\ 0,& 2\le \left|r\right|\end{array}\right)---\left(10\right)$

流体网格节点x_j速度u(x_j,t)可由下式更新

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

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

图3、图4显示了n＝10000步，即计算时间t＝20时流固边界的作用力密度沿流固边界曲 线的分布。图3为流固界面力密度在x方向的分量沿导叶翼型边界曲线的分布(s_i＝1为 起点，位于翼型的后缘点，沿逆时针方向)，图4为流固界面力密度在y方向的分量沿 导叶翼型边界曲线的分布。

图5为本发明中计算得到的n＝10000步，即计算时间t＝20时流场的速度分布等值线图。

S6：时间推进

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

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