CN111859646A - 一种基于b样条映射函数物质点法的冲击波变步长求解方法 - Google Patents

一种基于b样条映射函数物质点法的冲击波变步长求解方法 Download PDF

Info

Publication number
CN111859646A
CN111859646A CN202010658129.7A CN202010658129A CN111859646A CN 111859646 A CN111859646 A CN 111859646A CN 202010658129 A CN202010658129 A CN 202010658129A CN 111859646 A CN111859646 A CN 111859646A
Authority
CN
China
Prior art keywords
shock wave
material point
flow field
grid
solving
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.)
Granted
Application number
CN202010658129.7A
Other languages
English (en)
Other versions
CN111859646B (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 CN202010658129.7A priority Critical patent/CN111859646B/zh
Publication of CN111859646A publication Critical patent/CN111859646A/zh
Application granted granted Critical
Publication of CN111859646B publication Critical patent/CN111859646B/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
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/12—Timing analysis or timing optimisation
    • 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

本发明提出了一种基于B样条映射函数物质点法的冲击波变步长求解方法，首先建立B样条映射函数；然后建立冲击波流场模型，对冲击波流场进行离散化处理；再给定冲击波流场的初始条件与边界条件，设置计算总长；接着根据稳定性条件计算当前计算步的时间步长；利用物质点法求解冲击波问题；最后对冲击波流场进行可视化处理。本发明提出的方法能避免欧拉法中对流项处理困难的问题和拉格朗日法中的网格畸变问题，可有效地减少物质点穿越网格时的噪声，相比于传统的标准物质点法和GIMP法，本发明能更好地抑制振荡，稳定性和鲁棒性更强，计算效率更高，可为工程中冲击波的求解提供一种新方法。

Description

一种基于B样条映射函数物质点法的冲击波变步长求解方法
技术领域
本发明属于冲击波流场求解方法技术领域，特别是一种基于B样条映射函数物质点法的冲击波变步长求解方法。
背景技术
冲击波是一种不连续的峰在流体介质中的传播，导致流场的压力、温度、密度等物理性质发生阶跃式变化。其在工程中广泛地存在，例如高速列车运行时形成列车冲击波、矿山爆破、防爆车辆等，可能造成人员的伤亡和设备的损坏。随着计算机技术的发展，CFD技术逐渐被应用于冲击波的仿真模拟，通过数值计算冲击波的传播过程，得到其对流场的扰动作用，进而可为其作用对象的气动载荷计算提供计算参数。此外，冲击波问题是CFD的经典算例之一，波前和波后流场的密度、压力等参数都呈现间断性，这些间断问题的求解一直是CFD发展的难点和核心问题。因此，探究一种新型的冲击波求解方法具有广泛的工程价值。
目前，广泛地采用有限体积法求解冲击波流场，该方法严重地依赖网格，网格划分占整体工作量的比重过大，当计算物体的形状较为复杂时会造成网格扭曲、负体积网格等，造成计算出现奇异解进而发散。且在求解波面两侧的间断时，波面内的网格数量对计算精度的影响较大，若网格过少，则会产生耗散，若增加网格数量，一方面冲击波的传播轨迹是未知的，无法在前处理阶段确定网格加密区，另一方面，网格数量的增加势必造成计算耗时的延长。
物质点法是一种混合了拉格朗日法和欧拉法的数值计算方法，已经被广泛地应用于固体结构中的高速冲击问题和爆炸问题求解。该方法用物质点来离散物体，相关物理参数均由物质点携带，背景网格仅用于积分动量方程和求解梯度，因而只需简单地生成覆盖计算域的结构网格即可。物质点法在一定程度上摆脱了网格的约束，避免了大变形造成的网格畸变问题，能很好地处理材料断裂失效、结构破碎等现象，并且不会造成质量损失，适用于求解极端工况。然而，在利用物质点法求解冲击波问题时，由于传统的物质点法简单地采用线性函数来描述物质点和网格节点之间的关系，该映射函数下物质点的影响域小，函数导数不连续，当物质点穿越网格的时候，会产生较为严重的振荡现象。冲击波波面两侧的物理参数呈现阶跃式变化，若数值模拟方法产生的振荡过大，很容易出现负压力、负密度等非物理解。且仅在计算的初始阶段基于初始条件计算稳定步长，对于冲击波问题来说，流场各处的当地声速受冲击波的影响会发生较大的变化，依照初始条件计算得到的步长在计算过程中很可能不再满足稳定性条件，进而导致计算发散。
发明内容
本发明的目的在于提供一种基于B样条映射函数物质点法的冲击波变步长求解方法，为实际工程中预测冲击波的影响提供方法。
实现本发明目的的技术解决方案为：
一种基于B样条映射函数物质点法的冲击波变步长求解方法，包括以下步骤：
步骤1、建立B样条映射函数：推导出B样条插值函数及其导数；
步骤2、建立冲击波流场模型，对冲击波流场进行离散化处理：划分背景网格并布设物质点；
步骤3、给定冲击波流场的初始条件与边界条件，设置计算总长；
步骤4、计算当前计算步的时间步长；
步骤5、利用物质点法求解冲击波问题；
步骤6、对冲击波流场进行可视化处理：输出流场的密度、压力参数。
一种基于B样条映射函数物质点法的冲击波变步长求解系统，包括B样条映射函数建立模块、离散化处理模块、计算总长设置模块、时间步长计算模块、冲击波问题求解模块、可视化处理模块；
所述B样条映射函数建立模块用于建立B样条映射函数，推导出B样条插值函数及其导数；所述离散化处理模块用于建立冲击波流场模型，对冲击波流场进行离散化处理；所述计算总长设置模块用于给定冲击波流场的初始条件与边界条件，设置计算总长；所述时间步长计算模块用于计算当前计算步的时间步长；所述冲击波问题求解模块用于利用物质点法求解冲击波问题；所述可视化处理模块用于对冲击波流场进行可视化处理，输出流场的密度、压力参数。
本发明与现有技术相比，其显著优点是：
(1)利用物质点法求解冲击波问题，生成覆盖计算域的规则结构网格，大大减少了网格划分的工作量和时间；单个物质点的质量在计算中始终不变，自动满足质量守恒定律，无需求解质量方程；压力、密度等物理参数均储存在物质点上，降低其数值耗散；在单个计算步中采用拉格朗日法求解，避免了欧拉法中对流项处理困难的问题，在两个计算步之间抛弃变形的背景网格，重新建立物质点和网格之间的映射关系，避免了拉格朗日法中的网格畸变问题。
(2)利用B样条基函数推导得到B样条映射函数，在B样条映射函数下，物质点的影响半径为两个网格长度，即每个物质点都和4个网格节点存在映射关系，高于标准物质点法中线性映射函数下的2个存在映射关系的网格节点和GIMP法中3个存在映射关系的网格节点。有效减少了物质点穿越背景网格时产生的噪声，增加了单个物质点的影响域，可更精确地捕捉冲击波波面前后流场的间断特征，并降低波面处计算结果的色散和耗散，有效地抑制冲击波流场振荡，提高了数值模拟的计算精度。
(3)采用变步长计算方法，在每一个计算步的开始时基于当前流场的参数值计算时间步长，避免了定步长法时计算过程中冲击波流场变化剧烈导致的计算发散，当流场趋于稳定时，步长会随之增大，更有利于冲击波流场细节的捕捉，并提高了收敛速度。
附图说明
图1为本发明冲击波求解方法的流程图。
图2为Riemann问题示意图。
图3为本发明和GIMP法对Riemann问题的求解结果与其解析解的对比图。
具体实施方式
下面结合附图及具体实施例对本发明做进一步的介绍。
步骤1、建立B样条插值函数及其导数为：
Figure BDA0002577510030000031
Figure BDA0002577510030000041
Figure BDA0002577510030000042
其中Ni,p为物质点p和网格节点i之间的映射函数，xp为物质点的全局坐标，xi为网格节点的全局坐标，L为网格长度。在B样条映射函数下，物质点p的影响半径为2L，即每个物质点都和4个网格节点存在映射关系，高于标准物质点法中线性映射函数下的2个存在映射关系的网格节点和GIMP法中3个存在映射关系的网格节点。
步骤2、建立冲击波流场模型，对冲击波流场进行离散化处理，对流场的边界进行几何清洗并建立拓扑，用均匀的结构网格覆盖整个计算域生成背景网格，在边界外侧至少生成一层网格作为虚网格，并在冲击波的计算范围内均匀地布设物质点。
步骤3、给定冲击波流场的初始条件与边界条件，设置计算总长，给定初始时刻各个物质点的密度
Figure BDA0002577510030000043
压力
Figure BDA0002577510030000044
速度
Figure BDA0002577510030000045
计算出初始时刻各个物质点的体积
Figure BDA0002577510030000046
质量mp、内能
Figure BDA0002577510030000047
Figure BDA0002577510030000048
Figure BDA0002577510030000049
Figure BDA00025775100300000410
其中VΩ是整个冲击波流场的体积，np为物质点总数，γ为气体比热比，空气取1.4。
在背景网格上施加边界条件，对于固壁边界来说，边界处的网格和边界外的虚网格速度、动量、节点力始终为0。
步骤4、根据稳定性条件计算时间步长为：
Figure BDA0002577510030000051
Figure BDA0002577510030000052
其中Δtt为当前计算步的时间步长，CCFL为库朗数，
Figure BDA0002577510030000053
为物质点p当前时刻的声速，
Figure BDA0002577510030000054
分别为物质点p在当前计算步中沿x、y、z方向的速度，
Figure BDA0002577510030000055
为物质点p在当前计算步中的密度，
Figure BDA0002577510030000056
为物质点p在当前计算步中的压力。
步骤5、利用物质点法求解冲击波问题，包括：
5.1：将物质点参数映射到网格节点上，得到网格节点的质量和动量为：
Figure BDA0002577510030000057
Figure BDA0002577510030000058
其中
Figure BDA0002577510030000059
为网格节点i在当前计算步中的质量，
Figure BDA00025775100300000510
为网格节点i在当前计算步中的动量，
Figure BDA00025775100300000511
为物质点p在当前计算步中的速度，
Figure BDA00025775100300000512
为当前计算步中网格节点i与物质点p之间的映射函数。
5.2：计算网格节点力为：
Figure BDA00025775100300000513
其中fi t为网格节点i在当前计算步中的力值，
Figure BDA00025775100300000514
为当前计算步中网格节点i与物质点p之间的映射函数导数。
5.3：积分动量方程为：
Figure BDA00025775100300000515
其中
Figure BDA00025775100300000516
为网格节点i在当前计算步中的动量，
Figure BDA00025775100300000517
为网格节点i在下一个计算步中的动量。
5.4：计算得到下一时刻物质点的位置和速度为：
Figure BDA0002577510030000061
Figure BDA0002577510030000062
其中
Figure BDA0002577510030000063
为物质点p在当前计算步中的速度，
Figure BDA0002577510030000064
为物质点p在下一个计算步中的速度，ni为网格节点数，
Figure BDA0002577510030000065
为物质点p在当前计算步中的位置，
Figure BDA0002577510030000066
为物质点p在下一个计算步中的位置。
5.5：修正网格节点的动量为：
Figure BDA0002577510030000067
5.6：更新物质点的密度为：
Figure BDA0002577510030000068
Figure BDA0002577510030000069
Figure BDA00025775100300000610
其中
Figure BDA00025775100300000611
为物质点p在下一个计算步中的应变增量，它是一个方阵，行列为j和k，
Figure BDA00025775100300000612
为物质点p在当前计算步中的密度，
Figure BDA00025775100300000613
为物质点p在下一个计算步中的密度。
5.7：利用理想气体状态方程更新物质点的压力为：
Figure BDA00025775100300000614
其中
Figure BDA00025775100300000615
为物质点p在下一个计算步中的内能。
5.8：判断是否已计算至终点时刻，若否则转入4，若是则结束。
步骤6、对冲击波流场进行可视化处理，存储网格节点上的密度等参数的时程数据，得到流场内冲击波的传播过程。
基于上述的基于B样条映射函数物质点法的冲击波变步长求解方法和计算机程序，本发明还提出了一种基于B样条映射函数物质点法的冲击波变步长求解系统，具体包括B样条映射函数建立模块、离散化处理模块、计算总长设置模块、时间步长计算模块、冲击波问题求解模块、可视化处理模块；
所述B样条映射函数建立模块用于建立B样条映射函数，推导出B样条插值函数及其导数，具体过程见上述步骤1；所述离散化处理模块用于建立冲击波流场模型，对冲击波流场进行离散化处理，具体过程见上述步骤2；所述计算总长设置模块用于给定冲击波流场的初始条件与边界条件，设置计算总长，具体过程见上述步骤3；所述时间步长计算模块用于计算当前计算步的时间步长，具体过程见上述步骤4；所述冲击波问题求解模块用于利用物质点法求解冲击波问题，具体过程见上述步骤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 BDA0002577510030000071
步骤1、建立B样条插值函数及其导数为：
Figure BDA0002577510030000081
Figure BDA0002577510030000082
Figure BDA0002577510030000083
其中Ni,p为物质点p和网格节点i之间的映射函数，xp为物质点的全局坐标，xi为网格节点的全局坐标，L为网格长度。
步骤2、建立激波管模型，将其视为长度为1的一维线模型，对激波管内流场进行离散化处理，在管内均匀地生成1000个网格，在管的左右端点外各设置1个虚网格，背景网格总数为1002，在除虚网格外的每个网格内设置2个物质点，物质点总数为2000。
步骤3、给定激波管的初始条件与边界条件，设置管内物质点在初始时刻的密度、速度、压力为：
Figure BDA0002577510030000084
计算出各个物质点的体积、质量、内能：
Figure BDA0002577510030000085
Figure BDA0002577510030000091
Figure BDA0002577510030000092
这里VΩ＝1，np＝2000。
对激波管的左右两端施加边界条件，令端点及其外侧的网格节点速度、动量、力始终为0。
步骤4、根据稳定性条件计算时间步长为：
Figure BDA0002577510030000093
Figure BDA0002577510030000094
其中Δtt为当前计算步的时间步长，CCFL＝0.8，
Figure BDA0002577510030000095
为物质点p当前时刻的声速，
Figure BDA0002577510030000096
为物质点p在当前计算步中沿x方向的速度，
Figure BDA0002577510030000097
为物质点p在当前计算步中的密度，
Figure BDA0002577510030000098
为物质点p在当前计算步中的压力。
步骤5、利用物质点法求解Riemann问题，包括：
5.1：将物质点参数映射到网格节点上，得到网格节点的质量和动量为：
Figure BDA0002577510030000099
Figure BDA00025775100300000910
其中
Figure BDA00025775100300000911
为网格节点i在当前计算步中的质量，
Figure BDA00025775100300000912
为网格节点i在当前计算步中的动量，mp为物质点p的质量，
Figure BDA00025775100300000913
为物质点p在当前计算步中的速度，
Figure BDA00025775100300000914
为当前计算步中网格节点i与物质点p之间的映射函数。
5.2：计算网格节点力为：
Figure BDA00025775100300000915
其中fi t为网格节点i在当前计算步中的力值，
Figure BDA0002577510030000101
为当前计算步中网格节点i与物质点p之间的映射函数导数。
5.3：积分动量方程为：
Figure BDA0002577510030000102
其中
Figure BDA0002577510030000103
为网格节点i在当前计算步中的动量，
Figure BDA0002577510030000104
为网格节点i在下一个计算步中的动量。
5.4：计算得到下一时刻物质点的位置和速度为：
Figure BDA0002577510030000105
Figure BDA0002577510030000106
其中
Figure BDA0002577510030000107
为物质点p在当前计算步中的速度，
Figure BDA0002577510030000108
为物质点p在下一个计算步中的速度，ni为网格节点数，
Figure BDA0002577510030000109
为物质点p在当前计算步中的位置，
Figure BDA00025775100300001010
为物质点p在下一个计算步中的位置。
5.5：修正网格节点的动量为：
Figure BDA00025775100300001011
5.6：更新物质点的密度为：
Figure BDA00025775100300001012
Figure BDA00025775100300001013
Figure BDA00025775100300001014
其中
Figure BDA00025775100300001015
为物质点p在下一个计算步中的应变增量，这里它是一个数，
Figure BDA00025775100300001016
为物质点p在当前计算步中的密度，
Figure BDA00025775100300001017
为物质点p在下一个计算步中的密度。
5.7：利用理想气体状态方程更新物质点的压力为：
Figure BDA0002577510030000111
其中
Figure BDA0002577510030000112
为物质点p在下一个计算步中的内能。
5.8：判断计算时间是否已达到0.2，若否则转入步骤4，若是则结束进入下一步。
步骤6、对激波管内流场进行可视化处理，输出各个网格节点上的密度，得到流场的密度分布特征。
利用标准物质点法计算Riemann问题至t＝0.03时出现非物理解，无法继续求解；本发明方法和GIMP法计算得到的Riemann问题的结果如图3所示，可以看出相比于GIMP法，本发明的计算结果与解析解更接近，得到的密度曲线振荡较小，色散、耗散都低于GIMP法。此外，对于实施例1，GIMP法需要计算约1090个计算步，本发明方法仅需计算约1020个计算步，说明本发明方法的计算效率更高。
本发明的基于B样条映射函数物质点法的冲击波变步长求解方法，应用范围广，不仅可以应用于实施例1中的经典Riemann问题，还可以应用于其他冲击波问题的求解，上述实施例和说明书中描述的只是说明本发明的原理及技术效果，在不脱离本发明的精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围有所附权利要求书及其等效物界定。

Claims (10)

1.一种基于B样条映射函数物质点法的冲击波变步长求解方法，其特征在于，包括以下步骤：
步骤1、建立B样条映射函数：推导出B样条插值函数及其导数；
步骤2、建立冲击波流场模型，对冲击波流场进行离散化处理：划分背景网格并布设物质点；
步骤3、给定冲击波流场的初始条件与边界条件，设置计算总长；
步骤4、计算当前计算步的时间步长；
步骤5、利用物质点法求解冲击波问题；
步骤6、对冲击波流场进行可视化处理：输出流场的密度、压力参数。
2.根据权利要求1所述的一种基于B样条插值物质点法的冲击波变步长求解方法，其特征在于，步骤1建立B样条插值函数及其导数分别为：
Figure FDA0002577510020000011
Figure FDA0002577510020000012
其中Ni,p为物质点p和网格节点i之间的映射函数，
Figure FDA0002577510020000013
xp为物质点的全局坐标，xi为网格节点的全局坐标，L为网格长度。
3.根据权利要求1所述的一种基于B样条插值物质点法的冲击波变步长求解方法，其特征在于，步骤2建立冲击波流场模型，对冲击波流场进行离散化处理，对流场的边界进行几何清洗并建立拓扑，用均匀的结构网格覆盖整个计算域生成背景网格，在边界外侧至少生成一层网格作为虚网格，并在冲击波的计算范围内均匀地布设物质点。
4.根据权利要求1所述的一种基于B样条插值物质点法的冲击波变步长求解方法，其特征在于，步骤3给定冲击波流场的初始条件与边界条件，设置计算总长，给定初始时刻各个物质点的密度
Figure FDA0002577510020000021
压力
Figure FDA0002577510020000022
速度
Figure FDA0002577510020000023
计算出初始时刻各个物质点的体积
Figure FDA0002577510020000024
质量mp、内能
Figure FDA0002577510020000025
Figure FDA0002577510020000026
Figure FDA0002577510020000027
Figure FDA0002577510020000028
其中VΩ是整个冲击波流场的体积，np为物质点总数，γ为气体比热比；
在背景网格上施加边界条件，对于固壁边界来说，边界处的网格和边界外的虚网格速度、动量、节点力始终为0；设置冲击波计算的总时长。
5.根据权利要求1所述的一种基于B样条插值物质点法的冲击波变步长求解方法，其特征在于，步骤4计算当前计算步的时间步长，基于稳定性条件和设置的库朗数，计算出当前计算步的时间步长为：
Figure FDA0002577510020000029
Figure FDA00025775100200000210
其中Δtt为当前计算步的时间步长，CCFL为库朗数，
Figure FDA00025775100200000211
为物质点p当前时刻的声速，
Figure FDA00025775100200000212
分别为物质点p在当前计算步中沿x、y、z方向的速度，
Figure FDA00025775100200000213
为物质点p在当前计算步中的密度，
Figure FDA00025775100200000214
为物质点p在当前计算步中的压力，γ为气体比热比，np为物质点总数。
6.根据权利要求1所述的一种基于B样条插值物质点法的冲击波变步长求解方法，其特征在于，步骤5利用物质点法求解冲击波问题，包括：
步骤5.1、将物质点参数映射到网格节点上，得到网格节点的质量和动量为：
Figure FDA0002577510020000031
Figure FDA0002577510020000032
其中
Figure FDA0002577510020000033
为网格节点i在当前计算步中的质量，
Figure FDA0002577510020000034
为网格节点i在当前计算步中的动量，
Figure FDA0002577510020000035
为物质点p在当前计算步中的速度，
Figure FDA0002577510020000036
为当前计算步中网格节点i与物质点p之间的映射函数；
步骤5.2：计算网格节点力为：
Figure FDA0002577510020000037
其中fi t为网格节点i在当前计算步中的力值，
Figure FDA0002577510020000038
为当前计算步中网格节点i与物质点p之间的映射函数导数；
步骤5.3、积分动量方程为：
Figure FDA0002577510020000039
其中
Figure FDA00025775100200000310
为网格节点i在当前计算步中的动量，
Figure FDA00025775100200000311
为网格节点i在下一个计算步中的动量。
步骤5.4、计算得到下一时刻物质点的位置和速度为：
Figure FDA00025775100200000312
Figure FDA00025775100200000313
其中
Figure FDA00025775100200000314
为物质点p在当前计算步中的速度，
Figure FDA00025775100200000315
为物质点p在下一个计算步中的速度，ni为网格节点数，
Figure FDA00025775100200000316
为物质点p在当前计算步中的位置，
Figure FDA00025775100200000317
为物质点p在下一个计算步中的位置；
步骤5.5：修正网格节点的动量为：
Figure FDA0002577510020000041
步骤5.6：更新物质点的密度为：
Figure FDA0002577510020000042
Figure FDA0002577510020000043
Figure FDA0002577510020000044
其中
Figure FDA0002577510020000045
为物质点p在下一个计算步中的应变增量，它是一个方阵，行列为j和k，
Figure FDA0002577510020000046
为物质点p在当前计算步中的密度，
Figure FDA0002577510020000047
为物质点p在下一个计算步中的密度；
步骤5.7：利用理想气体状态方程更新物质点的压力为：
Figure FDA0002577510020000048
其中
Figure FDA0002577510020000049
为物质点p在下一个计算步中的内能；
步骤5.8：判断是否已计算至终点时刻，若否则转入步骤4，若是则进入下一步。
7.根据权利要求1中所述的一种基于B样条插值物质点法的冲击波变步长求解方法，其特征在于，步骤6对冲击波流场进行可视化处理，存储网格节点上的密度等参数的时程数据，得到流场内冲击波的传播过程。
8.一种基于B样条映射函数物质点法的冲击波变步长求解系统，其特征在于，包括B样条映射函数建立模块、离散化处理模块、计算总长设置模块、时间步长计算模块、冲击波问题求解模块、可视化处理模块；
所述B样条映射函数建立模块用于建立B样条映射函数，推导出B样条插值函数及其导数；所述离散化处理模块用于建立冲击波流场模型，对冲击波流场进行离散化处理；所述计算总长设置模块用于给定冲击波流场的初始条件与边界条件，设置计算总长；所述时间步长计算模块用于计算当前计算步的时间步长；所述冲击波问题求解模块用于利用物质点法求解冲击波问题；所述可视化处理模块用于对冲击波流场进行可视化处理，输出流场的密度、压力参数。
9.根据权利要求8所述的冲击波变步长求解系统，其特征在于，所述冲击波流场进行离散化处理具体处理过程为：对冲击波流场进行离散化处理，对流场的边界进行几何清洗并建立拓扑，用均匀的结构网格覆盖整个计算域生成背景网格，在边界外侧至少生成一层网格作为虚网格，并在冲击波的计算范围内均匀地布设物质点。
10.根据权利要求8所述的冲击波变步长求解系统，其特征在于，所述时间步长计算模块计算出当前计算步的时间步长为：
Figure FDA0002577510020000051
Figure FDA0002577510020000052
其中Δtt为当前计算步的时间步长，CCFL为库朗数，
Figure FDA0002577510020000053
为物质点p当前时刻的声速，
Figure FDA0002577510020000054
Figure FDA0002577510020000055
分别为物质点p在当前计算步中沿x、y、z方向的速度，
Figure FDA0002577510020000056
为物质点p在当前计算步中的密度，
Figure FDA0002577510020000057
为物质点p在当前计算步中的压力，γ为气体比热比，np为物质点总数。
CN202010658129.7A 2020-07-09 2020-07-09 一种基于b样条映射函数物质点法的冲击波变步长求解方法 Active CN111859646B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010658129.7A CN111859646B (zh) 2020-07-09 2020-07-09 一种基于b样条映射函数物质点法的冲击波变步长求解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010658129.7A CN111859646B (zh) 2020-07-09 2020-07-09 一种基于b样条映射函数物质点法的冲击波变步长求解方法

Publications (2)

Publication Number Publication Date
CN111859646A true CN111859646A (zh) 2020-10-30
CN111859646B CN111859646B (zh) 2023-06-09

Family

ID=73152091

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010658129.7A Active CN111859646B (zh) 2020-07-09 2020-07-09 一种基于b样条映射函数物质点法的冲击波变步长求解方法

Country Status (1)

Country Link
CN (1) CN111859646B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115329646A (zh) * 2022-10-12 2022-11-11 国家超级计算天津中心 冲击波仿真方法、装置、设备及介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100250205A1 (en) * 2009-03-31 2010-09-30 Airbus Espana, S.L. Method and system for a quick calculation of aerodynamic forces on an aircraft in transonic conditions
CN107766287A (zh) * 2017-10-26 2018-03-06 哈尔滨工程大学 一种应用于爆炸冲击工程中的基于物质点法的随机动力学分析方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100250205A1 (en) * 2009-03-31 2010-09-30 Airbus Espana, S.L. Method and system for a quick calculation of aerodynamic forces on an aircraft in transonic conditions
CN107766287A (zh) * 2017-10-26 2018-03-06 哈尔滨工程大学 一种应用于爆炸冲击工程中的基于物质点法的随机动力学分析方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115329646A (zh) * 2022-10-12 2022-11-11 国家超级计算天津中心 冲击波仿真方法、装置、设备及介质
CN115329646B (zh) * 2022-10-12 2023-03-10 国家超级计算天津中心 冲击波仿真方法、装置、设备及介质

Also Published As

Publication number Publication date
CN111859646B (zh) 2023-06-09

Similar Documents

Publication Publication Date Title
US10296672B2 (en) Generating inviscid and viscous fluid-flow simulations over an aircraft surface using a fluid-flow mesh
CN111859645B (zh) 冲击波求解的改进musl格式物质点法
US9348956B2 (en) Generating a simulated fluid flow over a surface using anisotropic diffusion
CN110795869B (zh) 流场数据的数值计算方法和装置
Parsani et al. High-order accurate entropy-stable discontinuous collocated Galerkin methods with the summation-by-parts property for compressible CFD frameworks: Scalable SSDC algorithms and flow solver
CN112069689A (zh) 一种航空发动机燃油雾化特性的仿真方法及系统
CN111859646A (zh) 一种基于b样条映射函数物质点法的冲击波变步长求解方法
Zangeneh Development of a new algorithm for modeling viscous transonic flow on unstructured grids at high Reynolds numbers
Borisov et al. Parallel implementation of an implicit scheme based on the LU-SGS method for 3D turbulent flows
CN113378440A (zh) 一种流固耦合数值模拟计算方法、装置及设备
Coirier et al. A Cartesian, cell-based approach for adaptively-refined solutions of the Euler and Navier-Stokes equations
Zhang et al. A dual interpolation boundary face method for 3D elasticity
Walsh et al. A preconditioned lattice Boltzmann flux solver for steady flows on unstructured hexahedral grids
Ji Comparison and analysis of different numerical schemes in sod’s one-dimensional shock tube problems
Chen et al. Subdivision Surface-Based Isogeometric Boundary Element Method for Steady Heat Conduction Problems with Variable Coefficient.
CN110909511B (zh) 一种无曲面体积分的无粘低速绕流数值模拟方法
Isaev et al. Testing of numerical methods, convective schemes, algorithms for approximation of flows, and grid structures by the example of a supersonic flow in a step-shaped channel with the use of the CFX and fluent packages
Kim et al. Three-dimensional building-cube method for inviscid compressible flow computations
Jiang et al. A study of three-dimensional acoustic scattering by arbitrary distribution multibodies using extended immersed boundary method
Assonitis et al. A shock-fitting technique for 2D/3D flows with interactions using structured grids
Tong et al. High-order methods for three-dimensional strand-cartesian grids
Bigarella et al. A study of convective flux computation schemes for aerodynamic flows
Yoshimura Parallel partitioned simulations of real world’s coupled problems
Korneev et al. Effective solving of three-dimensional gas dynamics problems with the Runge-Kutta discontinuous Galerkin method
Ban et al. Fundamental Design Optimization of an Innovative Supersonic Transport Configuration and Its Design Knowledge Extraction

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