CN111859645B - 冲击波求解的改进musl格式物质点法 - Google Patents

冲击波求解的改进musl格式物质点法 Download PDF

Info

Publication number
CN111859645B
CN111859645B CN202010657277.7A CN202010657277A CN111859645B CN 111859645 B CN111859645 B CN 111859645B CN 202010657277 A CN202010657277 A CN 202010657277A CN 111859645 B CN111859645 B CN 111859645B
Authority
CN
China
Prior art keywords
point
momentum
flow field
shock wave
grid
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
CN202010657277.7A
Other languages
English (en)
Other versions
CN111859645A (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 Science and Technology
Original Assignee
Nanjing University of Science and Technology
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 Science and Technology filed Critical Nanjing University of Science and Technology
Priority to CN202010657277.7A priority Critical patent/CN111859645B/zh
Publication of CN111859645A publication Critical patent/CN111859645A/zh
Application granted granted Critical
Publication of CN111859645B publication Critical patent/CN111859645B/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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • 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)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提出了一种冲击波求解的改进MUSL格式物质点法，首先改进传统MUSL格式中的动量修正步骤；然后建立冲击波流场模型，对冲击波流场进行离散化处理；再给定冲击波流场的初始条件与边界条件，设置计算总长；接着基于稳定性条件计算当前计算步的时间步长；基于冲击波求解的改进MUSL格式物质点法求解冲击波流场；最后对冲击波流场进行可视化处理。本发明提出的方法能避免欧拉法中对流项处理困难的问题和拉格朗日法中的网格畸变问题，相比于传统的MUSL格式，本发明可有效地减少物质点穿越网格时的噪声，对形函数的精度要求更低，鲁棒性和稳定性更强，计算效率更高，计算能耗更少，可为工程中冲击波的求解提供一种新方法。

Description

冲击波求解的改进MUSL格式物质点法
技术领域
本发明属于冲击波流场求解方法技术领域，特别是一种冲击波求解的改进MUSL格式物质点法。
背景技术
冲击波是一种不连续的峰在流体介质中的传播，导致流场的压力、温度、密度等物理性质发生阶跃式变化。其在工程中广泛地存在，例如高速列车运行时形成列车冲击波、矿山***、防爆车辆等，可能造成人员的伤亡和设备的损坏。随着计算机技术的发展，CFD技术逐渐被应用于冲击波的仿真模拟，通过数值计算冲击波的传播过程，得到其对流场的扰动作用，进而可为其作用对象的气动载荷计算提供计算参数。此外，冲击波问题是CFD的经典算例之一，波前和波后流场的密度、压力等参数都呈现间断性，这些间断问题的求解一直是CFD发展的难点和核心问题。因此，探究一种新型的冲击波求解方法具有广泛的工程价值。
物质点法是一种混合了拉格朗日法和欧拉法的数值计算方法，已经被广泛地应用于固体结构中的高速冲击问题和***问题求解。该方法用物质点来离散物体，相关物理参数均由物质点携带，背景网格仅用于积分动量方程和求解梯度，因而只需简单地生成覆盖计算域的结构网格即可。物质点法在一定程度上摆脱了网格的约束，避免了大变形造成的网格畸变问题，能很好地处理材料断裂失效、结构破碎等现象，并且不会造成质量损失，适用于求解极端工况。MUSL格式是一种带有动量修正的计算格式，物质点法中传统的MUSL格式中基于当前时间步的形函数将下一时间步的物质点动量映射到网格节点上，得到修正后的下一时刻网格节点动量，用修正后的网格节点动量和修正前的网格节点质量求出网格节点速度。由于此时物质点位置已发生改变，其形函数也应随之变化，然传统的MUSL格式未考虑该变化，当物质点穿越网格时，该修正方法就会产生较大误差，进而影响计算的精度甚至导致计算不收敛。
发明内容
本发明的目的在于提供一种冲击波求解的改进MUSL格式物质点法，为实际工程中预测冲击波的影响提供方法。
实现本发明目的的技术解决方案为：
一种冲击波求解的改进MUSL格式物质点法，包括以下步骤：
步骤1、改进传统MUSL格式中的动量修正步骤：基于当前时间步的形函数将下一时间步的物质点动量映射到网格节点上，得到修正后的下一时刻网格节点动量，取消传统MUSL格式中的动量修正步，基于下一时间步的物质点位置计算形函数，同时将下一时间步的物质点质量和动量映射到网格节点上；
步骤2、建立冲击波流场模型，对冲击波流场进行离散化处理：划分背景网格并布设物质点；
步骤3、给定冲击波流场的初始条件与边界条件，设置计算总长；
步骤4、基于稳定性条件计算当前计算步的时间步长；
步骤5、基于冲击波求解的改进MUSL格式物质点法求解冲击波流场；
步骤6、对冲击波流场进行可视化处理：输出流场的密度、压力参数。
一种基于改进MUSL格式物质点法的冲击波求解***，包括物质点动量修正模块、冲击波流场模型建立模块、计算总长设置模块、时间步长计算模块、冲击波流场计算模块、可视化处理模块；
所述物质点动量修正模块用于基于当前时间步的形函数将下一时间步的物质点动量映射到网格节点上，得到修正后的下一时刻网格节点动量，取消传统MUSL格式中的动量修正步，基于下一时间步的物质点位置计算形函数，同时将下一时间步的物质点质量和动量映射到网格节点上；所述冲击波流场模型建立模块用于、建立冲击波流场模型，对冲击波流场进行离散化处理；所述计算总长设置模块用于给定冲击波流场的初始条件与边界条件，设置计算总长；所述时间步长计算模块用于计算当前计算步的时间步长；冲击波流场计算模块用于基于冲击波求解的改进MUSL格式物质点法求解冲击波流场；所述可视化处理模块用于对冲击波流场进行可视化处理：输出流场的密度、压力参数。
本发明与现有技术相比，其显著优点是：
(1)利用物质点法求解冲击波问题，生成覆盖计算域的规则结构网格，大大减少了网格划分的工作量和时间；单个物质点的质量在计算中始终不变，自动满足质量守恒定律，无需求解质量方程；压力、密度等物理参数均储存在物质点上，降低其数值耗散；在单个计算步中采用拉格朗日法求解，避免了欧拉法中对流项处理困难的问题，在两个计算步之间抛弃变形的背景网格，重新建立物质点和网格之间的映射关系，避免了拉格朗日法中的网格畸变问题。
(2)该修正方法改进了传统MUSL格式中的动量修正步，避免了传统MUSL格式中形函数比节点动量落后一个时间步的弊端，不仅修正网格节点动量，还修正网格节点质量，减少了物质点穿越网格造成的误差，可更精确地捕捉冲击波波面前后的间断特征，能很好地抑制振荡，并能提高计算收敛速度。
附图说明
图1为本发明冲击波物质点求解方法的流程图。
图2为Riemann问题示意图。
图3为本发明和传统MUSL格式对Riemann问题的求解结果图。
图4为本发明和传统MUSL格式对Riemann问题的解析解的对比图。
具体实施方式
下面结合附图及具体实施例对本发明做进一步的介绍。
步骤1、改进传统MUSL格式中的动量修正步骤，基于当前时间步的形函数将下一时间步的物质点动量映射到网格节点上，得到修正后的下一时刻网格节点动量，取消传统MUSL格式中的动量修正步，基于下一时间步的物质点位置计算形函数，同时将下一时间步的物质点质量和动量映射到网格节点上，得到修正后的网格节点质量和动量为：
Figure BDA0002577215680000031
Figure BDA0002577215680000032
其中
Figure BDA0002577215680000033
为网格节点i在下一计算步中的质量，/>
Figure BDA0002577215680000034
为网格节点i在下一计算步中的动量，mp为物质点质量，/>
Figure BDA0002577215680000035
为物质点p在下一计算步中的速度，/>
Figure BDA0002577215680000036
为下一计算步中网格节点i与物质点p之间的映射函数。
步骤2、建立冲击波流场模型，对冲击波流场进行离散化处理，对流场的边界进行几何清洗并建立拓扑，用均匀的结构网格覆盖整个计算域生成背景网格，在边界外侧至少生成一层网格作为虚网格，并在冲击波的计算范围内均匀地布设物质点。
步骤3、给定冲击波流场的初始条件与边界条件，设置计算总长，给定初始时刻各个物质点的密度
Figure BDA0002577215680000037
压力/>
Figure BDA0002577215680000038
速度/>
Figure BDA0002577215680000039
计算出初始时刻各个物质点的体积/>
Figure BDA00025772156800000310
质量mp、内能/>
Figure BDA0002577215680000041
Figure BDA0002577215680000042
Figure BDA0002577215680000043
Figure BDA0002577215680000044
其中VΩ是整个冲击波流场的体积，np为物质点总数，γ为气体比热比，空气取1.4。
在背景网格上施加边界条件，对于固壁边界来说，边界处的网格和边界外的虚网格速度、动量、节点力始终为0。
步骤4、基于稳定性条件计算当前计算步的时间步长，设置的库朗数，计算出当前计算步的时间步长为：
Figure BDA0002577215680000045
Figure BDA0002577215680000046
其中Δtt为当前计算步的时间步长，CCFL为库朗数，
Figure BDA0002577215680000047
为物质点p当前时刻的声速，
Figure BDA0002577215680000048
分别为物质点p在当前计算步中沿x、y、z方向的速度，/>
Figure BDA0002577215680000049
为物质点p在当前计算步中的密度，/>
Figure BDA00025772156800000410
为物质点p在当前计算步中的压力，L为网格长度。
步骤5、基于冲击波求解的改进MUSL格式物质点法求解冲击波流场，包括：
5.1：将物质点参数映射到网格节点上，得到网格节点的质量和动量为：
Figure BDA00025772156800000411
Figure BDA00025772156800000412
其中
Figure BDA00025772156800000413
为网格节点i在当前计算步中的质量，/>
Figure BDA00025772156800000414
为网格节点i在当前计算步中的动量，/>
Figure BDA0002577215680000051
为物质点p在当前计算步中的速度，/>
Figure BDA0002577215680000052
为当前计算步中网格节点i与物质点p之间的映射函数。
5.2：计算网格节点力为：
Figure BDA0002577215680000053
其中fi t为网格节点i在当前计算步中的力值，
Figure BDA0002577215680000054
为当前计算步中网格节点i与物质点p之间的映射函数导数。
5.3：积分动量方程为：
Figure BDA0002577215680000055
其中
Figure BDA0002577215680000056
为网格节点i在当前计算步中的动量，/>
Figure BDA0002577215680000057
为网格节点i在下一个计算步中的动量。
5.4：计算得到下一时刻物质点的位置和速度为：
Figure BDA0002577215680000058
Figure BDA0002577215680000059
其中
Figure BDA00025772156800000510
为物质点p在当前计算步中的速度，/>
Figure BDA00025772156800000511
为物质点p在下一个计算步中的速度，ni为网格节点数，/>
Figure BDA00025772156800000512
为物质点p在当前计算步中的位置，/>
Figure BDA00025772156800000513
为物质点p在下一个计算步中的位置。
5.5：取消传统MUSL格式中修正网格节点动量的步骤，再次将物质点的质量和动量映射到网格节点上，修正网格节点的质量和动量为：
Figure BDA00025772156800000514
Figure BDA00025772156800000515
5.6、更新物质点的密度为：
Figure BDA0002577215680000061
/>
Figure BDA0002577215680000062
Figure BDA0002577215680000063
其中
Figure BDA0002577215680000064
为物质点p在下一个计算步中的应变增量，它是一个方阵，行列为j和k，
Figure BDA0002577215680000065
为物质点p在当前计算步中的密度，/>
Figure BDA0002577215680000066
为物质点p在下一个计算步中的密度。
5.7：利用理想气体状态方程更新物质点的压力为：
Figure BDA0002577215680000067
其中
Figure BDA0002577215680000068
为物质点p在下一个计算步中的内能。
5.8：判断是否已计算至终点时刻，若否则令t＝t+1并转入步骤4，再次计算时间步长，若是则结束本步骤进入下一步。
步骤6、对冲击波流场进行可视化处理，存储网格节点上的密度等参数的时程数据，得到流场内冲击波的传播过程。
基于上述的冲击波求解的改进MUSL格式物质点法和计算机程序，本发明还提出了一种基于改进MUSL格式物质点法的冲击波求解***，具体包括物质点动量修正模块、冲击波流场模型建立模块、计算总长设置模块、时间步长计算模块、冲击波流场计算模块、可视化处理模块；
所述物质点动量修正模块用于基于当前时间步的形函数将下一时间步的物质点动量映射到网格节点上，得到修正后的下一时刻网格节点动量，取消传统MUSL格式中的动量修正步，基于下一时间步的物质点位置计算形函数，同时将下一时间步的物质点质量和动量映射到网格节点上，具体过程见上述步骤1；所述冲击波流场模型建立模块用于、建立冲击波流场模型，对冲击波流场进行离散化处理，具体过程见上述步骤2；所述计算总长设置模块用于给定冲击波流场的初始条件与边界条件，设置计算总长，具体过程见上述步骤3；所述时间步长计算模块用于计算当前计算步的时间步长，具体过程见上述步骤4；冲击波流场计算模块用于基于冲击波求解的改进MUSL格式物质点法求解冲击波流场，具体过程见上述步骤5；所述可视化处理模块用于对冲击波流场进行可视化处理：输出流场的密度、压力参数，具体过程见上述步骤6；由于程序处理的过程同上述方法的处理过程，本发明不再进行赘述。
实例1
以上冲击波求解方法的过程没有指定任何具体对象，也就是说对于工程中任何冲击波传播问题的求解均适用。为详细说明方法的操作流程以及应用效果，下面以经典Riemann问题为例详细说明本发明的应用。Riemann问题即初始间断的分解问题，是CFD的经典算例之一，它包含了CFD中的各种间断，这些间断问题的求解一直是CFD发展的难点和核心问题。Riemann问题自身存在解析解，可以用来校验数值模拟方法的正确性和精度，正因如此，几乎所有的流场求解方法都是建立在求解Riemann问题基础上的，而后向多维拓展。
Riemann问题描述的是一根封闭的管道内左右两侧有不同的气体，中间以一薄膜将两侧气体隔开，在计算开始的时刻，突然撤去中间的隔膜，两侧的气体将在压力的作用下互相混合，最终达到平衡状态,也称为激波管问题。本算例计算总时长为t＝0.2,算例中所有参数均为无量纲参数，因使用B样条映射函数时需要在边界处生成虚拟网格，本案例生成的背景网格总数为1002，在激波管左右两端外侧各设置1个虚网格，初始时刻除虚网格外每个网格体内包含2个物质点，因此物质点数量为2000，设置库朗数为0.8，气体比热比γ＝1.4，计算域为x∈(0,1)，初始时刻的间断分布为：
Figure BDA0002577215680000071
为了验证本发明方法的鲁棒性，分别采用精度较低的2次B样条基函数(3点插值)和精度较高的3次B样条基函数(4点插值)作为物质点与网格节点间的形函数，并将本发明方法的计算结果与传统MUSL格式、Riemann问题解析解作对比。
步骤1、改进传统MUSL格式中的动量修正步骤，基于当前时间步的形函数将下一时间步的物质点动量映射到网格节点上，得到修正后的下一时刻网格节点动量，取消传统MUSL格式中的动量修正步，基于下一时间步的物质点位置计算形函数，同时将下一时间步的物质点质量和动量映射到网格节点上得到修正后的网格节点质量和动量为：
Figure BDA0002577215680000072
Figure BDA0002577215680000081
其中
Figure BDA0002577215680000082
为网格节点i在下一计算步中的质量，/>
Figure BDA0002577215680000083
为网格节点i在下一计算步中的动量，mp为物质点质量，/>
Figure BDA0002577215680000084
为物质点p在下一计算步中的速度，/>
Figure BDA0002577215680000085
为下一计算步中网格节点i与物质点p之间的映射函数。
步骤2、建立激波管模型，将其视为长度为1的一维线模型，对激波管内流场进行离散化处理，在管内均匀地生成1000个网格，在管的左右端点外各设置1个虚网格，背景网格总数为1002，在除虚网格外的每个网格内设置2个物质点，物质点总数为2000。
步骤3、给定激波管的初始条件与边界条件，设置管内物质点在初始时刻的密度、速度、压力为：
Figure BDA0002577215680000086
计算出各个物质点的体积、质量、内能：
Figure BDA0002577215680000087
Figure BDA0002577215680000088
Figure BDA0002577215680000089
这里VΩ＝1，np＝2000。
对激波管的左右两端施加边界条件，令端点及其外侧的网格节点速度、动量、力始终为0。
步骤4、根据稳定性条件计算时间步长为：
Figure BDA00025772156800000810
Figure BDA00025772156800000811
其中Δtt为当前计算步的时间步长，CCFL＝0.8，
Figure BDA0002577215680000091
为物质点p当前时刻的声速，/>
Figure BDA0002577215680000092
为物质点p在当前计算步中沿x方向的速度，/>
Figure BDA0002577215680000093
为物质点p在当前计算步中的密度，/>
Figure BDA0002577215680000094
为物质点p在当前计算步中的压力。
步骤5、基于冲击波求解的改进MUSL格式物质点法求解冲击波流场，包括：
5.1：将物质点参数映射到网格节点上，得到网格节点的质量和动量为：
Figure BDA0002577215680000095
Figure BDA0002577215680000096
其中
Figure BDA0002577215680000097
为网格节点i在当前计算步中的质量，/>
Figure BDA0002577215680000098
为网格节点i在当前计算步中的动量，mp为物质点p的质量，/>
Figure BDA0002577215680000099
为物质点p在当前计算步中的速度，/>
Figure BDA00025772156800000910
为当前计算步中网格节点i与物质点p之间的映射函数。
5.2：计算网格节点力为：
Figure BDA00025772156800000911
其中fi t为网格节点i在当前计算步中的力值，
Figure BDA00025772156800000912
为当前计算步中网格节点i与物质点p之间的映射函数导数。
5.3：积分动量方程为：
Figure BDA00025772156800000913
其中
Figure BDA00025772156800000914
为网格节点i在当前计算步中的动量，/>
Figure BDA00025772156800000915
为网格节点i在下一个计算步中的动量。
5.4：计算得到下一时刻物质点的位置和速度为：
Figure BDA00025772156800000916
Figure BDA00025772156800000917
/>
其中
Figure BDA0002577215680000101
为物质点p在当前计算步中的速度，/>
Figure BDA0002577215680000102
为物质点p在下一个计算步中的速度，ni为网格节点数，/>
Figure BDA0002577215680000103
为物质点p在当前计算步中的位置，/>
Figure BDA0002577215680000104
为物质点p在下一个计算步中的位置。
5.5：取消传统MUSL格式中修正网格节点动量的步骤，再次将物质点的质量和动量映射到网格节点上，修正网格节点的质量和动量为：
Figure BDA0002577215680000105
Figure BDA0002577215680000106
5.6、更新物质点的密度为：
Figure BDA0002577215680000107
Figure BDA0002577215680000108
Figure BDA0002577215680000109
其中
Figure BDA00025772156800001010
为物质点p在下一个计算步中的应变增量，它是一个方阵，行列为j和k，
Figure BDA00025772156800001011
为物质点p在当前计算步中的密度，/>
Figure BDA00025772156800001012
为物质点p在下一个计算步中的密度。
5.7：利用理想气体状态方程更新物质点的压力为：
Figure BDA00025772156800001013
其中
Figure BDA00025772156800001014
为物质点p在下一个计算步中的内能。
5.8：判断计算时间是否已达到0.2，若否则转入步骤4，若是则结束进入下一步。
步骤6、对激波管内流场进行可视化处理，输出各个网格节点上的密度，得到流场的密度分布特征。
以2次B样条函数为形函数时，计算结果如图3所示，可以看出传统MUSL格式计算结果与解析解相差较大，密度曲线的色散、耗散都较大，且结果振荡严重；而本发明方法计算结果与解析解十分接近，且传统MUSL格式需要计算约1220步才能完成计算，本发明仅需1040步。以3次B样条函数为形函数时，计算结果如图4所示，随着映射函数精度的提高，本发明和传统MUSL格式的计算精度都有所提高，但本发明的色散和耗散更小，尤其是在间断处，本发明的密度曲线更接近于解析解，且收敛速度更快。
综上所述，在同一形函数下，本发明方法计算速度和计算精度要高于传统MUSL格式；本发明方法的鲁棒性较好，对形函数精度要求较低，在3点插值格式下也能取得较好的计算结果，而传统的MUSL格式则需在4点插值格式下才能较好地捕捉到冲击波特征，说明本发明方法的计算能耗更小。

Claims (6)

1.一种冲击波求解的改进MUSL格式物质点法，其特征在于，包括以下步骤：
步骤1、改进传统MUSL格式中的动量修正步骤：基于当前时间步的形函数将下一时间步的物质点动量映射到网格节点上，得到修正后的下一时刻网格节点动量，取消传统MUSL格式中的动量修正步，基于下一时间步的物质点位置计算形函数，同时将下一时间步的物质点质量和动量映射到网格节点上；
步骤2、建立冲击波流场模型，对冲击波流场进行离散化处理：划分背景网格并布设物质点；
步骤3、给定冲击波流场的初始条件与边界条件，设置计算总长；
步骤4、基于稳定性条件计算当前计算步的时间步长；设置的库朗数，计算出当前计算步的时间步长为：
Figure FDA0004006512500000011
Figure FDA0004006512500000012
其中Δtt为当前计算步的时间步长，CCFL为库朗数，
Figure FDA0004006512500000013
为物质点p当前时刻的声速，
Figure FDA0004006512500000014
分别为物质点p在当前计算步中沿x、y、z方向的速度，/>
Figure FDA0004006512500000015
为物质点p在当前计算步中的密度，/>
Figure FDA0004006512500000016
为物质点p在当前计算步中的压力，γ为气体比热比，L为网格长度，np为物质点总数；
步骤5、基于冲击波求解的改进MUSL格式物质点法求解冲击波流场；包括：
步骤5.1、将物质点参数映射到网格节点上，得到网格节点的质量和动量为：
Figure FDA0004006512500000017
Figure FDA0004006512500000018
其中
Figure FDA0004006512500000019
为网格节点i在当前计算步中的质量，/>
Figure FDA00040065125000000110
为网格节点i在当前计算步中的动量，
Figure FDA00040065125000000111
为物质点p在当前计算步中的速度，/>
Figure FDA00040065125000000112
为当前计算步中网格节点i与物质点p之间的映射函数，mp为物质点质量；
步骤5.2、计算网格节点力为：
Figure FDA0004006512500000021
其中fi t为网格节点i在当前计算步中的力值，
Figure FDA0004006512500000022
为当前计算步中网格节点i与物质点p之间的映射函数导数；
步骤5.3、积分动量方程为：
Figure FDA0004006512500000023
/>
其中
Figure FDA0004006512500000024
为网格节点i在当前计算步中的动量，/>
Figure FDA0004006512500000025
为网格节点i在下一个计算步中的动量；
步骤5.4、计算得到下一时刻物质点的位置和速度为：
Figure FDA0004006512500000026
Figure FDA0004006512500000027
其中
Figure FDA0004006512500000028
为物质点p在当前计算步中的速度，/>
Figure FDA0004006512500000029
为物质点p在下一个计算步中的速度，ni为网格节点数，/>
Figure FDA00040065125000000210
为物质点p在当前计算步中的位置，/>
Figure FDA00040065125000000211
为物质点p在下一个计算步中的位置；
步骤5.5、取消传统MUSL格式中修正网格节点动量的步骤，再次将物质点的质量和动量映射到网格节点上，修正网格节点的质量和动量为：
Figure FDA00040065125000000212
Figure FDA00040065125000000213
步骤5.6、更新物质点的密度为：
Figure FDA00040065125000000214
Figure FDA0004006512500000031
Figure FDA0004006512500000032
其中
Figure FDA0004006512500000033
为物质点p在下一个计算步中的应变增量，它是一个方阵，行列为j和k，/>
Figure FDA0004006512500000034
为物质点p在当前计算步中的密度，/>
Figure FDA0004006512500000035
为物质点p在下一个计算步中的密度；
步骤5.7、利用理想气体状态方程更新物质点的压力为：
Figure FDA0004006512500000036
其中
Figure FDA0004006512500000037
为物质点p在下一个计算步中的内能，γ为气体比热比；
步骤5.8、判断是否已计算至终点时刻，若否则令t＝t+1并转入步骤4，若是则结束进入下一步；
步骤6、对冲击波流场进行可视化处理：输出流场的密度、压力参数。
2.根据权利要求1所述的冲击波求解的改进MUSL格式物质点法，其特征在于，步骤1改进传统MUSL格式中的动量修正步骤，得到修正后的网格节点质量和动量为：
Figure FDA0004006512500000038
Figure FDA0004006512500000039
其中
Figure FDA00040065125000000310
为网格节点i在下一计算步中的质量，/>
Figure FDA00040065125000000311
为网格节点i在下一计算步中的动量，mp为物质点质量，/>
Figure FDA00040065125000000312
为物质点p在下一计算步中的速度，/>
Figure FDA00040065125000000313
为下一计算步中网格节点i与物质点p之间的映射函数。
3.根据权利要求1所述的冲击波求解的改进MUSL格式物质点法，其特征在于，步骤2建立冲击波流场模型，对冲击波流场进行离散化处理，对流场的边界进行几何清洗并建立拓扑，用均匀的结构网格覆盖整个计算域生成背景网格，在边界外侧至少生成一层网格作为虚网格，并在冲击波的计算范围内均匀地布设物质点。
4.根据权利要求1所述的冲击波求解的改进MUSL格式物质点法，其特征在于，步骤3给定冲击波流场的初始条件与边界条件，设置计算总长，给定初始时刻各个物质点的密度
Figure FDA0004006512500000041
压力/>
Figure FDA0004006512500000042
速度/>
Figure FDA0004006512500000043
计算出初始时刻各个物质点的体积/>
Figure FDA0004006512500000044
物质点质量mp、内能/>
Figure FDA0004006512500000045
Figure FDA0004006512500000046
Figure FDA0004006512500000047
Figure FDA0004006512500000048
其中VΩ是整个冲击波流场的体积，np为物质点总数，γ为气体比热比；
在背景网格上施加边界条件，对于固壁边界，边界处的网格和边界外的虚网格速度、动量、节点力始终为0。
5.根据权利要求1所述的冲击波求解的改进MUSL格式物质点法，其特征在于，步骤6对冲击波流场进行可视化处理，存储网格节点上的密度参数的时程数据，得到流场内冲击波的传播过程。
6.一种基于改进MUSL格式物质点法的冲击波求解***，其特征在于，包括物质点动量修正模块、冲击波流场模型建立模块、计算总长设置模块、时间步长计算模块、冲击波流场计算模块、可视化处理模块；
所述物质点动量修正模块用于基于当前时间步的形函数将下一时间步的物质点动量映射到网格节点上，得到修正后的下一时刻网格节点动量，取消传统MUSL格式中的动量修正步，基于下一时间步的物质点位置计算形函数，同时将下一时间步的物质点质量和动量映射到网格节点上；所述冲击波流场模型建立模块用于、建立冲击波流场模型，对冲击波流场进行离散化处理；所述计算总长设置模块用于给定冲击波流场的初始条件与边界条件，设置计算总长；所述时间步长计算模块用于计算当前计算步的时间步长；冲击波流场计算模块用于基于冲击波求解的改进MUSL格式物质点法求解冲击波流场；所述可视化处理模块用于对冲击波流场进行可视化处理：输出流场的密度、压力参数；
所述物质点动量修正模块处理得到修正后的网格节点质量和动量为：
Figure FDA0004006512500000049
Figure FDA0004006512500000051
其中
Figure FDA0004006512500000052
为网格节点i在下一计算步中的质量，/>
Figure FDA0004006512500000053
为网格节点i在下一计算步中的动量，mp为物质点质量，/>
Figure FDA0004006512500000054
为物质点p在下一计算步中的速度，/>
Figure FDA0004006512500000055
为下一计算步中网格节点i与物质点p之间的映射函数；
所述时间步长计算模块处理得到当前计算步的时间步长为：
Figure FDA0004006512500000056
Figure FDA0004006512500000057
其中Δtt为当前计算步的时间步长，CCFL为库朗数，
Figure FDA0004006512500000058
为物质点p当前时刻的声速，
Figure FDA0004006512500000059
分别为物质点p在当前计算步中沿x、y、z方向的速度，/>
Figure FDA00040065125000000510
为物质点p在当前计算步中的密度，/>
Figure FDA00040065125000000511
为物质点p在当前计算步中的压力，γ为气体比热比，L为网格长度，np为物质点总数。/>
CN202010657277.7A 2020-07-09 2020-07-09 冲击波求解的改进musl格式物质点法 Active CN111859645B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010657277.7A CN111859645B (zh) 2020-07-09 2020-07-09 冲击波求解的改进musl格式物质点法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010657277.7A CN111859645B (zh) 2020-07-09 2020-07-09 冲击波求解的改进musl格式物质点法

Publications (2)

Publication Number Publication Date
CN111859645A CN111859645A (zh) 2020-10-30
CN111859645B true CN111859645B (zh) 2023-03-31

Family

ID=73151958

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010657277.7A Active CN111859645B (zh) 2020-07-09 2020-07-09 冲击波求解的改进musl格式物质点法

Country Status (1)

Country Link
CN (1) CN111859645B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112364362B (zh) * 2020-11-16 2023-12-29 宁波九寰适创科技有限公司 一种面向流体仿真方向的并行多层自适应局部加密方法
CN113673185B (zh) * 2021-08-25 2024-01-30 北京理工大学 一种加权双向映射的精确捕捉激波间断面方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107766287A (zh) * 2017-10-26 2018-03-06 哈尔滨工程大学 一种应用于***冲击工程中的基于物质点法的随机动力学分析方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8666715B2 (en) * 2009-03-31 2014-03-04 Airbus Operations S.L. Method and system for a quick calculation of aerodynamic forces on an aircraft in transonic conditions

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107766287A (zh) * 2017-10-26 2018-03-06 哈尔滨工程大学 一种应用于***冲击工程中的基于物质点法的随机动力学分析方法

Also Published As

Publication number Publication date
CN111859645A (zh) 2020-10-30

Similar Documents

Publication Publication Date Title
Farhat et al. Robust and provably second‐order explicit–explicit and implicit–explicit staggered time‐integrators for highly non‐linear compressible fluid–structure interaction problems
Wang Adaptive high-order methods in computational fluid dynamics
CN111859645B (zh) 冲击波求解的改进musl格式物质点法
CN108446445B (zh) 一种基于气动力降阶模型的复合材料机翼优化设计方法
Eymann et al. Cartesian adaptive mesh refinement with the HPCMP CREATE™-AV kestrel solver
CN108563843B (zh) 定常可压缩流动的扰动区域更新方法
US8775140B2 (en) Time and space scaled S-model for turbulent fluid flow simulations
CN113673027B (zh) 一种基于代理模型的高超声速飞行器气动载荷优化设计方法
CN110502788B (zh) 一种考虑密封条非线性压缩特性的车门变形获取方法
CN116738891B (zh) 一种增强飞行器流场模拟稳定性的lu-sgs改进方法
CN115496006A (zh) 一种适用于高超声速飞行器的高精度数值模拟方法
CN114492250A (zh) 基于递归分解的曲面网格生成方法及***、计算机设备
Hou et al. Discrete shape sensitivity equations for aerodynamic problems
Elmiligui et al. USM3D Simulations for Third Sonic Boom Workshop
CN110321588B (zh) 基于数值模拟的轨道车辆空气阻力计算方法
CN105718619A (zh) 一种基于有限元法的飞行器燃油质量特性确定方法
Darmofal et al. Eigenmode analysis of boundary conditions for the one-dimensional preconditioned Euler equations
CN113378440A (zh) 一种流固耦合数值模拟计算方法、装置及设备
CN116341421B (zh) 高超声速流场数值模拟方法、***、电子设备及存储介质
CN111859646B (zh) 一种基于b样条映射函数物质点法的冲击波变步长求解方法
Da Ronch et al. Modeling of unsteady aerodynamic loads
CN110909511B (zh) 一种无曲面体积分的无粘低速绕流数值模拟方法
Murayama et al. Summary of JAXA studies for the fifth AIAA CFD drag prediction workshop using UPACS and FaSTAR
Righi et al. Uncertainties Quantification of CFD-Based Flutter Prediction
CN104778321B (zh) 一种复杂形状声学空间的分析预报方法

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