WO2014043847A1 - 求解亚音速流动的反问题的数值方法 - Google Patents

求解亚音速流动的反问题的数值方法 Download PDF

Info

Publication number
WO2014043847A1
WO2014043847A1 PCT/CN2012/081519 CN2012081519W WO2014043847A1 WO 2014043847 A1 WO2014043847 A1 WO 2014043847A1 CN 2012081519 W CN2012081519 W CN 2012081519W WO 2014043847 A1 WO2014043847 A1 WO 2014043847A1
Authority
WO
WIPO (PCT)
Prior art keywords
variable
flow
solving
solid wall
function
Prior art date
Application number
PCT/CN2012/081519
Other languages
English (en)
French (fr)
Inventor
路明
Original Assignee
Lu Ming
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Lu Ming filed Critical Lu Ming
Priority to PCT/CN2012/081519 priority Critical patent/WO2014043847A1/zh
Priority to US13/985,047 priority patent/US9384169B2/en
Publication of WO2014043847A1 publication Critical patent/WO2014043847A1/zh

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/06—Multi-objective optimisation, e.g. Pareto optimisation using simulated annealing [SA], ant colony algorithms or genetic algorithms [GA]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/28—Fuselage, exterior or interior
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Definitions

  • the invention belongs to the field of computational fluid dynamics, and is specifically a numerical method for solving a class of inverse problems.
  • This type of inverse problem refers to the design problem of the geometry of the wall in the subsonic flow field based on the known pressure distribution on the wall.
  • a typical inverse problem is to specify the solid wall profile to meet the pressure distribution requirements by giving a pressure distribution across the solid wall.
  • the traditional methods for solving such problems are based on the Euler plane method.
  • the conjugate equation method (adjoint method) the process is to first estimate the geometry of the undetermined solid wall, then generate a computational grid around it, solve the flow field, and obtain the wall pressure distribution. Next, solve the conjugate equation and obtain information for correcting the geometry of the solid wall. This process is repeated until the final goal is found.
  • the solid wall For a non-viscous flow, regardless of the viscosity of the fluid, the solid wall must be a complete streamline. Each streamline in the flow field has its own flow function. The stream function on a streamline is a constant. So no matter how the geometry of the solid wall changes, its flow function is constant. According to this principle, if such an inverse problem can be solved on the plane of the stream function, it is not necessary to repeatedly generate the calculation grid as in the Euler plane, so that the calculation time can be greatly reduced. In addition, since the flow function represents different flow lines on the flow function plane, there is no fluid passing between the two flow lines. In the direction perpendicular to the streamline, there is no convective flux, which avoids any form of error in the approximation of the convection term.
  • the present invention first provides a coordinate transformation from the Eulerian variable (t, x, y) of the Euler plane to the variable ( ⁇ , ⁇ ) on the plane of the flow function.
  • the flow function dimension [ 2 ⁇ - 1 ]
  • the variable r and the time item t have the same and are introduced as time history entries.
  • the invention begins with an Euler equation on a Euler plane describing a two-dimensional, non-viscous, compressible fluid flow.
  • the variables t, v, p and p are time, flow velocity V in two Cartesian coordinates (- direction of fluid density and pressure.
  • Total energy E and total ⁇ H are
  • Equation (7) A two-dimensional unstable Euler equation, called a flow function, is more simplified than equation (1).
  • the first and fourth terms of the sum vector are 0, which means that the continuous equation and the energy equation have no convection term crossing the mesh interface. This feature not only reduces the popular approximation error, but also makes a solution to the inverse problem of the geometric design of a given solid wall pressure more accurate and faster.
  • the characteristic line equation of equation (7) can be written.
  • the solid wall must be a complete streamline.
  • This streamline can be represented by a value.
  • the idea of a numerical method for solving the inverse problem of geometric design is: Solving the Riemann problem directly on the solid wall of the geometry to be sought. At this time, the solid wall, that is, the local streamline is also the calculation network. Grid interface. The so-called Riemann problem is well known and will not be described here.
  • the direction along the ⁇ direction is the direction across the streamline.
  • Figure 1 shows the interpretation of the Riemann problem of wall flow lines.
  • the left variable 0 £ or the right variable is a shock wave or an expansion wave, which can be approximated by the characteristic line (10, 12).
  • the contact wave represented by the characteristic line (11) is divided into a left intermediate variable Q and a right intermediate variable Q. Since the slope of the characteristic line (11) is 0, it means that the contact wave always develops along the streamline over time.
  • the following task is to solve the Riemann problem and find the original variable ⁇ , ⁇ , ⁇ on the streamline (streamline), and then find the 0 here.
  • the original variable is also an intermediate variable in the Riemann problem.
  • a Riemann problem solver is derived as follows. First, along the characteristic equations (14) and (15), the left and right variables are connected to the intermediate variables, so
  • Equation (17)_(20) is the value of the intermediate variable in the Riemann problem
  • Equations (17) - (25) are called to solve the Riemann problem on the solid wall.
  • the Riemann problem on the solid wall also includes left or right variables.
  • one side is the mirrored variable in the virtual grid for the solid wall on the other side.
  • the subscripts indicate the physical variables p, u, V, ⁇ in the virtual grid, and the subscripts indicate that these physical variables are in the wall grid cell. It is the angle of the wall of the object.
  • Figure 2 shows the relationship between the wall mesh and its mirrored variables, specifically:
  • the above formula is a combined function of the known pressure distribution, that is, after the wall pressure distribution ⁇ is specified, the flow angle in the virtual grid can be solved.
  • the geometric design problem means finding the angle of the solid wall based on a given solid wall pressure distribution.
  • the Newton-Raphson calendar method can be used to find the root of the equation (35) at the given time.
  • the updated value of £ is expressed as
  • the flow angle in the virtual grid can be obtained from t 1 ⁇ .
  • the solid wall angle ⁇ is obtained from the following equation according to the mirror boundary condition.
  • Figure 3 shows an algorithmic flow chart for solving the inverse problem of the geometry design. This method starts with the assumed parameter value of the flow field physical quantity and the solid wall angle. According to the specified pressure distribution on the solid wall, the Riemann problem on the solid wall is solved, and the flow angle is obtained in the virtual grid unit. Then the wall shape is obtained and brought back to the flow field solver to solve the updated flow field. Physical quantity. This process continues until the angle of the solid wall of the calendar has reached convergence.
  • the convergence criterion is the wall angle rather than any physical variable, which means that the wall angle also becomes a flow field variable.
  • the change in the wall boundary does not change the characteristics of the equation.
  • the equation is still a strong hyperbolic type and is still an initial value problem.
  • the current method does not need to generate a mesh around the updated geometry, nor does it need to solve the adjoint equation, which simultaneously obtains a converged flow field physical variable. Solutions and solutions to geometric shapes. For this type of problem, the overall CPU time is reduced by an order of magnitude.
  • This example is about the design problem of solving the solid wall geometry by solving the two-dimensional Euler equation and obtaining the subsonic and non-viscous flow conditions.
  • Distribution a type of inverse problem that solves the geometry of a solid wall.
  • the geometry is first defined, the solid wall pressure distribution is obtained, and then a pressure distribution is used as input.
  • the geometry of the solid wall is obtained by the method of the present invention, and the result should be consistent with the original geometry.
  • the solid wall profile of the nozzle is given by the following function

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Computing Systems (AREA)
  • Operations Research (AREA)
  • Fluid Mechanics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)

Abstract

一种在给定固体壁面压力情况下进行固体壁面几何形状设计的反问题的数值方法，通过坐标变换，将二维欧拉平面上的欧拉方程变换到流函数平面，求解流函数形式的欧拉方程，同时求解固体壁面上的黎曼问题，同时给出流场的物理参数的解和固体壁面的几何形状的解，从而提高计算效率。

Description

说 明 书
求解亚音速流动的反问题的数值方法
1. 技术领域
本发明属于计算流体力学领域， 具体是一种求解一类反问题的数值方法。 这类反问题 是指在亚音速流场中， 根据已知的壁面上的压力分布， 求得壁面的几何形状的设计问题。
2. 背景技术
在空气动力学中， 一类典型的反问题是通过给定固体壁面上的压力分布， 然后设计固 体壁面形状以符合压力分布的要求。 求解这类问题的传统方法都基于欧拉平面的方法。 例 如共轭方程法 (adjoint 法）， 其流程是首先估计这个未被确定的固体壁面的几何形状， 然 后在其周围生成计算网格， 对流场求解,获得壁面压力分布。 下一步， 求解共轭方程， 获 得修正固体壁面几何形状的信息。 这个过程反复进行， 直到找到最终目标为止。
上述对流场的求解， 和目前大多数计算流体力学的计算中一样， 是在笛卡尔坐标下的 欧拉平面中求解欧拉形式的流体控制方程。 在这种框架下， 计算网格必须根据物体的形状 事先生成。 由于该类反问题中的固体几何形状需要不断修正， 计算网格也要随之不断重新 生成。 再加上必须的网格光顺、 正交化处理， 整体求解的计算通常持续较长时间。
对于无粘流， 不考虑流体的粘性， 固体壁面一定是一条完整的流线。 而流场中的每一 条流线都有对应其各自的流函数。 在一条流线上的流函数是常数。 所以无论固体壁面的几 何形状怎样变化， 它的流函数是常数。 根据这个原理， 如果能够在流函数的平面上求解该 类反问题，则不需要像在欧拉平面上那样反复生成计算网格， 因而可以大量减少计算时间。 此外， 由于在流函数平面上， 流函数代表不同的流线， 而二条流线之间是没有流体穿过的。 在与流线垂直方向， 没有对流项通量， 则避免了任何形式对对流项的逼近产生的误差。
尽管在流函数平面上求解该类反问题时表现了如此卓越的优势， 但目前只是应用于超 音速流动中。 此种方法求解超音速流动时， 即方程在空间上向下游方向发展， 不需要考虑 任何下游对上游的影响， 这一切完美地符合超音速流动的特点。 以前， 此种方法成功地应 用于二维超音速流动中的固体壁面的几何形状的设计中。 到目前为止， 应用流函数作为坐 标进行求解几何形状设计的反问题的例子仅限于有势流和线性化的可压缩流。
在工业界有明显的对该类反问题求解的需求。 例如在产品设计的初级阶段， 需要对产 品进行数值模拟研究和快速的几何形状设计， 如机翼、 喷管的外形的设计等。 亚音速也是 工业界最常见的流动状态。 本发明的目的在于提供一个在亚音速、 无粘、 可压缩流场中进 行几何形状设计的反问题的最有效数值方法。 3. 发明内容
本发明首先提供一个从欧拉平面的欧拉变量 (t,x,y)到流函数平面上的变量 (Γ,^ ）的 坐标变换。 流函数平面上， 两个独立的变量， 是流函数 （量纲 [ 2 ^-1]), 代表了不同 流体颗粒在流线上的位置 （量纲 [ ])。 变量 r和时间项 t有相同， 被引入作为时间历遍项。 本发明开始于描述二维、 无粘、 可压缩流体流动的欧拉平面上的欧拉方程
6F_ dG
(1) dt dx dy 其中 f是守恒变量矢量； F, G 是对流项通量在笛卡尔方向 -y上的两个分量。 具体地,
(2)
Figure imgf000003_0004
变量 t, v, p和 p分别是时间、 流动速度 V在两个笛卡尔坐标 （ - 方向上的分 流体密度和压力。 总能 E 和总焓 H分别是
E (3)
Figure imgf000003_0001
上式中的 ； κ是气体比热比。 从坐标 (t,x,j)到 的坐标变换的雅各比矩阵可写为
Figure imgf000003_0002
且它的绝对值 /成为
Figure imgf000003_0003
其中， ^是流动方向角； U，V 被称作流函数几何状态变量， 单位量纲是 [ ^ g- ^， 被 定义为
U = 5x_ V = dy_
(6)
~ θξ' ~ θξ
最终， 二维不稳定欧拉方程 )被转换成流函数平面上的形式， df, d¥. 5G (
■+ - +■ = 0 (7) δτ δλ δξ 其中 f ， 和（^与（2)中的 f， F， G称谓相同， 只是在流函数平面。 具体地，
(8)
Figure imgf000004_0003
加上（5)中 J的定义和 （6 ) 中的定义， 且
V
e = tg (9) 方程 （7) 被称作成流函数形式的二维不稳定欧拉方程， 比方程（1 )更加简化。 和 向量的第一、 四项为 0， 意味着连续方程和能量方程没有对流项穿过网格交界面。 这个 特点不仅减小了对流行的逼近误差， 也使得一类关于给定固体壁面压力的几何形状设计的 反问题的求解更加精确、 快捷。
按照公知的偏微分方程的特征分析方法， 可以写出方程 (7)的的特征线方程
Figure imgf000004_0001
其中， 《为声速。 沿着特征线（10、 12)的特征方程分别为 d ,p +士 pa Vu 2 +ν 2 a ,u = η ϋ； (14)
V
Figure imgf000004_0002
如前所述， 对于无粘流的流动， 固体壁面一定是一条完整的流线。 该流线可以用 的 某个值表示。 对于二维空间的例子， 一个求解几何设计反问题的数值方法的思路是： 在待 求几何形状的固体壁面上直接求解黎曼问题。 此时， 固体壁面， 即当地的流线也是计算网 格交界面。 所谓黎曼问题是公知的， 这里不再叙述
本发明中， 沿着 ^方向是跨过流线的方向。 将方程 (7)写成仅沿着 ^方向的形式， df_ dG(f) , ξ<0,
(16) a + δξ ξ> , 其中， 1^和^分别是一条流线（网格交界面）左侧和右侧变量中的守恒变量。从现在开始， 下标 Z、 7?和 *分别代表黎曼问题中的左、 右和中间变量。 图 1表示了壁面流线的黎曼问题 的解释。 图中表示， 按照黎曼问题的定义， 左侧变量0£或右侧变量 是激波或者是膨胀 波， 可用特征线 （10、 12) 近似表示。 在其中间是中间变量， 被用特征线 （11) 表示的接 触波分为左中间变量 Q 和右中间变量 Q 。 由于特征线 （11) 的斜率是 0， 意味着接触 波一直沿着流线随时间发展。
下面的任务是求解黎曼问题找到网格交界面上 （流线上） 的原始变量 ρ,ί,ν， 再求 得此处的0。 所述的原始变量也是黎曼问题中的中间变量。 一个黎曼问题求解器按照以下 流程导出。 首先沿着特征方程（14)和（15)积分， 将左侧变量和右侧变量与中间变量连接， 于是有
p„+CL- f(¾*,v*) = pL + CL- \uL (17)
P*_CR · {u„v,) = pR—CR · ^{uR,vR (18)
Figure imgf000005_0001
ptR - f(w»,v, PR \ = pR - ^uR,vR PR (20) 其中， = pLaL， CR = pRaR。 并且
1 U + V
f( U,V) = +丄 w lnfv + u + v (21)
2 、 上式中， f(M,v)被称作组合函数， 它可以被认为是速度分量 在流函数平面上的组 合， 与速度有同样的量纲 [ /^|。 公式 (21)有它的极坐标形式 A(V,0)^ (tg0 + cos ^ln((l + sin θ)) + ^ - VlnV , (22) 其中， V = yju2 +v2 , = Vcos0, v = Vsm0 方程（17)_(20) 的解是黎曼问题中的中间 变量的值
Figure imgf000006_0001
对于本发明中涉及的一类几何形状设计的反问题， 需要在固体壁面上求解黎曼问题。 首先， 处理固体壁面边界需要围绕在计算域周围再加上两层网格， 被称为虚拟网格， 它可 以使空间离散算法扩展至边界。 公式 （17) - (25) 被称为求解固体壁面上的黎曼问题。
本发明中， 固体壁面上的黎曼问题同样包括左侧或右侧变量。 但由于固体壁面的特殊 性， 可以假设其中一侧是另一侧针对固体壁面的、 处于虚拟网格中的镜像变量。 下标 表 示物理变量 p, u, V, ^在虚拟网格中， 下标 表示这些物理变量在壁面网格单元中。 是物体壁面的角度。 图 2 表示了壁面网格和其镜像变量之间的关系， 具体是：
V
Figure imgf000006_0002
根据上式以及 (22)— (25)给出的黎曼问题的解， 在固体壁面边界上有
Figure imgf000006_0003
中间变量中的压力即可被认为是物体壁面上的压力 pw， 因而 = *。 （29) 如果物体壁面上的压力已经被给定， 这个给定压力可以直接被用做黎曼问题中的中间 变量中的压力， 所以
P* = PW ° (30) 从上式和（22)， 可以获得
PW =PR + (VR,eL)-f,{VR,0R)), (31)
Figure imgf000007_0001
中 ^,^,^,^Α是边界上已知的量， f是公式（21)和（22)给出的组合函数。 从（31) 导出
2( V
(32)
P 上式是已知压力分布的组合函数， 即在规定了壁面压力分布 ^后， 可以求解出虚拟网 格中流动角度 的大小。
几何形状的设计问题意味着要根据给定的固体壁面压力分布， 找到固体壁面的角度。 这个过程按照下面的步骤建立， 首先定义
tg0L = ε (33)
Figure imgf000007_0002
然后带入到组合 (22)中， 并被进一步 ε的函数，
V。
Figure imgf000007_0003
下面的任务是求解在规定的物理参数 和 ^条件下的方程 F(s) = 0。 解得 后， 再求 。 首先， 方程 (35) 是可微分的， F针对 £的一阶导数是
Figure imgf000007_0004
数值试验表明， 对于无量纲速度 VR≤5 (亚音速区间） 而且 1范围内，
:'(ε)>0 , 这意味着函数 F(f)在此范围内是单调的。 同时， F(-f).F(f)<0， 所以
Newton-Raphson 历遍方法可用来找到方程（35) 在给定 时的根 £。 更新的 £值表示为 虚拟网格中的流动角度可由 t 1^ 求得。 固体壁面角度 ^按镜像边界条件从 下式求得，
Figure imgf000008_0001
图 3 给出了求解几何外形的设计的反问题的算法流程图。这个方法起始于假设的流场 物理量的参数值和固体壁面角度。 按固体壁面上按照的指定的压力分布， 求解固体壁面上 的黎曼问题， 在虚拟网格单元中得出流动角度， 然后壁面形状可获得并被带回到流场求解 器中求解更新流场物理量。这个过程持续进行，直到历遍的固体壁面的角度达到收敛为止。 很明显， 收敛标准是壁面角度而不是任何一个物理变量， 这意味着， 壁面角度也变成了一 个流场变量。 壁面边界的变化并不改变方程的特征， 方程仍然是强双曲线型， 仍然是一个 初值问题。 与应用在欧拉平面上的方法 （如 adjoint法） 相比， 现在的方法既不需要在更 新的几何形状周围生成网格， 也不需要求解 adjoint方程， 它同时获得收敛的流场物理变 量的解和几何形状的解。 针对这类问题， 总体的 CPU时间降低一个量级。
4. 附图说明
图 1 流函数平面上沿 ^方向跨过流线的黎曼问题
图 2 壁面网格和其镜像变量之间的关系
图 3 几何形状设计的反问题的计算流程图
图 4 根据给定压力分布计算的固体壁面的几何形状
(a)欧拉平面上的网格
(b) JST方法计算的固体壁面上压力分布及精解
(c)流函数平面上的网格
(d) 用本发明的方法计算的固体壁面形状
5. 具体实施方式
下面以一个具体实施例子进一步说明本发明的原理。 该例子是关于求解二维欧拉方程 获得亚音速、 无粘流条件下去解决固体壁面几何形状的设计问题， 即已知固体壁面的压力 分布， 求解固体壁面的几何形状的一类反问题。 例子中是一个收缩 -扩张喷管。 例子中先 规定几何形状， 求得固体壁面压力分布， 然后用着个压力分布作为输入， 用本发明的方法 求得固体壁面的几何形状， 结果应该与最初的几何形状一致。 喷管固体壁面外形由下面的 函数给出
y = Hthe l3x2 , _l≤x≤l, (39) 其中， Hi¾取 0.1， 表明收缩喉口高度和入口高度之比为 10% 。 取 6.5以控制喉口长度约 占喷管长度的三分之一。 入口马赫数为 0.5， 是纯亚音速流动。 喷管的几何形状和欧拉平 面上的计算网格如图 4 (a)所示。 计算采用两组网格数目 60x20和 90X30个计算网格单元 在流函数平面。 首先按照公知的 JST方法计算流程完成计算， 获得固体壁面压力分布， 如 图 4(b)所示， 并作为求解几何形状设计的反问题的压力分布的输入。 同时， adjoint方法 将同时应用这个例子以检验本发明的方法的计算效率。 本发明提出的计算一类反问题的具 体步骤是：
(1)初始化流场和固体壁面角度。 流场初始变量的设置 Q° =[ ,^， , ]^及 °,V。均为 初始值， 其中 Q°取无穷远出的值， U °=0; V °=l。 上标 0 表示流动问题在初始时刻， 在 x- 平面和 - 平面， t = 0 ( =0)。 进而可求得守恒变量 f°。 固体壁面初始角度 假设 =0°。 因为流场中的流动角度设为 0， 所以虚拟网格单元中的流动角度 =0°。 流函数平面上的计算网格如图 4 (c) 所示。 无论固体壁面角度如何变化， 计算网格形 状不变。
(2)输入规定的固体壁面的压力分布 pw。
(3)在固体壁面上求解基于给定的压力分布 pw， 按照公式 (30)-(37)，获得 。
(4)计算固体壁面角度。 固体壁面角度计算式根据镜像条件有，
Figure imgf000009_0001
(5)检验固体壁面角度的收敛性。 固体壁面角度 的变化量是两个时刻 ^ 7和 7的差值，
Figure imgf000009_0002
如果 小于收敛标准，流场物理量和固体壁面的角度 +1均达到最终的解;否则， 继 续下一步。 (6)更新流场变量。 这一步的任务是根据更新的固体壁面角度 +1计算流场的物理参数从 [Q 更新到 [Q] 。 这一步可以采用很多公知的方法， 这里不再叙述。 完成这一步后， 获得流场新的物理参数变量 p PR,VR eR：)。
(7)重复步骤 (2)。 图 4(d) 表示了计算的固体几何形状和原始的形状的比较。从结果可以估计计算精度， 计算的固体几何形状和原始的形状吻合良好， 最大误差在之内 0.5%。 最后， 在计算效率方 面， 本发明的方法使用的计算的时间是 adjoint方法的十分之一。 而且同时给出流场的物 理参数的解和固体壁面的几何形状的解， 从而使这类问题达到了最高效率。

Claims

权利要求
1. 一种求解亚音速流动的反问题的数值方法， 其特征是， 包括以下步骤:
1 0 0
(1) 将二维欧拉平面上的欧拉方程经过以雅各比矩阵 J: u cos0 U
Figure imgf000011_0001
V sin^ V 变换，变成为时间 r表示的时间方向、 由流函数 ^和流体颗粒运行的距离 I表示的 两个空间方向所决定的流函数平面上的流函数形式的欧拉方程， 它的形式为： 5fc 5FC dG
■+- +^=o,其中， f;是守恒变量矢量； 和^分别是流函数平面上 δτ δλ δξ 的 I 和 方向上的对流项通量， 而且，
e V
= tg-x\-\, 以上公式中 和 E
Figure imgf000011_0002
分别是流体密度、 压力、 总能； V是流动速度的笛卡尔坐标分量； U，V 是流 函数几何状态变量；
(2) 建立计算网格；
(3) 求解流函数形式的欧拉方程，同时求解固体壁面上的黎曼问题。
2. 根据权利要求 1所述的一种求解亚音速流动的反问题的数值方法， 其中所述的计算网 格是流函数平面上以 和 ^为两个方向的二维直角网格。
3. 根据权利要求 1所述的一种求解亚音速流动的反问题的数值方法， 其中所述的求解流 函数形式的欧拉方程是沿着时间 r方向进行守恒变量 fs的历遍， 以找到其稳定解。
4. 根据权利要求 1所述的一种求解亚音速流动的反问题的数值方法， 其中所述的固体壁 面上的黎曼问题， 其特征是： 以固体壁面为对称， 形成计算域一侧的壁面网格单元中 的变量及其镜像变量， 其中一侧形成以激波或者是膨胀波形式存在的左侧变量或右侧 变量， 在其中间是中间变量； 中间变量又可被分为左中间变量和右中间变量。
5. 根据权利要求 1所述的一种求解亚音速流动的反问题的数值方法， 其中所述的求解固 体壁面上的黎曼问题， 包括以下步骤：
(1) 初始化流场变量参数和固体壁面角度；
(2) 录入规定的分布在物体壁面上的压力值；
(3) 求解黎曼问题并找到物体壁面边界的镜像变量；
(4) 计算固体壁面角度；
(5) 检测固体壁面角度的收敛性， 如果收敛则计算结束；
(6) 更新流场变量参数；
(7) 重复步骤 （3)。
6. 根据权利要求 5所述的求解黎曼问题并找到物体壁面边界的镜像变量，包括以下步骤:
(1)沿着流函数形式的欧拉方程的特征方程积分， 将左侧变量或右侧变量与中间变量 连接， 其中， 左侧变量、 右侧变量、 中间变量由权利要求 4给出；
(2)在中间变量中恢复流动速度的幅值；
(3)求解组合函数 f(M, v)以找到中间变量的流动角；
(4)在中间变量中求解速度分量。
7. 根据权利要求 5所述的求解黎曼问题并找到物体壁面边界的镜像变量， 其中所述的在 中间变量中恢复流动速度的幅值，是利用通过跨过激波的 Rankine-Hugoniot关系和跨 过膨胀波的焓不变关系来获得的。
8. 根据权利要求 5所述的求解黎曼问题并找到物体壁面边界的镜像变量， 其中所述的组 合函数表示为 f(M，v) = iw¾ + +1Μιη ν + Λ/Μ 2 + ν2 在已知压力分布的情况下， 该组合函数表 2 u 2 . 示为 i(VR,0L ) ^ i(VR, 0R) + 2^ ~ PR 其中， pR ,pR ,aR ,VR,eR 是在固体壁面边界上的已知
P 的流动参数； pw是规定的固体壁面的压力分布； 是镜像中的流动角度。
PCT/CN2012/081519 2012-09-18 2012-09-18 求解亚音速流动的反问题的数值方法 WO2014043847A1 (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
PCT/CN2012/081519 WO2014043847A1 (zh) 2012-09-18 2012-09-18 求解亚音速流动的反问题的数值方法
US13/985,047 US9384169B2 (en) 2012-09-18 2012-09-18 Numerical method for solving an inverse problem in subsonic flows

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2012/081519 WO2014043847A1 (zh) 2012-09-18 2012-09-18 求解亚音速流动的反问题的数值方法

Publications (1)

Publication Number Publication Date
WO2014043847A1 true WO2014043847A1 (zh) 2014-03-27

Family

ID=50340511

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2012/081519 WO2014043847A1 (zh) 2012-09-18 2012-09-18 求解亚音速流动的反问题的数值方法

Country Status (2)

Country Link
US (1) US9384169B2 (zh)
WO (1) WO2014043847A1 (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112580436A (zh) * 2020-11-25 2021-03-30 重庆邮电大学 一种基于黎曼流形坐标对齐的脑电信号域适应方法
CN113095006A (zh) * 2021-03-30 2021-07-09 湖南科技大学 一种营造宽薄水幕的缝隙喷口边界流线确定方法
CN117494322A (zh) * 2024-01-02 2024-02-02 中国人民解放军国防科技大学 亚跨超声速流场可控喷管的设计方法、装置、设备和介质

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109241585B (zh) * 2018-08-20 2021-03-23 西北工业大学 一种高低压涡轮过渡流道型面反问题设计方法
CN115122667B (zh) * 2022-01-28 2023-11-10 北京宇航系统工程研究所 自动成型冯.卡门型复合材料斜置网格结构及其成型工装
CN115238397B (zh) * 2022-09-15 2022-12-02 中国人民解放军国防科技大学 高超声速飞行器热环境计算方法、装置和计算机设备
CN117034470B (zh) * 2023-09-08 2024-03-29 北京流体动力科学研究中心 基于高性能数值计算的飞行器外形快速反设计方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20080023946A (ko) * 2006-09-12 2008-03-17 한국지질자원연구원 오일러 해를 이용한 지하공동의 3차원 중력역산 방법 및이를 이용한 3차원 영상화 방법
CN102203782A (zh) * 2010-09-09 2011-09-28 天津空中代码工程应用软件开发有限公司 求解拉格朗日形式的欧拉方程的数值方法
CN102890733A (zh) * 2012-09-18 2013-01-23 天津空中代码工程应用软件开发有限公司 求解亚音速流动的反问题的数值方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5926399A (en) * 1996-03-04 1999-07-20 Beam Technologies, Inc. Method of predicting change in shape of a solid structure

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20080023946A (ko) * 2006-09-12 2008-03-17 한국지질자원연구원 오일러 해를 이용한 지하공동의 3차원 중력역산 방법 및이를 이용한 3차원 영상화 방법
CN102203782A (zh) * 2010-09-09 2011-09-28 天津空中代码工程应用软件开发有限公司 求解拉格朗日形式的欧拉方程的数值方法
CN102890733A (zh) * 2012-09-18 2013-01-23 天津空中代码工程应用软件开发有限公司 求解亚音速流动的反问题的数值方法

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112580436A (zh) * 2020-11-25 2021-03-30 重庆邮电大学 一种基于黎曼流形坐标对齐的脑电信号域适应方法
CN112580436B (zh) * 2020-11-25 2022-05-03 重庆邮电大学 一种基于黎曼流形坐标对齐的脑电信号域适应方法
CN113095006A (zh) * 2021-03-30 2021-07-09 湖南科技大学 一种营造宽薄水幕的缝隙喷口边界流线确定方法
CN113095006B (zh) * 2021-03-30 2022-02-01 湖南科技大学 一种营造宽薄水幕的缝隙喷口边界流线确定方法
CN117494322A (zh) * 2024-01-02 2024-02-02 中国人民解放军国防科技大学 亚跨超声速流场可控喷管的设计方法、装置、设备和介质
CN117494322B (zh) * 2024-01-02 2024-03-22 中国人民解放军国防科技大学 亚跨超声速流场可控喷管的设计方法、装置、设备和介质

Also Published As

Publication number Publication date
US20150186333A1 (en) 2015-07-02
US9384169B2 (en) 2016-07-05

Similar Documents

Publication Publication Date Title
WO2014043847A1 (zh) 求解亚音速流动的反问题的数值方法
WO2012031398A1 (en) Numerical method for simulating subsonic flows based on euler equations in lagrangian formulation
Bourguet et al. Anisotropic Organised Eddy Simulation for the prediction of non-equilibrium turbulent flows around bodies
Hosters et al. Fluid–structure interaction with NURBS-based coupling
CN105677994A (zh) 流体－固体耦合传热的松耦合建模方法
WO2014043846A1 (zh) 求解二维黎曼问题模拟亚音速无粘流的数值方法
CN114112286B (zh) 一种高超声速风洞轴对称型面喷管拟合喉道段设计方法
Jin et al. A unified moving grid gas-kinetic method in Eulerian space for viscous flow computation
Hou et al. Smoothed particle hydrodynamics simulations of flow separation at bends
WO2013078627A1 (zh) 可压缩旋流场的数值模拟方法
Li et al. Coupled simulation of fluid flow and propellant burning surface regression in a solid rocket motor
CN116502338A (zh) 一种通用化的基于线性稳定性理论的工程转捩预测方法
He et al. An efficient nonlinear reduced-order modeling approach for rapid aerodynamic analysis with OpenFOAM
Páscoa et al. A fast iterative inverse method for turbomachinery blade design
Francescatto et al. A semi‐coarsening strategy for unstructured multigrid based on agglomeration
Arias et al. Finite Volume Simulation Of A Flow Over A Naca 0012 Using Jameson, Maccormack, Shu And Tvd Esquemes.
Halanay et al. Existence and approximation for a steady fluid-structure interaction problem using fictitious domain approach with penalization
CN102890733A (zh) 求解亚音速流动的反问题的数值方法
Qin et al. Flow feature aligned grid adaptation
JP2004012248A (ja) 翼形性能の推定方法および翼形性能の推定プログラム
WO2013078626A2 (zh) 不可压缩旋流场的数值模拟中使用的涡量保持技术
Kuwahara et al. Direct simulation of a flow around a subsonic airfoil
Li et al. The Simulation of Wraparound Fins Aerodynamic Characteristics
Orang et al. Computational study of incompressible turbulent flows with method of characteristics
Yu et al. Homotopy continuation of the high-order flux reconstruction/correction procedure via reconstruction (FR/CPR) method for steady flow simulation

Legal Events

Date Code Title Description
WWE Wipo information: entry into national phase

Ref document number: 13985047

Country of ref document: US

121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 12884921

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 12884921

Country of ref document: EP

Kind code of ref document: A1