基于求解定压力反问题的流动几何边界的设计方法与流程

1.技术领域

本发明属于计算流体力学领域，具体是一种求解一类反问题的数值方法。这类反问题是指在亚音速管道流场中，根据已知的壁面上的压力分布，求得管道壁面的几何形状的设计问题。

2.

背景技术：

在空气动力学中，一类典型的反问题是通过给定固体壁面上的压力分布，然后设计固体壁面形状以符合压力分布的要求。求解这类问题的传统方法都基于欧拉平面的方法。例如共轭方程法(adjoint法)，其流程是首先估计这个未被确定的固体壁面的几何形状，然后在其周围生成计算网格，对流场求解,获得壁面压力分布。下一步，求解共轭方程，获得修正固体壁面几何形状的信息。这个过程反复进行，直到找到最终目标为止。

上述对流场的求解，和目前大多数计算流体力学的计算中一样，是在笛卡尔坐标下的欧拉平面中求解欧拉形式的流体控制方程。在这种框架下，计算网格必须根据物体的形状事先生成。由于该类反问题中的固体几何形状需要不断修正，计算网格也要随之不断重新生成。再加上必须的网格光顺、正交化处理，整体求解的计算通常持续较长时间。

对于无粘流，不考虑流体的粘性，固体壁面一定是一条完整的流线。而流场中的每一条流线都有对应其各自的流函数。在一条流线上的流函数是常数。所以无论固体壁面的几何形状怎样变化，它的流函数是常数。根据这个原理，如果能够在流函数的平面上求解该类反问题，则不需要像在欧拉平面上那样反复生成计算网格，因而可以大量减少计算时间。此外，由于在流函数平面上，流函数代表不同的流线，而二条流线之间是没有流体穿过的。在与流线垂直方向，没有对流项通量，则避免了任何形式对对流项的逼近产生的误差。

尽管在流函数平面上求解该类反问题时表现了如此卓越的优势，但目前只是应用于超音速流动中。此种方法求解超音速流动时，即方程在空间上向下游方向发展，不需要考虑任何下游对上游的影响，这一切完美地符合超音速流动的特点。以前，此种方法成功地应用于二维超音速流动中的固体壁面的几何形状的设计中。到目前为止，应用流函数作为坐标进行求解几何形状设计的反问题的例子仅限于有势流和线性化的可压缩流。

在工业界有明显的对该类反问题求解的需求。例如在产品设计的初级阶段，需要对产品进行数值模拟研究和快速的几何形状设计，如机翼、喷管的外形的设计等。亚音速也是工业界最常见的流动状态。本发明的目的在于提供一个在亚音速、无粘、可压缩流场中进 行几何形状设计的反问题的最有效数值方法。

3.

技术实现要素：

本发明首先提供一个从欧拉平面的欧拉变量(t,x,y)到流函数平面上的变量(τ,λ,ξ)的坐标变换。流函数平面上，两个独立的变量，ξ是流函数(量纲[m²·s^-1])，λ代表了不同流体颗粒在流线上的位置(量纲[m])。变量τ和时间项t有相同，被引入作为时间历遍项。本发明开始于描述二维、无粘、可压缩流体流动的欧拉平面上的欧拉方程

其中f是守恒变量矢量；f,g是对流项通量在笛卡尔方向x-y上的两个分量。具体地,

变量t,u,v,ρ和p分别是时间、流动速度v在两个笛卡尔坐标(x-y)方向上的分量、流体密度和压力。总能e和总焓h分别是

和

上式中的γ是气体比热比。从坐标(t,x,y)到(τ,λ,ξ)的坐标变换的雅各比矩阵可写为

且它的绝对值j成为

其中，θ是流动方向角；被称作流函数几何状态变量，单位量纲是[m²s·kg^-1]，被定义为

最终，二维不稳定欧拉方程(1)被转换成流函数平面上的形式，

其中fl,fl和gl与(2)中的f,f,g称谓相同，是在拉格朗日平面。具体地，

h＝0.99，加上(5)中j的定义和(6)中的定义，且

方程(7)被称作成流函数形式的二维不稳定欧拉方程，比方程(1)更加简化。fl和gl向量的第一、四项为0，意味着连续方程和能量方程没有对流项穿过网格交界面。这个特点不仅减小了对流行的逼近误差，也使得一类关于给定固体壁面压力的几何形状设计的反问题的求解更加精确、快捷。

按照公知的偏微分方程的特征分析方法，可以写出方程(7)的的特征线方程

dλ/dτ＝0；(11)

dξ/dτ＝±a/j；(12)

dξ/dτ＝0，(13)

其中，a为声速。沿着特征线(10、12)的特征方程分别为

如前所述，对于无粘流的流动，固体壁面一定是一条完整的流线。该流线可以用ξ的某个值表示。对于二维空间的例子，一个求解几何设计反问题的数值方法的思路是：在待求几何形状的固体壁面上直接求解黎曼问题。此时，固体壁面，即当地的流线也是计算网格交界面。所谓黎曼问题是公知的，这里不再叙述。

本发明中，沿着ξ方向是跨过流线的方向。将方程(7)写成仅沿着ξ方向的形式，

其中,fl和fr分别是一条流线(网格交界面)左侧和右侧变量中的守恒变量。从现在开始，下标l、r和*分别代表黎曼问题中的左、右和中间变量。图1表示了壁面流线的黎曼问题的解释。图中表示，按照黎曼问题的定义，左侧变量ql或右侧变量qr是激波或者是膨胀波，可用特征线(10、12)近似表示。在其中间是中间变量，被用特征线(11)表示的接触波分为左中间变量q*l和右中间变量q*r。由于特征线(11)的斜率是0，意味着接触波一直沿着流线随时间发展。

下面的任务是求解黎曼问题找到网格交界面上(流线上)的原始变量p,ρ,u,v，再求得此处的g。所述的原始变量也是黎曼问题中的中间变量。一个黎曼问题求解器按照以下流程导出。首先沿着特征方程(14)和(15)积分，将左侧变量和右侧变量与中间变量连接，于是有

其中，cl＝ρlal，cr＝ρrar。并且

上式中，f(u,v)被称作组合函数，它可以被认为是速度分量(u,v)在流函数平面上的组 合，与速度有同样的量纲[m/s]。公式(21)有它的极坐标形式

其中，u＝vcosθ，v＝vsinθ。方程(17)-(20)的解是黎曼问题中的中间变量的值

对于本发明中涉及的一类几何形状设计的反问题，需要在固体壁面上求解黎曼问题。首先，处理固体壁面边界需要围绕在计算域周围再加上两层网格，被称为虚拟网格，它可以使空间离散算法扩展至边界。公式(17)-(25)被称为求解固体壁面上的黎曼问题。

本发明中，固体壁面上的黎曼问题同样包括左侧或右侧变量。但由于固体壁面的特殊性，可以假设其中一侧是另一侧针对固体壁面的、处于虚拟网格中的镜像变量。下标l表示物理变量ρ,p,u,v,θ在虚拟网格中，下标r表示这些物理变量在壁面网格单元中。θw是物体壁面的角度。图2表示了壁面网格和其镜像变量之间的关系，具体是：

根据上式以及(22)—(25)给出的黎曼问题的解，在固体壁面边界上有

p*＝pr,(27)

中间变量中的压力即可被认为是物体壁面上的压力pw，因而

pw＝p*。(29)

如果物体壁面上的压力已经被给定，这个给定压力可以直接被用做黎曼问题中的中间变量中的压力，所以

p*＝pw。(30)

从上式和(22)，可以获得

其中pr,ρr,ar,vr,θr是边界上已知的量，f是公式(21)和(22)给出的组合函数。从(31)推导出

上式是已知压力分布的组合函数，即在规定了壁面压力分布pw后，可以求解出虚拟网格中流动角度θl的大小。

几何形状的设计问题意味着要根据给定的固体壁面压力分布，找到固体壁面的角度。这个过程按照下面的步骤建立，首先定义

tgθl＝ε。(33)

虚拟网格中的速度分量按下式求得

然后带入到组合函数(22)中，并被进一步简化为ε的函数，

下面的任务是求解在规定的物理参数ρr,pr,vr,θr和pw条件下的方程f(ε)＝0。解得ε后，再求θl。首先，方程(35)是可微分的，f针对ε的一阶导数是

数值试验表明，对于无量纲速度vr≤5(亚音速区间)而且ε∈[-1,1]的范围内,f′(ε)>0，这意味着函数f(ε)在此范围内是单调的。同时，f(-ε)·f(ε)<0，所以newton-raphson历遍方法可用来找到方程(35)在给定vr时的根ε。更新的ε值表示为

虚拟网格中的流动角度可由θl＝tg^-1(εⁿ⁺¹)求得。固体壁面角度θw按镜像边界条件从下式求得，

图3给出了求解几何外形的设计的反问题的算法流程图。这个方法起始于假设的流场物理量的参数值和固体壁面角度。按固体壁面上按照的指定的压力分布，求解固体壁面上的黎曼问题，在虚拟网格单元中得出流动角度，然后壁面形状可获得并被带回到流场求解器中求解更新流场物理量。这个过程持续进行，直到历遍的固体壁面的角度达到收敛为止。很明显，收敛标准是壁面角度而不是任何一个物理变量，这意味着，壁面角度也变成了一个流场变量。壁面边界的变化并不改变方程的特征，方程仍然是强双曲线型，仍然是一个初值问题。与应用在欧拉平面上的方法(如adjoint法)相比，现在的方法既不需要在更新的几何形状周围生成网格，也不需要求解adjoint方程，它同时获得收敛的流场物理变量的解和几何形状的解。针对这类问题，总体的cpu时间降低一个量级。

4.附图说明

图1流函数平面上沿ξ方向跨过流线的黎曼问题

图2壁面网格和其镜像变量之间的关系

图3几何形状设计的反问题的计算流程图

图4根据给定压力分布计算的固体壁面的几何形状

(a)欧拉平面上的网格

(b)jst方法计算的固体壁面上压力分布及精解

(c)流函数平面上的网格

(d)用本发明的方法计算的固体壁面形状

5.具体实施方式

下面以一个具体实施例子进一步说明本发明的原理。该例子是关于求解二维欧拉方程获得亚音速、无粘流条件下去解决固体壁面几何形状的设计问题，即已知固体壁面的压力分布，求解固体壁面的几何形状的一类反问题。例子中是一个收缩-扩张喷管。例子中先规定几何形状，求得固体壁面压力分布，然后用着个压力分布作为输入，用本发明的方法求得固体壁面的几何形状，结果应该与最初的几何形状一致。喷管固体壁面外形由下面的函数给出

其中，hth取0.1，表明收缩喉口高度和入口高度之比为10％。β取6.5以控制喉口长度约占喷管长度的三分之一。入口马赫数为0.5，是纯亚音速流动。喷管的几何形状和欧拉平面上的计算网格如图4(a)所示。计算采用两组网格数目60x20和90x30个计算网格单元在流函数平面。首先按照公知的jst方法计算流程完成计算，获得固体壁面压力分布，如图4(b)所示，并作为求解几何形状设计的反问题的压力分布的输入。同时，adjoint方法将同时应用这个例子以检验本发明的方法的计算效率。本发明提出的计算一类反问题的具体步骤是：

(1)初始化流场和固体壁面角度。流场初始变量的设置q⁰＝[ρ⁰,u⁰,v⁰,p⁰]^t及均为初始值，其中q⁰取无穷远出的值，上标0表示流动问题在初始时刻，在x-y平面和λ-ξ平面，t＝0(τ＝0)。进而可求得守恒变量f⁰。固体壁面初始角度假设θw＝0°。因为流场中的流动角度设为0，所以虚拟网格单元中的流动角度θl＝0°。流函数平面上的计算网格如图4(c)所示。无论固体壁面角度如何变化，计算网格形状不变。

(2)输入规定的固体壁面的压力分布pw。

(3)在固体壁面上求解基于给定的压力分布pw，按照公式(30)-(37),获得θl。

(4)计算固体壁面角度。固体壁面角度计算式根据镜像条件有，

(5)检验固体壁面角度的收敛性。固体壁面角度θw的变化量是两个时刻n+1和n的差值，

如果δθw小于收敛标准，流场物理量和固体壁面的角度均达到最终的解；否则，继续下一步。

(6)更新流场变量。这一步的任务是根据更新的固体壁面角度计算流场的物理参数从更新到这一步可以采用很多公知的方法，这里不再叙述。完成这一步后，获得流场新的物理参数变量(ρr,pr,vr,θr,)。

(7)重复步骤(2)。

图4(d)表示了计算的固体几何形状和原始的形状的比较。从结果可以估计计算精度，计算的固体几何形状和原始的形状吻合良好，最大误差在之内0.5％。最后，在计算效率方面，本发明的方法使用的计算的时间是adjoint方法的十分之一。而且同时给出流场的物理参数的解和固体壁面的几何形状的解，从而使这类问题达到了最高效率。