CN109726433B - 基于曲面边界条件的三维无粘低速绕流的数值模拟方法 - Google Patents

基于曲面边界条件的三维无粘低速绕流的数值模拟方法 Download PDF

Info

Publication number
CN109726433B
CN109726433B CN201811454186.2A CN201811454186A CN109726433B CN 109726433 B CN109726433 B CN 109726433B CN 201811454186 A CN201811454186 A CN 201811454186A CN 109726433 B CN109726433 B CN 109726433B
Authority
CN
China
Prior art keywords
point
boundary
gaussian
plane
unit
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
CN201811454186.2A
Other languages
English (en)
Other versions
CN109726433A (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.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
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 University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN201811454186.2A priority Critical patent/CN109726433B/zh
Publication of CN109726433A publication Critical patent/CN109726433A/zh
Application granted granted Critical
Publication of CN109726433B publication Critical patent/CN109726433B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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

  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明属于三维流体力学数值求解技术领域，涉及一种基于曲面边界条件的三维无粘低速绕流的数值模拟方法。本发明提出的一种基于曲面边界条件的三维无粘低速绕流的数值模拟方法，通过把平面高斯积分点投影到物面，然后通过物面法向计算固壁边界条件来提高边界条件的精度，从而在较少的网格情况下，达到精度要求，减少计算量。

Description

基于曲面边界条件的三维无粘低速绕流的数值模拟方法
技术领域
本发明属于三维流体力学数值求解技术领域，涉及一种基于曲面边界条件的三维无粘低速绕流的数值模拟方法。
背景技术
飞机在飞行过程中会受到空气压力的作用，船舶在航行过程中会受到水的作用力，房屋、桥梁在大风中会受到风力的作用等等，这些可以统称为流体动力学。人们一般采用理论计算法、实验测量法和数值计算法三种方法对上述现象进行分析。理论计算法只能对简单的结构进行分析，实验测量法成本太高，随着计算机技术的发展，数值模拟方法显得越来越重要，尤其是复杂结构，如飞机、桥梁、船舶等结构三维绕流的分析。
随着科学技术的进步，流体动力学对数值算法的精度提出了更高的要求，因此需要高精度的数值模拟方法。相比传统的有限差分、有限体积、有限元法，间断Galerkin有限元法具有以下优点：易于实现高精度、灵活处理间断问题、有利于实现并行算法，因此有很好的应用前景和工程实用价值。然而高精度的数值算法需要高精度的边界条件支撑，尽管有些研究者们采用非结构网格的方式来提高边界的拟合度，如采用四面体网格来更好的拟合曲面边界。然而这始终是一种近似的方法，边界的拟合程度取决于四面体网格的大小，网格越小拟合的越好，但会带来大量的网格，无疑会导致巨大的计算量。因此需要从算法本身出发，提出一种高精度的曲面边界计算方法，在较少的网格情况下，达到精度要求，从而减少计算量。
发明内容
针对上述存在问题或不足，为解决边界拟合导致的巨大计算量的问题，从算法本身出发获得一种高精度三维无粘低速绕流的数值模拟方法；本发明提供了一种基于曲面边界条件的三维无粘低速绕流的数值模拟方法。
一种基于曲面边界条件的三维无粘低速绕流的数值模拟方法，包括以下步骤：
A.将球型结构进行建模，然后建立流体计算域；
B.对步骤A所建流体计算域采用四面体网格进行剖分，转化为离散空间模型；
C.利用间断Galerkin有限元法，将三维无粘低速绕流的控制方程在步骤B所得每一个网格上进行空间离散，获得一个关于时间微分的有限元方程；
D.采用物面高斯积分点上的物面法向计算固壁边界条件，结合HLLC数值通量和出入口边界条件计算步骤C获得的有限元方程中的各积分项，得到时间微分的有限元方程；
E.对D步骤所得有限元方程进行时间离散，获得迭代方程；
F.对每一个由四面体网格剖分后得到的四面体单元，先给定初值，然后对步骤E获得的迭代方程进行循环迭代，直至满足迭代终止条件，获得整个计算域的场分布。
本发明提出的一种基于曲面边界条件的三维无粘低速绕流的数值模拟方法，通过把平面高斯积分点投影到物面，然后通过物面法向计算固壁边界条件来提高边界条件的精度，从而在较少的网格情况下，达到精度要求，减少计算量。
附图说明
图1是本发明流程图；
图2是四面体网格非曲面边界条件高斯积分点法向示意图；
图3是平面点投影到物面点示意图；
图4是四面体网格曲面边界条件高斯积分点法向示意图；
图5是实施例的流体计算域模型剖面图；
图6是实施例场分布剖视图；
图7是非曲面边界条件三维无粘低速绕流数值模拟方法场分布剖视图。
具体实施方式
下面结合附图和实施例来详细说明本发明的技术方案。
参照附图1，一种基于曲面边界条件的三维无粘低速绕流的数值模拟方法，包括以下步骤：
A.将球型结构进行建模，然后建立流体计算域；
建立球结构的几何模型，然后建立流体计算域，其结构剖视图如图5所示。
B.对步骤A所建流体计算域采用四面体网格进行剖分，转化为离散空间模型；
采用四面体网格剖分步骤A所建流体计算域，剖分后的计算域被人为分割为多个三维四面体网格，从而将连续的几何空间转化为离散的网格空间。
C.利用间断Galerkin有限元法，将三维无粘低速绕流的控制方程在步骤B所得每一个网格上进行空间离散，获得一个关于时间微分的有限元方程；
对于无粘绕流问题，我们需要求解的是如下三维守恒形式的欧拉方程
Figure BDA0001887368020000021
式中U为守恒变量，
Figure BDA0001887368020000022
为无粘通量张量，其具体形式如下：
Figure BDA0001887368020000031
式中ρ,p分别为密度和压力；x,y,z,分别为笛卡尔坐标系下的坐标分量；u,v,w分别为笛卡尔坐标系下的速度分量；E为总能。(1)式实则是5个方程的组合，为了表示方便我们用h(h＝1,2,3,4,5)表示(1)式中第几个方程以及U和
Figure BDA0001887368020000032
中的第几个分量。
对于间断Galerkin有限元法，变量在单元内的分布采用下面的多项式近似表达
Figure BDA0001887368020000033
φj(x,y,z)表示插值基函数，在这里我们选择正交基函数。N表示插值基函数的个数。(1)式两端乘以测试函数φi(x,y,z)(i＝0,…N)，然后在单元Ω内积分并且把(3)式代入，可以得到(1)式中第h个方程的Galerkin方法弱形式
Figure BDA0001887368020000034
其中，
Figure BDA0001887368020000035
为单元Ω的边界，
Figure BDA0001887368020000036
为单元边界的外单位法向量。
令
Figure BDA0001887368020000037
Figure BDA0001887368020000038
(4)式最后简化为
Figure BDA0001887368020000039
其中
Figure BDA00018873680200000310
Mh为质量矩阵，其元素为mij，只与单元坐标及单元类型相关。对(1)式中每一个方程执行上述(3)-(7)式同样的操作。最后得到一个关于时间微分的有限元方程
Figure BDA00018873680200000311
其中u＝[u1 u2 u3 u4 u5]T，RHS＝[RHS1 RHS2 RHS3 RHS4 RHS5]T
Figure BDA0001887368020000041
D.采用物面高斯积分点上的物面法向计算固壁边界条件，结合HLLC数值通量和出入口边界条件计算步骤C获得的有限元方程中的各积分项，得到时间微分的有限元方程；
对于(6)式中的面积分和体积分采用高斯数值积分进行计算。对于面积分项
Figure BDA0001887368020000042
由于单元边界左右两边的值不同，所以
Figure BDA0001887368020000043
通量需要采用数值通量格式代替，这里我们采用HLLC通量，该通量格式是一种已经公开的格式，这里不再详细阐述，下面给出具体形式
Figure BDA0001887368020000044
其中l,r分别表示本单元和相邻单元。SL，SM，SR表示黎曼边界的不同位置。
从(11)式可以看出在进行面的高斯积分时，需要用到本单元和相邻单元在面上高斯积分点的值。对于内部单元可以直接计算，但是对于边界单元，其相邻单元是不存在的因此我们需要构建相应的边界条件。最终根据边界条件给出虚构相邻单元在边界面高斯点上的值。对于出入口边界，我们采用远场边界条件，这是一种公知的边界条件，这里不再详细阐述。固壁边界条件一般是附着在物面上的，对于固壁边界，在无粘情况下，速度为无穿透条件，其表达式为
Figure BDA0001887368020000045
其中
Figure BDA0001887368020000046
表示笛卡尔坐标系下单位矢量分量，固壁边界单元面上每一个高斯积分点上的值都满足上述方程。由于流体计算域是采用四面体网格离散的，因此一般在固壁边界上，单元边界面是一个平面。如果采用平面边界条件，单元边界面外单位法向量
Figure BDA0001887368020000047
在各个高斯积分点的方向是一样的，如图2所示。
本发明提出曲面边界条件来代替平面边界条件，提高边界条件的精度。首先需要通过坐标投影把平面高斯点投影到物面上，获得物面高斯点的坐标。如图3所示，A点是平面上一点，B点是A点投影到物面上所对应的点。点2，3，4为四面体网格在物面上的三个顶点，点5，6，7为四面体边界边在物面上的中点，这些点的坐标信息(x,y,z)可以通过一般网格生成程序获得，这里不再详细阐述。我们采用下式来计算投影点B的坐标
Figure BDA0001887368020000051
Figure BDA0001887368020000052
Figure BDA0001887368020000053
其中，ξ2(xA,yA,zA)，ξ3(xA,yA,zA)，ξ4(xA,yA,zA)是点2，3，4所对应的体积坐标。xA,yA,zA是A点的坐标分量。对于平面高斯点，我们采用上式投影到物面从而获得物面高斯点的坐标，最终可以根据该点物面曲率等信息获得物面高斯点的单位法向量。因此面外单位法向量
Figure BDA0001887368020000056
在每一个高斯积分点的方向是不一样的，如图4所示。在计算(12)式时，在每个高斯积分点上，需要采用物面外单位法向量。因此(12)式变成
Figure BDA0001887368020000054
其中上下标i表示第i个高斯积分点，
Figure BDA0001887368020000055
表示第i个高斯积分点所对应的物面外单位法向量。
E.对D步骤所得有限元方程进行时间离散，获得迭代方程；
时间离散上我们采用大家所熟知的二阶龙格库塔方法，显示二阶龙格库塔方法如下：
Figure BDA0001887368020000061
其中k表示时间步。上述方程，是一个随时间步迭代的方程，可以通过前一个时刻k的值计算下一个时刻k+1的值。
F.对每一个由四面体网格剖分后得到的四面体单元，先给定初值，然后对步骤E获得的迭代方程进行循环迭代，直至满足迭代终止条件，获得整个计算域的场分布。
根据实际问题给定所有四面体单元的初值，按照(17)式计算所有四面体单元当前时刻的值，然后再把当前时刻的值当成初值，继续计算下一时刻的值，以此反复循环迭代，直至计算结果收敛。最后根据(3)式计算每个单元的场分布，给出整个计算域的场分布。
图5示出了实施例的流体计算域模型剖视图；图6示出了实施例场分布剖视图；图7示出了非曲面边界条件三维无粘低速绕流数值模拟方法场分布剖视图。对比图6、图7的场分布可以看出具体实施例的场分布更加具有对称性，这更加接近真实情况。

Claims (1)

1.一种基于曲面边界条件的三维无粘低速绕流的数值模拟方法，包括以下步骤：
A.将球型结构进行建模，然后建立流体计算域；
B.对步骤A所建流体计算域采用四面体网格进行剖分，转化为离散空间模型；
C.利用间断Galerkin有限元法，将三维无粘低速绕流的控制方程在步骤B所得每一个网格上进行空间离散，获得一个关于时间微分的有限元方程；
D.采用物面高斯积分点上的物面法向计算固壁边界条件，结合HLLC数值通量和出入口边界条件计算步骤C获得的有限元方程中的各积分项，得到时间微分的有限元方程；
对于面积分和体积分采用高斯数值积分进行计算；对于面积分项
Figure FDA0003681740390000011
其中，
Figure FDA0003681740390000012
为单元Ω的边界，
Figure FDA0003681740390000013
为单元边界的外单位法向量，φj表示插值基函数；由于单元边界左右两边的值不同，所以
Figure FDA0003681740390000014
通量需要采用数值通量格式代替，这里采用HLLC通量，具体形式：
Figure FDA0003681740390000015
其中l,r分别表示本单元和相邻单元；SL，SM，SR表示黎曼边界的不同位置；
在进行面的高斯积分时，需要用到本单元和相邻单元在面上高斯积分点的值；对于内部单元直接计算；对于边界单元，其相邻单元不存在，因此构建相应的边界条件，最终根据构建的相应边界条件给出虚构相邻单元在边界面高斯点上的值；对于出入口边界，采用远场边界条件；对于固壁边界，在无粘情况下，速度为无穿透条件，其表达式为
Figure FDA0003681740390000016
其中ρ为密度，p为压力；x,y,z,分别为笛卡尔坐标系下的坐标分量；u,v,w分别为笛卡尔坐标系下的速度分量，
Figure FDA0003681740390000017
表示笛卡尔坐标系下单位矢量分量，固壁边界单元面上每一个高斯积分点上的值都满足上述方程；由于流体计算域是采用四面体网格离散的，因此在固壁边界上，单元边界面是一个平面；如果采用平面边界条件，单元边界面外单位法向量
Figure FDA0003681740390000018
在各个高斯积分点的方向相同；
首先通过坐标投影把平面高斯点投影到物面上，获得物面高斯点的坐标，A点是平面上一点，B点是A点投影到物面上所对应的点；点2，3，4为四面体网格在物面上的三个顶点，点5为四面体边界边2-3在物面上的中点，点6为四面体边界边2-4在物面上的中点，点7为四面体边界边3-4在物面上的中点，这些点的坐标信息x,y,z通过网格生成程序获得，采用下式来计算投影点B的坐标；
Figure FDA0003681740390000021
Figure FDA0003681740390000022
Figure FDA0003681740390000023
其中，ξ2(xA,yA,zA)，ξ3(xA,yA,zA)，ξ4(xA,yA,zA)是点2，3，4所对应的体积坐标；xA,yA,zA是A点的坐标分量；对于平面高斯点，采用上式投影到物面从而获得物面高斯点的坐标，最终根据该点物面曲率信息获得物面高斯点的单位法向量；面外单位法向量
Figure FDA0003681740390000026
在每一个高斯积分点的方向不同，在计算(12)式时，在每个高斯积分点上，采用物面外单位法向量，因此(12)式变成
Figure FDA0003681740390000024
上下标i表示第i个高斯积分点，
Figure FDA0003681740390000025
表示第i个高斯积分点所对应的物面外单位法向量；
E.对D步骤所得有限元方程进行时间离散，获得迭代方程；
F.对每一个由四面体网格剖分后得到的四面体单元，先给定初值，然后对步骤E获得的迭代方程进行循环迭代，直至满足迭代终止条件，获得整个计算域的场分布。
CN201811454186.2A 2018-11-30 2018-11-30 基于曲面边界条件的三维无粘低速绕流的数值模拟方法 Active CN109726433B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811454186.2A CN109726433B (zh) 2018-11-30 2018-11-30 基于曲面边界条件的三维无粘低速绕流的数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811454186.2A CN109726433B (zh) 2018-11-30 2018-11-30 基于曲面边界条件的三维无粘低速绕流的数值模拟方法

Publications (2)

Publication Number Publication Date
CN109726433A CN109726433A (zh) 2019-05-07
CN109726433B true CN109726433B (zh) 2022-10-14

Family

ID=66295208

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811454186.2A Active CN109726433B (zh) 2018-11-30 2018-11-30 基于曲面边界条件的三维无粘低速绕流的数值模拟方法

Country Status (1)

Country Link
CN (1) CN109726433B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110377959B (zh) * 2019-06-21 2023-07-07 上海数设科技有限公司 确定飞机加筋平板对应的有限元及稳定性校核方法和装置
CN110909511B (zh) * 2019-11-22 2022-10-14 电子科技大学 一种无曲面体积分的无粘低速绕流数值模拟方法
CN112330679B (zh) * 2020-11-03 2023-08-01 中国科学院大学 一种面向数控加工的自由曲面的分割方法及系统
CN113723017B (zh) * 2021-07-14 2024-01-02 中国航发沈阳发动机研究所 一种确定航空发动机旋转盘壁面温度场的方法及装置
CN117252131B (zh) * 2023-11-20 2024-03-01 深圳十沣科技有限公司 一种适用于薄壁结构物的数值仿真方法及装置

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107045488A (zh) * 2017-01-21 2017-08-15 三峡大学 一种基于Bezier曲面四边形参数方程的曲面边界元计算方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2378444B1 (de) * 2010-04-13 2015-07-29 CST-Computer Simulation Technology AG Verfahren, Vorrichtung und Computerprogrammprodukt zur Bestimmung eines elektromagnetischen Nahfeldes einer Feldanregungsquelle eines elektrischen Systems
CN107273625A (zh) * 2017-06-22 2017-10-20 西南科技大学 一种三维不可压缩非定常n‑s方程有限元数值求解方法
CN107515982B (zh) * 2017-08-22 2020-08-11 电子科技大学 一种三维力学有限元模态分析中的接触分析方法
CN107577857B (zh) * 2017-08-28 2020-12-29 电子科技大学 一种基于热辐射边界条件的三维有限元模拟方法
CN108052738B (zh) * 2017-12-13 2021-10-15 电子科技大学 色散媒质的高阶局部无条件稳定时域间断伽辽金分析方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107045488A (zh) * 2017-01-21 2017-08-15 三峡大学 一种基于Bezier曲面四边形参数方程的曲面边界元计算方法

Also Published As

Publication number Publication date
CN109726433A (zh) 2019-05-07

Similar Documents

Publication Publication Date Title
CN109726433B (zh) 基于曲面边界条件的三维无粘低速绕流的数值模拟方法
US10296672B2 (en) Generating inviscid and viscous fluid-flow simulations over an aircraft surface using a fluid-flow mesh
CN112016167B (zh) 基于仿真和优化耦合的飞行器气动外形设计方法及系统
CN109726465B (zh) 基于非结构曲边网格的三维无粘低速绕流的数值模拟方法
CN113505443B (zh) 一种任意外形的三维绕流问题自适应笛卡尔网格生成方法
CN112100835B (zh) 一种适用于复杂流动的高效高精度翼型绕流数值模拟方法
CN108763683B (zh) 一种三角函数框架下新weno格式构造方法
CN112613206A (zh) 一种基于各向异性体调和场的附面层网格生成方法
CN111079326B (zh) 二维各向异性网格单元度量张量场光滑化方法
Zhang et al. Rotor airfoil aerodynamic design method and wind tunnel test verification
CN115859781A (zh) 基于注意力和卷积神经网络编解码器的流场预测方法
US20220414282A1 (en) Boundary layer mesh generation method based on anisotropic volume harmonic field
Lu et al. Flow simulation system based on high order space-time extension of flux reconstruction method
CN110765694B (zh) 一种基于简化型浅水方程组的城市地表水流数值模拟方法
Murman et al. A vortex wake capturing method for potential flow calculations
Coirier et al. A Cartesian, cell-based approach for adaptively-refined solutions of the Euler and Navier-Stokes equations
CN110909511B (zh) 一种无曲面体积分的无粘低速绕流数值模拟方法
CN111008492B (zh) 一种基于无雅克比矩阵的高阶单元欧拉方程数值模拟方法
Duan et al. High order FR/CPR method for overset meshes
CN117436317B (zh) 基于海上风电桩基的波流载荷仿真计算方法、系统和设备
Sadrehaghighi Dynamic & Adaptive Meshing
Wang et al. Meshfree simulation of flow around airfoil using different turbulent models
Tarsia Morisco et al. Extension of the Vertex-Centered Mixed-Element-Volume MUSCL scheme to mixed-element meshes
Ghoreyshi et al. Grid Quality and Resolution Effects on the Aerodynamic Modeling of Parachute Canopies
Makino et al. Aerodynamic analysis of NASA common research model by block-structured cartesian mesh

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