CN114117968B - 水气两相流全耦合求解的飞行器水上起降载荷分析方法 - Google Patents

水气两相流全耦合求解的飞行器水上起降载荷分析方法 Download PDF

Info

Publication number
CN114117968B
CN114117968B CN202210080907.8A CN202210080907A CN114117968B CN 114117968 B CN114117968 B CN 114117968B CN 202210080907 A CN202210080907 A CN 202210080907A CN 114117968 B CN114117968 B CN 114117968B
Authority
CN
China
Prior art keywords
equation
water
time
fluid
control body
Prior art date
Legal status (The legal status 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 status listed.)
Active
Application number
CN202210080907.8A
Other languages
English (en)
Other versions
CN114117968A (zh
Inventor
肖天航
支豪林
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
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 Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN202210080907.8A priority Critical patent/CN114117968B/zh
Publication of CN114117968A publication Critical patent/CN114117968A/zh
Application granted granted Critical
Publication of CN114117968B publication Critical patent/CN114117968B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • 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

Landscapes

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

Abstract

本发明公开了一种水气两相流全耦合求解的飞行器水上起降载荷分析方法，首先应用有限体积法将不可压N‑S方程与VOF界面输运方程在任意网格控制单元上进行离散，再根据n时刻的流场值和边界条件，计算对流通量、扩散通量和体积力及其雅克比矩阵，然后通过全隐式数值迭代计算方法得到第n+1时刻的流场值，最后判断当前时间步是否达到计算时长。本发明中，压力方程与动量方程采用虚拟压缩法耦合求解，保证质量和动量守恒，在计算、编程兼容性上有很大优势，同时不可压N‑S方程与VOF方程耦合在单个方程系统中进行全隐式迭代计算，能够大大提高水气两相流数值求解的收敛性和计算效率。

Description

水气两相流全耦合求解的飞行器水上起降载荷分析方法
技术领域
本发明涉及流体力学技术领域，具体涉及一种水气两相流全耦合求解的飞行器水上起降载荷分析方法，基于水气两相流全耦合虚拟压缩有限体积全隐式求解。
背景技术
随着新时代海洋资源的不断开发以及海战的立体化、多层次发展，我国近年来开始从国家层面上重视海洋权益，积极开展水面飞行器的研制工作。水面飞行器需要在水和空气两相介质情况下起降飞行，其对传统的飞行器性能分析方法提出了挑战。此外，通用飞行器水面迫降适航标准的制定也需要行之有效的分析水气两相流的方法。目前针对水面飞行器水面起降、通用飞行器的水上迫降的性能评估主要依靠物理水池模型试验实现，然而由于设备和场地的限制，试验模型的尺寸与实机相比差异过大，该尺寸效应使得试验数据无法直接应用于实机的指导上，同时，试验也面临着人力、财力资源消耗巨大、周期过长等问题。
近年来，随着计算机技术的飞速发展，计算流体力学(CFD)得到了蓬勃发展，CFD技术已成为水动力、气动力学分析的基础性关键技术。由于CFD技术的巨大应用潜力，比拟于物理水池，数值水池的概念应运而生，取代传统水池模型试验以分析水面飞行器水面起降、通用飞行器的水上迫降性能。与传统水池模型试验相比，数值水池具有用时短、成本低、获取信息丰富、无尺度限制等相当突出的优势，在航空航天、船舶海洋等领域有着广泛的应用前景。
数值水池主要是求解水气两相不可压流Navier-Stokes(N-S)方程来获取流场信息。由于不可压缩流体的假设，N-S方程中的质量方程和动量方程是解耦的，一般要采用压力-速度修正方式建立起两者的联系，因此无法保证质量和动量同时守恒。其次，该方式无法采用成熟的双曲型方程时间推进方式求解，在编程、计算兼容性等方面带来了不便。另一方面，目前成熟的两相流数值模拟方法通过引入一个VOF输运方程实现对水气交界面的捕捉来实现，其求解的是每个网格控制单元内水所占的体积分数，N-S方程中流体属性根据该体积分数来确定。如今常用的求解方式是将N-S方程与VOF方程进行解耦求解，这样一定程度上割裂了两者的关系，影响了计算精度和收敛性。
发明内容
为了克服现有技术的不足，本发明提出一种水气两相流全耦合求解的飞行器水上起降载荷分析方法，通过建立压强和密度的虚拟关系，在质量方程中引入压力的虚拟时间导数项，将原本是椭圆型的连续方程变成双曲型方程，即可通过时间推进方法迭代求解，保证质量和动量守恒。同时将N-S方程与VOF方向相耦合，构建统一的全隐式数值求解方程组，提高收敛性和精度，从而指导飞行器在水上起降或迫降。
本发明的技术方案为：一种水气两相流全耦合求解的飞行器水上起降载荷分析方法，包括以下步骤：
S10、基于工程对象几何模型，生成流体计算网格；
S20、设置数值计算的流体属性、边界条件、最大模拟步数和时间步长；
S30、基于流体计算网格控制单元的几何信息(面积、体积等)；
S40、采用虚拟压缩法在N-S方程中建立密度与压力的联系：
Figure GDA0003624281350000021
ρ为密度，p为压力，τ为虚拟时间项，β为虚拟压缩参数，一般取5～15。这样可以将原本椭圆型的连续方程变成双曲型方程，即可通过时间推进法求解，使得质量和动量都能守恒。
此时引入虚拟时间导数项后的不可压N-S方程和VOF方程具体表达式为：
Figure GDA0003624281350000022
Figure GDA0003624281350000023
Figure GDA0003624281350000024
Figure GDA0003624281350000025
u为流体速度向量，μ为流体粘度，fb为体积力向量，fst为表面张力向量，T为流体温度，k为流体热传导系数，Re为雷诺数，Pr为普朗特数，α为水体积分数，t为物理时间。
S50、运用有限体积法将上述N-S方程与VOF方程写成任意欧拉-拉格朗日的积分形式：
Figure GDA0003624281350000031
采用有限体积法在任意网格控制单元上离散，同时物理时间方向上采用k阶后向差分离散，每个物理时间步伪时间推进求解，伪时间用一阶后向差分，可以得到任意控制体i处离散后的方程为：
Figure GDA0003624281350000032
Figure GDA0003624281350000033
S60、根据流体属性、边界条件和n时刻的流场值计算对流通量、扩散通量、体积力及其雅克比矩阵；
S70、采用全隐式数值迭代计算方法，得到n+1时刻的流场值；
S80、判断当前时间步数是否达到最大模拟步数，若是则结束计算，反之返回步骤S60。
进一步地，步骤S50中：
Q为流场变量，F为对流通量，Fv为扩散通量，Sb为体积力源项，SST为表面张力源项，Ω(t)为t时刻的控制体体积，Γ为预处理矩阵，避免两相流时密度差异过大使特征系统变的刚性，同时也是为了将积分方程表达为守恒形式。
上述变量表达式为：
Figure GDA0003624281350000034
Figure GDA0003624281350000035
Figure GDA0003624281350000041
其中，
Figure GDA0003624281350000042
表示控制体边界
Figure GDA0003624281350000043
的法向运动速度，
Figure GDA0003624281350000044
和n分别为控制体边界
Figure GDA0003624281350000045
的运动速度和单位法向矢量，n＝(nx ny nz)。θ为控制体边界流体法向速度，θ＝u·n＝unx+vny+wnz。Δρ＝ρwater-ρair表示水和气体的密度差。Fr为弗劳德数。
其中j为控制体i的邻居控制体，ij表示控制体i和控制体j的交接面，nface为当前控制体边界面的数量。离散方程中包括两个时间层的推进，一个是物理时间层t、n，是真实物理意义层面的时间；另一个是伪时间层τ、m，其作用在于将每个物理时间步当作一个伪时间层面的定常问题来推进求解。因此，对非定常问题的求解，在每个物理时间步内都嵌套一个伪时间步的迭代循环；伪时间步迭代求解收敛后的结果，即认为是第n个物理时间层的非定常流场，随即进入下一个物理时间步的求解循环。
此外，φ为后向差分系数，RESGCL为几何守恒定律残差项，
Figure GDA0003624281350000046
进一步地，步骤S60具体为：
不可压N-S方程的对流通量采用Roe迎风格式进行计算，扩散通量采用中心差分格式计算；VOF方程的对流通量采用Modified HRIC格式计算。
进一步地，步骤S70具体为：
将离散方程构造成全隐式线性方程组：
Figure GDA0003624281350000047
该方程可以近似为大型稀疏系数矩阵的线性方程组Ax＝b，可采用成熟且高效的迭代求解方法，如LU-SGS、广义最小残差法(GMRES)求解得到n+1时刻的流场值。
进一步地，步骤S20中的边界条件包括入口边界条件和出口边界条件，基于特征线理论，由对流通量雅克比矩阵的特征值正负情况判定流场信息传入或传出计算域，并由特征变量判定边界上的值，具体为：
对流通量雅克比矩阵可以表示成特征矩阵的形式：
Γ-1A＝XΛX-1
定义流场变量转换为特征变量的形式为：
C＝X-1Q，Γ-1AδQ＝XΛX-1δQ＝XΛδC
边界面上的值根据特征变量不变的原则进行确定，即：
Figure GDA0003624281350000051
其中下标b表示边界值，in表示内场值，∞表示来流值。
本发明相对于现有技术的有益效果是：本发明在不可压流N-S方程的连续性方程中，借助虚拟压缩法建立密度和压力的联系，将原本是椭圆型的连续方程变成双曲型方程，即可通过时间推进方法迭代求解，需要引入压力-速度修正，保证质量和动量守恒，以及编程和计算的统一性。在此基础上，又将N-S方程与VOF方向相耦合，构建统一的全隐式数值求解方程组，大大提高收敛性和精度，从而指导飞行器在水上起降或迫降。
附图说明
图1是本发明一实施例提供的流程图；
图2是本发明一实施例提供的VOF方程对流通量格式示意图；
图3是本发明一实施例提供的精度验证算例；
图4是发明一实施例提供的水球入水算例过程；
图5是发明一实施例提供的两栖飞机波浪水面高速滑行算例姿态变化过程。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明提供一种水气两相流全耦合求解的飞行器水上起降载荷分析方法，如图1所示，包括以下步骤S10-S80：
S10、基于工程对象几何模型，生成流体计算网格。本发明中，采用无量纲形式处理网格节点坐标，计算网格单元可以是任意类型，如四面体、六面体、多面体和切割体等。
S20、设置数值计算需要的相关参数和条件，具体包括：
(1)流体属性
本发明中，流体属性包括密度ρ、粘度μ、热传导系数k、定压比热容Cp，重力加速度g，并根据流体属性确定常规的无量纲参数，具体为：
雷诺数：
Figure GDA0003624281350000061
普朗特数：
Figure GDA0003624281350000062
弗劳德数：
Figure GDA0003624281350000063
韦伯数：
Figure GDA0003624281350000064
其中，L为参考长度，U为参考速度，都由计算工况确定。σ为表面张力系数，一般取常温下水的表面张力系数为：σ＝0.072N/m。
(2)边界条件
本发明中，边界条件包括入口边界、出口边界、壁面和对称面，入口边界需要设置来流速度、压力和温度、水平面位置，出口边界需要设置压力和水平面位置。
(3)最大模拟步数和时间步长；
模拟时长由最大模拟步数和时间步长确定。
S30、基于流体计算网格控制单元的几何信息，包括控制体体积、质心坐标，控制体边界面面积、面法矢、面左右侧单元距离，以备后续通量计算中使用。
S40、采用虚拟压缩法在N-S方程中建立密度与压力的联系：
Figure GDA0003624281350000065
ρ为密度，p为压力，τ为虚拟时间项，β为虚拟压缩参数，一般取5～15。这样可以将原本椭圆型的连续方程变成双曲型方程，即可通过时间推进法求解，使得质量和动量都能守恒。
此时引入虚拟时间导数项后的不可压N-S方程和VOF方程具体表达式为：
Figure GDA0003624281350000071
Figure GDA0003624281350000072
Figure GDA0003624281350000073
Figure GDA0003624281350000074
u为流体速度向量，fb为体积力向量，fst为表面张力向量，T为流体温度，α为水体积分数，t为物理时间。
S50、运用有限体积法将上述N-S方程与VOF方程写成任意欧拉-拉格朗日的积分形式：
Figure GDA0003624281350000075
其中，Q为流场变量，F为对流通量，Fv为扩散通量，Sb为体积力源项，SST为表面张力源项，Ω(t)为t时刻的控制体体积，Γ为预处理矩阵，避免两相流时密度差异过大使特征系统变的刚性，同时也是为了将积分方程表达为守恒形式。
上述变量表达式为：
Figure GDA0003624281350000076
Figure GDA0003624281350000077
Figure GDA0003624281350000081
其中，
Figure GDA0003624281350000082
表示控制体边界
Figure GDA0003624281350000083
的法向运动速度，
Figure GDA0003624281350000084
和n分别为控制体边界
Figure GDA0003624281350000085
的运动速度和单位法向矢量，n＝(nx ny nz)。θ为控制体边界流体法向速度，θ＝u·n＝unx+vny+wnz。Δρ＝ρwater-ρair表示水和气体的密度差。Fr为弗劳德数。
采用有限体积法在任意网格控制单元上离散，同时物理时间方向上采用k阶后向差分离散，每个物理时间步伪时间推进求解，伪时间用一阶后向差分，可以得到任意控制体i处离散后的方程为：
Figure GDA0003624281350000086
Figure GDA0003624281350000087
其中j为控制体i的邻居控制体，ij表示控制体i和控制体j的交接面，nface为当前控制体边界面的数量。离散方程中包括两个时间层的推进，一个是物理时间层t、n，是真实物理意义层面的时间；另一个是伪时间层τ、m，其作用在于将每个物理时间步当作一个伪时间层面的定常问题来推进求解。因此，对非定常问题的求解，在每个物理时间步内都嵌套一个伪时间步的迭代循环；伪时间步迭代求解收敛后的结果，即认为是第n个物理时间层的非定常流场，随即进入下一个物理时间步的求解循环。
此外，φ为后向差分系数，RESGCL为几何守恒定律残差项，
Figure GDA0003624281350000088
S60、根据流体属性、边界条件和n时刻的流场值计算对流通量、扩散通量、体积力及其雅克比矩阵，具体包括：
(1)不可压N-S方程的对流通量通过Roe迎风格式计算：
Figure GDA0003624281350000089
式中L、R分别为控制体交界ij面两侧的流场值，Aij为面ij对流通量雅可比矩阵：
Figure GDA0003624281350000091
Γ为预处理矩阵。将Γ-1A写成，Γ-1A＝XΛX-1，Λ为对角矩阵，对角元素为特征值，X-1,X分别为方程系统的左右特征向量，具体为：
Λ＝diag(λ1 λ2 λ3 λ4 λ5 λ6)，
Figure GDA0003624281350000092
Figure GDA0003624281350000093
Figure GDA0003624281350000101
其中，Θ＝2cc+c-，a和b为两单位矢量，满足：a×b＝n，a·b＝n·a＝n·b＝0。
进一步地，在计算入口、出口边界面上的对流通量值时，基于特征线理论，由对流通量雅克比矩阵的特征值正负情况判定流场信息传入或传出计算域，并由特征变量判定边界上的值，具体为：
对流通量雅克比矩阵可以表示成特征矩阵的形式：
Γ-1A＝XΛX-1
定义流场变量转换为特征变量的形式为：
C＝X-1Q，Γ-1AδQ＝XΛX-1δQ＝XΛδC
边界面上的值根据特征变量不变的原则进行确定，即：
Figure GDA0003624281350000102
其中下标b表示边界值，in表示内场值，∞表示来流值。这样可以保证全隐式数值计算的收敛性。
(2)不可压N-S方程的扩散通量通过中心差分格式计算：
扩散通量的计算，控制体交界面处的速度、温度的梯度通过将左右控制体的梯度平均并沿rij修正得到，
Figure GDA0003624281350000103
其中
Figure GDA0003624281350000111
平均梯度值，
Figure GDA0003624281350000112
为由控制体i的中心指向控制体j的中心的单位方向矢量，左右控制体内速度、温度的梯度由最小二乘法或高斯格林法重构得到。
(3)VOF方程的对流通量通过Modified HRIC格式计算：
Modified High Resolution Interface Capturing(HRIC)格式相比于传统的中心差分或迎风格式，能够更为准确的进行自由液面的捕捉。该格式是一种复合NVD格式，包含了迎风与顺风差分的非线性组合。首先，定义一个体积函数的标准化的单元值
Figure GDA0003624281350000113
Figure GDA0003624281350000114
如图2所示，其中D为顺风单元，C为中心单元，U为迎风单元。得到
Figure GDA0003624281350000115
后，标准化的面值可以如下计算：
Figure GDA0003624281350000116
这里如果直接使用
Figure GDA0003624281350000117
会使得交界面出现褶皱，这里引入ULTIMATE QUICKEST格式，基于单元面法矢和交界面法矢的角度θ对
Figure GDA0003624281350000118
进行修正：
Figure GDA0003624281350000119
Figure GDA00036242813500001110
Figure GDA00036242813500001111
这里与常规HRIC格式不同，无需根据当地CFL数进行修正。最后，得到面f上的体积函数：
Figure GDA00036242813500001112
完成VOF对流通量的计算。
S70、将计算完通量的离散方程表示为如下全隐式格式：
Figure GDA0003624281350000121
该方程可以近似为大型稀疏系数矩阵的线性方程组Ax＝b，可采用成熟且高效的迭代求解方法，如LU-SGS、广义最小残差法(GMRES)求解得到n+1时刻的流场值。
S80、判断当前时间步数是否达到最大模拟步数，若是则结束计算，反之返回步骤S60，进入下一步时间循环计算，直到达到最大模拟步数，完成水气两相不可压流数值求解。
实施例1：
以一个三维楔形体水箱入水问题验证本发明的精度。本实施例为某一物理水池模型试验工况，水箱重600kg，触水速度为6m/s，水和空气的流体属性取常温时候的值，时间步长设为0.001s，最大模拟步数为5000步，模拟时长为5.0s，水箱的姿态受到六自由度模型控制。整个入水过程如图3所示，本发明提出的方法可以很好地捕捉物体与水面接触后引起的水面变化。从最后水箱受到的过载时历曲线与实验值和商用软件计算值的对比可知，本发明的精度较高，误差仅为1.8％。
实施例2：
如图4所示，本实施例模拟了一个水球入水过程，在初始时刻位于水面上方1m处释放一个半径为0.1m的水球，使其在重力场作用下落入水中。水和空气的流体属性取常温时候的值，时间步长取0.001s，最大模拟步数为1000步。图4很好地展示了水球在重力场下的变形、对水面的冲击、与水面融合、产生喷溅的整个过程，表明本发明提出的方法能够较好地模拟水气两相流问题。
实施例3：
本实施例模拟了一大型水上飞机全尺寸构型在波浪水面条件下水上滑行过程，飞机中49t，俯仰惯性矩为2.52×106kg·m2，运动过程受六自由度模型控制。滑行速度为50km/h，波长120m，波高2m。水和空气的流体属性取常温时候的值，时间步长取0.001s，最大模拟步数为5000步。图5展示了一个波浪周期内飞机姿态的变化情况，可以看到，随着波浪沿机尾方向传播，机头部分首先受到冲击，飞机受到抬头力矩，导致其开始抬头，同时，随着波浪传播造成的液面升高，飞机也开始抬升高度。直到波峰经过飞机重心位置之后，机头部分脱离水面，此时主要飞机后半段受力，飞机受到低头力矩，俯仰角开始变小。随着波面移动，飞机的升沉位移也变小。本实施例表明本发明提出的方法能够有效地处理复杂构型的水气两相流数值模拟问题。

Claims (4)

1.一种水气两相流全耦合求解的飞行器水上起降载荷分析方法，其特征在于，包括以下步骤：
S10、基于工程对象几何模型，生成流体计算网格；
S20、设置数值计算的流体属性、边界条件、最大模拟步数和时间步长；
S30、基于流体计算网格控制单元的几何信息；
S40、采用虚拟压缩法在N-S方程中建立密度与压力的联系；
Figure FDA0003624281340000011
ρ为密度，p为压力，τ为虚拟时间项，β为虚拟压缩参数，取5～15，将原本椭圆型的连续方程变成双曲型方程，即可通过时间推进法求解，使得质量和动量都能守恒；
此时引入虚拟时间导数项后的不可压N-S方程和VOF方程具体表达式为：
Figure FDA0003624281340000012
Figure FDA0003624281340000013
Figure FDA0003624281340000014
Figure FDA0003624281340000015
u为流体速度向量，fb为体积力向量，fst为表面张力向量，T为流体温度，α为水体积分数，t为物理时间，k为流体热传导系数，Re为雷诺数，Pr为普朗特数；
S50、运用有限体积法将N-S方程与VOF方程在网格控制单元上离散成代数方程组；
Figure FDA0003624281340000016
其中，Q为流场变量，F为对流通量，Fv为扩散通量，Sb为体积力源项，SST为表面张力源项，Ω(t)为t时刻的控制体体积，Γ为预处理矩阵；
上述变量表达式为：
Figure FDA0003624281340000021
Figure FDA0003624281340000022
Figure FDA0003624281340000023
其中，
Figure FDA0003624281340000024
表示控制体边界
Figure FDA0003624281340000025
的法向运动速度，
Figure FDA0003624281340000026
和n分别为控制体边界
Figure FDA0003624281340000027
的运动速度和单位法向矢量，n＝(nx ny nz)，θ为控制体边界流体法向速度，θ＝u·n＝unx+vny+wnz，Δρ＝ρwater-ρair表示水和气体的密度差，Fr为弗劳德数；
采用有限体积法在任意网格控制单元上离散，同时物理时间方向上采用k阶后向差分离散，每个物理时间步伪时间推进求解，伪时间用一阶后向差分，可以得到任意控制体i处离散后的方程为：
Figure FDA0003624281340000028
Figure FDA0003624281340000029
其中j为控制体i的邻居控制体，ij表示控制体i和控制体j的交接面，nface为当前控制体边界面的数量；离散方程中包括两个时间层的推进，一个是物理时间层t、n，是真实物理意义层面的时间；另一个是伪时间层τ、m，其作用在于将每个物理时间步当作一个伪时间层面的定常问题来推进求解；因此，对非定常问题的求解，在每个物理时间步内都嵌套一个伪时间步的迭代循环；伪时间步迭代求解收敛后的结果，即认为是第n个物理时间层的非定常流场，随即进入下一个物理时间步的求解循环；φ为后向差分系数，RESGCL为几何守恒定律残差项，
Figure FDA0003624281340000031
S60、根据流体属性、边界条件和n时刻的流场值计算对流通量、扩散通量、体积力及其雅克比矩阵；
S70、采用全隐式数值迭代计算方法，得到n+1时刻的流场值；
S80、判断当前时间步数是否达到最大模拟步数，若是则结束计算，根据得到的流场值指导飞行器起降或迫降，反之返回步骤S60。
2.根据权利要求1所述的水气两相流全耦合求解的飞行器水上起降载荷分析方法，其特征在于，所述步骤S20中的流体属性包括密度、粘度、热传导系数、重力加速度。
3.根据权利要求1述的水气两相流全耦合求解的飞行器水上起降载荷分析方法，其特征在于，步骤S60中不可压N-S方程的对流通量采用Roe迎风格式进行计算，扩散通量采用中心差分格式计算；VOF方程的对流通量采用Modified HRIC格式计算。
4.根据权利要求1述的水气两相流全耦合求解的飞行器水上起降载荷分析方法，其特征在于，步骤S70将离散方程构造成全隐式线性方程组：
Figure FDA0003624281340000032
将该方程近似为大型稀疏系数矩阵的线性方程组Ax＝b，采用迭代求解方法，得到n+1时刻的流场值。
CN202210080907.8A 2022-01-24 2022-01-24 水气两相流全耦合求解的飞行器水上起降载荷分析方法 Active CN114117968B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210080907.8A CN114117968B (zh) 2022-01-24 2022-01-24 水气两相流全耦合求解的飞行器水上起降载荷分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210080907.8A CN114117968B (zh) 2022-01-24 2022-01-24 水气两相流全耦合求解的飞行器水上起降载荷分析方法

Publications (2)

Publication Number Publication Date
CN114117968A CN114117968A (zh) 2022-03-01
CN114117968B true CN114117968B (zh) 2022-07-12

Family

ID=80361230

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210080907.8A Active CN114117968B (zh) 2022-01-24 2022-01-24 水气两相流全耦合求解的飞行器水上起降载荷分析方法

Country Status (1)

Country Link
CN (1) CN114117968B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114580196B (zh) * 2022-03-18 2022-08-16 中国科学院力学研究所 一种切割自由面水翼空泡长度预测方法
CN114912199B (zh) * 2022-05-18 2024-09-20 南京航空航天大学 一种基于解耦分析策略的翼伞最大开伞动载快速预测方法
CN114880971B (zh) * 2022-07-13 2022-09-20 中国空气动力研究与发展中心计算空气动力研究所 一种计算流体力学软件采用的隐式方法
CN117290642B (zh) * 2022-10-28 2024-09-17 国家电投集团科学技术研究院有限公司 基于Newton-Raphson求解器的热工水力模型的耦合方法、装置及设备
CN117010288B (zh) * 2023-06-06 2024-08-23 西安数峰信息科技有限责任公司 基于MHT框架的密度基LU-SGS和Krylov子空间求解器
CN116757123B (zh) * 2023-08-11 2023-11-07 南京航空航天大学 基于二维半理论的水面飞行器结构着水载荷预报方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108984874A (zh) * 2018-07-02 2018-12-11 中国航空发动机研究院 获取不可压缩流动的流场的数值模拟方法
CN111339714A (zh) * 2020-02-17 2020-06-26 珠江水利委员会珠江水利科学研究院 基于FVCOM和OpenFOAM模型的多尺度水动力耦合方法
CN113191095A (zh) * 2021-04-02 2021-07-30 南京航空航天大学 一种发动机吞水风险的评估方法及装置
CN113515805A (zh) * 2021-04-07 2021-10-19 南京航空航天大学 水上飞机波浪水面起降水气两相流数值仿真方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112989680B (zh) * 2021-05-14 2021-07-16 中国空气动力研究与发展中心计算空气动力研究所 减少网格使用量的fvfd远场积分边界条件计算方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108984874A (zh) * 2018-07-02 2018-12-11 中国航空发动机研究院 获取不可压缩流动的流场的数值模拟方法
CN111339714A (zh) * 2020-02-17 2020-06-26 珠江水利委员会珠江水利科学研究院 基于FVCOM和OpenFOAM模型的多尺度水动力耦合方法
CN113191095A (zh) * 2021-04-02 2021-07-30 南京航空航天大学 一种发动机吞水风险的评估方法及装置
CN113515805A (zh) * 2021-04-07 2021-10-19 南京航空航天大学 水上飞机波浪水面起降水气两相流数值仿真方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Numerical Simulation of Three Dimensional Tank Impacting on Wavy Water;Yutong Jin 等;《AIAA Aviation 2019 Forum》;20190621;第1-14页 *
飞行器水载荷结构完整性数值模拟现状与展望-Part I：水上迫降和水上漂浮;童明波 等;《航空学报》;20210525;第42卷(第5期);第524530-1-524530-35页 *

Also Published As

Publication number Publication date
CN114117968A (zh) 2022-03-01

Similar Documents

Publication Publication Date Title
CN114117968B (zh) 水气两相流全耦合求解的飞行器水上起降载荷分析方法
Fu et al. An assessment of computational fluid dynamics predictions of the hydrodynamics of high-speed planing craft in calm water and waves
CN110955991B (zh) 一种界面双向数据交换的流固耦合计算方法
Tahara et al. Computational fluid dynamics–based optimization of a surface combatant
Wang et al. Multi-body separation simulation with an improved general mesh deformation method
CN114168796B (zh) 一种建立飞行器高空气动力数据库的方法
Khazaee et al. Hydrodynamic evaluation of a planing hull in calm water using RANS and Savitsky's method
CN113092065A (zh) 一种潜降式网箱水动力特性计算的分析方法
CN115329689A (zh) 基于伪非定常时间推进的复杂湍流流动高效率计算方法
Guan et al. Automatic optimal design of self-righting deck of USV based on combined optimization strategy
Zhang et al. Moving particle semi-implicit method coupled with finite element method for hydroelastic responses of floating structures in waves
Zhang et al. Machine learning methods for turbulence modeling in subsonic flows over airfoils
Viola et al. A CFD-based wing sail optimisation method coupled to a VPP
CN110705184B (zh) 一种反应堆堆芯精细化数值求解的虚拟体积力动量源法
Anandhanarayanan et al. Development and applications of a gridfree kinetic upwind solver to multi-body configurations
CN116227073A (zh) 一种研究水下发射航行体肩空泡发展溃灭的方法
CN109635370A (zh) 开裂式阻力方向舵静气动弹性特性分析方法
CN116070424A (zh) 一种铅铋快堆铅池热分层降阶分析方法
CN113673007B (zh) 一种基于sph的双向波浪中船舶耐波性预报方法
Lamberson et al. High-Fidelity Aeroelastic Simulations with HPCMP CREATETM-AV Kestrel
Vedam et al. Evaluation of Gradient and Curvature-Based Adaptive Mesh Refinement for Viscous Transonic Flows
Sinha et al. Novel cfd techniques for in-cylinder flows on tetrahedral grids
Liu et al. Numerical analysis of wet-deck slamming characteristics for trimaran section with different main-hull profiles
Miyaji et al. On Accuracy of Prediction of Flutter Boundaries on Unstructured Grids
Li et al. A review of the numerical strategies for solving ship hydroelasticity based on CFD-FEM technology

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant