1.本发明属于计算流体力学技术领域，尤其涉及是一种多介质界面两侧物理量一阶空间导数预测值计算方法，所获得的物理量一阶空间导数预测值尤其适用于对一般状态方程下的带对称性源项多介质界面高阶处理方法中。
背景技术：2.在水下爆炸模拟等多介质耦合问题的数值模拟中，保持界面处速度压强平衡的物理特性是非常重要的技术难点，它表现为若在数值格式的设计过程中不满足界面处的压强与速度平衡关系，则对界面处的数值模拟结果会出现压强错位的非物理现象。
3.现有技术已经说明，通过建立和求解界面处具有线性初值分布的带对称性源项多介质广义riemann问题来预测界面状态和空间导数，从而对虚拟区域进行线性赋值的方法可以有效地消除一阶误差。然而实际上，精确求解多介质广义riemann问题存在很多困难。首先，精确求解该问题依赖于状态方程，对于较为复杂的状态方程，能否精确推导一阶空间导数的显示形式尚未可知。其次，即使存在精确解，由于其解方案的复杂性，应用于大规模计算时成本较高，不利于工程上的应用。
4.中国发明专利申请cn114254573a公开了一种带对称性源项多介质界面高阶处理方法，其采用了获取多介质界面处两侧的物理量的状态值及一阶空间导数，建立界面处的多介质riemann问题以及具有线性初值分布的带对称性源项多介质广义riemann问题，并分别获得界面两侧物理量的预测值和一阶空间导数的预测值，在界面两侧建立两个虚拟区域分别作为两种介质的边界区域，通过使用获得的预测值对两个虚拟区域进行线性赋值，将原多介质问题解耦为两个单介质问题，采用时间离散方法在时间上推进，分别计算。
5.其中，在其获得界面两侧物理量一阶空间导数的预测值方法上，其采用了针对特殊状态方程下广义黎曼不变量的高阶近似以及密度空间导数的高阶近似的预测方案，获得了较为真实的仿真结果。但其仍存在一定不足，即黎曼不变量高阶近似的表达形式较为复杂，求解方案局限于一般状态方程，不能广泛地应用于包含一般形势下状态方程的仿真模拟。
技术实现要素：6.针对以上问题，本发明旨在提出一种界面两侧物理量一阶空间导数预测值的获得方法，所获得的物理量一阶空间导数预测值尤其适用于对一般状态方程下的带对称性源项多介质界面高阶处理方法中。
7.本发明的具体技术方案如下：
8.一种多介质界面两侧物理量一阶空间导数预测值计算方法，通过求解具有线性初值分布的一般状态方程下的带对称性源项多介质广义riemann问题，获得界面两侧物理量一阶空间导数的预测值；包括以下步骤：
9.1)建立关于介质一的关系式：
[0010][0011]
其中，ua为界面处速度随时间的变化率，pa为界面处压强随时间的变化率，ρ
1*
为介质一在界面处的密度的预测值，c
1*
为介质一在界面处的声速的预测值，z
l
为仅和介质一初值相关的参数，且
[0012][0013]
其中ρ
l
为介质一在界面处的密度极限状态，p
′
l
为介质一在界面处的压强的一阶空间导数，c
l
为介质一在界面处的声速极限状态，u
′
l
为介质一在界面处的速度的一阶空间导数，m＝1,2,3是坐标维数，x
cd
为界面位置坐标，u
l
为介质一在界面处的速度极限状态；
[0014]
2)建立关于介质二的关系式：
[0015][0016]
其中，ρ
2*
为介质二在界面处的密度的预测值，c
2*
为介质二在界面处的声速的预测值，zr为仅和介质二初值相关的参数，且
[0017][0018]
其中ρr为介质二在界面处的密度极限状态，p
′r为介质二在界面处的压强的一阶空间导数，cr为介质二在界面处的声速极限状态，u
′r为介质二在界面处的速度的一阶空间导数，ur为介质二在界面处的速度极限状态；
[0019]
3)联立步骤1)和步骤2)的关系式，获取物质界面处速度随时间的变化率ua和压强随时间的变化率pa：
[0020][0021][0022]
4)获取界面处的介质一的速度的一阶空间导数预测值u
′
1*
、压强的一阶空间导数预测值p
′
1*
，界面处的介质二的速度的一阶空间导数预测值u
′
2*
、压强的一阶空间导数预测值p
′
2*
：
[0023][0024]
p
′
1*
＝-ρ
1*
ua[0025][0026]
p
′
2*
＝-ρ
2*
ua[0027]
其中，u
1*
为介质一在界面处的速度的预测值，u
2*
为介质二在界面处的速度的预测值；
[0028]
5)获取界面处的介质一的密度的一阶空间导数预测值ρ
′
1*
，界面处的介质二的密
度的一阶空间导数预测值ρ
′
2*
：
[0029]
建立介质一关于密度空间导数预测值ρ
′
1*
和介质一在界面处的密度的一阶空间导数ρ
′
l
的关系式：
[0030][0031]
其中和是介质一内能e
l
分别对压强p
l
和密度ρ
l
的偏导数，和是介质一预测内能e
1*
分别对预测压强p
1*
和预测密度ρ
1*
的偏导数，求得：
[0032][0033]
建立介质二关于密度空间导数预测值ρ
′
2*
和介质二在界面处的密度的一阶空间导数ρ
′r的关系式：
[0034][0035]
其中和是介质二内能er分别对压强pr和密度ρr的偏导数，和是介质二预测内能e
2*
分别对预测压强p
2*
和预测密度ρ
2*
的偏导数，求得：
[0036][0037]
6)获取界面两侧介质一的物理量的一阶空间导数预测值u
′
1*
和介质二的物理量的一阶空间导数预测值u
′
2*
：
[0038]u′
1*
:＝(ρ
′
1*
,u
′
1*
,p
′
1*
)
[0039]u′
2*
:＝(ρ
′
2*
,u
′
2*
,p
′
2*
)。
[0040]
具体地，将所获得的介质一的的物理量一阶空间导数预测值和介质二的物理量的一阶空间导数预测值用于对一般状态方程下的带对称性源项多介质界面高阶处理方法中。
[0041]
具体地，所述的一般状态方程下的带对称性源项多介质界面高阶处理方法包括：
[0042]
s1：读取流场数据，获取多介质界面处两侧的物理量的状态值及一阶空间导数，建立界面处的多介质riemann问题以及具有线性初值分布的一般状态方程下的带对称性源项多介质广义riemann问题；
[0043]
s2：使用已有多介质riemann问题求解器求解步骤s1中建立的界面处的多介质riemann问题，获得界面两侧物理量的预测值；
[0044]
s3：求解步骤s1中建立的具有线性初值分布的一般状态方程下的带对称性源项多介质广义riemann问题，获得界面两侧物理量一阶空间导数的预测值；
[0045]
s4：在界面两侧建立两个虚拟区域分别作为两种介质的边界区域，使用步骤s2中获得的界面两侧物理量的预测值和步骤s3中获得的界面两侧物理量一阶空间导数的预测值对两个虚拟区域分别进行线性赋值，将原多介质问题解耦为两个单介质问题，采用时间离散方法在时间上推进，分别计算，得到流场的仿真结果。
[0046]
本发明的有益效果在于：
[0047]
1.给出一种界面两侧物理量一阶空间导数预测值的获得方法，所获得的物理量一阶空间导数预测值尤其适用于对一般状态方程下的带对称性源项多介质界面高阶处理方法中；
[0048]
2.使用该方法能近似求解一般状态方程下的带对称性源项的多介质广义riemann问题，有效降低求解成本的同时，仍能保持理论精度，具有简单实用有效的特点。
[0049]
3.相对cn114254573a公开的处理方法，采用本发明所得到的界面两侧物理量一阶空间导数预测值，并进行带对称性源项多介质界面高阶处理，不再局限于指定的特殊状态方程，可以广泛应用于处理一般状态方程下的多介质耦合问题，具有很高的普适性。
附图说明
[0050]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，可以根据这些附图获得其他的附图。其中：
[0051]
图1为一维球对称空泡爆炸算例示意图；
[0052]
图2为本方法计算一维球对称空泡爆炸算例得到的压力曲线与未使用高阶方案、使用原始特殊状态方程下高阶方案的对比。
[0053]
图3为一维球对称膨胀活塞算例示意图；
[0054]
图4为本方法计算一维球对称膨胀活塞算例得到的密度曲线与未使用高阶方案、使用原始特殊状态方程下高阶方案的对比。
具体实施方式
[0055]
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在不冲突的情况下，本发明的实施例及实施例中的特征可以相互组合。
[0056]
在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。
[0057]
本发明涉及的一种界面两侧物理量一阶空间导数预测值的获得方法，所获得的物理量一阶空间导数预测值适用于对一般状态方程下的带对称性源项多介质界面高阶处理方法中。所述适用于一般状态方程下的带对称性源项多介质界面高阶处理方法基本流程与cn114254573a所使用的基本相同，均包括以下步骤：
[0058]
s1：读取流场数据，获取多介质界面处两侧的物理量的状态值及一阶空间导数，建立界面处的多介质riemann问题以及具有线性初值分布的一般状态方程下的带对称性源项多介质广义riemann问题；具体过程为：
[0059]
s1-1：带对称性源项的可压缩欧拉方程为：
[0060]
[0061][0062]
其中，ρ,u,p,e分别为密度，速度，压强和内能，x为径向坐标，t为时间，m＝1,2,3是坐标维数，u(ρ,u,p)为守恒变量，f(u)为通量，ψ(x,u)为对称性源项，e是单位体积总能量，一般状态方程可统一写成e＝e(ρ,p)，c是声速且大于0，
[0063]
s1-2：t0时刻界面处的多介质riemann问题描述为：
[0064][0065][0066]
其中，u(x,t)为t时刻坐标x处的状态表达式，x
cd
(t0)为t0时刻的物质界面的位置，x《x
cd
(t0)表示界面左侧，为介质一，x》x
cd
(t0)表示界面右侧，为介质二；
[0067]
ρ
l
,u
l
,p
l
,c
l
,e
l
是介质一在界面处的密度、速度、压强、声速和内能的极限状态，记u
l
＝u(ρ
l
,u
l
,p
l
)；
[0068]
ρr,ur,pr,cr,er是介质二在界面处的密度、速度、压强、声速和内能的极限状态，记ur＝u(ρr,ur,pr)；
[0069]
p
eos
(ρ,e)为状态方程，p
eos,1
(ρ,e)和p
eos,2
(ρ,e)为对应于介质一和介质二的一般状态方程；
[0070]
s1-3：t0时刻界面处的具有线性初值分布的一般状态方程下的带对称性源项多介质广义riemann问题描述为：
[0071][0072][0073]
其中，ρ
′
l
,u
′
l
,p
′
l
是介质一在界面处的密度、速度和压强的一阶空间导数，记u
′
l
:＝(ρ
′
l
,u
′
l
,p
′
l
)；
[0074]
ρ
′r,u
′r,p
′r是介质二在界面处的密度、速度和压强的一阶空间导数，记u
′r:＝(ρ
′r,u
′r,p
′r)。
[0075]
s2：使用已有多介质riemann问题求解器求解步骤s1中建立的界面处的多介质riemann问题，获得界面两侧物理量的预测值；ρ
1*
,u
1*
,p
1*
,c
1*
,e
1*
为介质一在界面处的密度、速度、压强、声速和内能的预测值，记u
1*
＝u(ρ
1*
,u
1*
,p
1*
)；ρ
2*
,u
2*
,p
2*
,c
2*
,e
2*
为介质二在界面处的密度、速度、压强、声速和内能的预测值，记u
2*
＝u(ρ
2*
,u
2*
,p
2*
)。
[0076]
s3：求解步骤s1中建立的具有线性初值分布的一般状态方程下的带对称性源项多介质广义riemann问题，获得界面两侧物理量一阶空间导数的预测值；包括以下步骤：
[0077]
s3-1：建立关于介质一的关系式：
[0078][0079]
其中，ua为界面处速度随时间的变化率，pa为界面处压强随时间的变化率，z
l
为仅和介质一初值相关的参数，
[0080]
s3-2：建立关于介质二的关系式：
[0081][0082]
其中，zr为仅和介质二初值相关的参数，
[0083]
s3-3：联立步骤s3-1和步骤s3-2的关系式，获取物质界面处速度随时间的变化率ua和压强随时间的变化率pa：
[0084][0085][0086]
s3-4：获取界面处的介质一的速度的一阶空间导数预测值u
′
1*
、压强的一阶空间导数预测值p
′
1*
，界面处的介质二的速度的一阶空间导数预测值u
′
2*
、压强的一阶空间导数预测值p
′
2*
：
[0087][0088]
p
′
1*
＝-ρ
1*
ua[0089][0090]
p
′
2*
＝-ρ
2*
ua[0091]
s3-5：获取界面处的介质一的密度的一阶空间导数预测值ρ
′
1*
，界面处的介质二的密度的一阶空间导数预测值ρ
′
2*
：
[0092]
建立介质一关于密度空间导数预测值ρ
′
1*
和ρ
′
l
的关系式：
[0093][0094]
其中和是介质一内能e
l
分别对压强p
l
和密度ρ
l
的偏导数，和是介质一预测内能e
1*
分别对预测压强p
1*
和预测密度ρ
1*
的偏导数，求得
[0095]
建立介质二关于密度空间导数预测值ρ
′
2*
和ρ
′r的关系式：
[0096][0097]
其中和是介质二内能er分别对压强pr和密度ρr的偏导数，和是介质二预测内能e
2*
分别对预测压强p
2*
和预测密度ρ
2*
的偏导数，求得
[0098]
s3-6：获取界面两侧介质一和介质二的物理量的一阶空间导数预测值：
[0099]u′
1*
:＝(ρ
′
1*
,u
′
1*
,p
′
1*
)，
[0100]u′
2*
:＝(ρ
′
2*
,u
′
2*
,p
′
2*
)。
[0101]
s4：在界面两侧建立两个虚拟区域分别作为两种介质的边界区域，使用步骤s2中获得的界面两侧物理量的预测值和步骤s3中获得的界面两侧物理量一阶空间导数的预测值对两个虚拟区域分别进行线性赋值，将原多介质问题解耦为两个单介质问题，采用时间离散方法在时间上推进，分别计算，得到流场的仿真结果。
[0102]
介质一的虚拟区域线性赋值为：
[0103]u1*
+(x-x
cd
)u
′
1*
:＝u(ρ
1*
+(x-x
cd
)ρ
′
1*
,u
1*
+(x-x
cd
)u
′
1*
,p
1*
+(x-xcdp1*
′
，x》xcd(t0)；
[0104]
介质二的虚拟区域线性赋值为：
[0105]u2*
+(x-x
cd
)u
′
2*
:＝u(ρ
2*
+(x-x
cd
)ρ
′
2*
,u
2*
+(x-x
cd
)u
′
2*
,p
2*
+(x-xcdp2*
′
，x《xcd(t0)。
[0106]
为了方便理解本发明的上述技术方案，以下通过具体实施例对本发明的上述技术方案进行详细说明。
[0107]
实施例1
[0108]
以一维球对称可压缩水中空泡爆炸过程为例，如图1所示，其中内部球为超高压气体，外部为常压可压缩液体水；对于气体和水，状态方程可以写成形式，γ为比热比，pb为压强常数，二者均和介质相关。
[0109]
总计算区域为12的一维球对称区域，设x为空间坐标，即x∈(0,12),网格为3000个点均匀分布，初始气泡半径x0＝0.401；
[0110]
气体速度初始值为u
gas
＝0，压强初始值为p
gas
＝8290.91，密度初始值为ρ
gas
＝1.27，γ
gas
＝1.4为理想气体比热比，p
b,gas
＝0为气体压强常数。
[0111]
可压液体水速度初始值为u
water
＝0，压强初始值为p
water
＝1.0，密度初始值为ρ
water
＝1.0，γ
water
＝7.0为水比热比，p
b,water
＝3000.0为水压强常数。
[0112]
为了说明本发明的方法的效果，对t＝0.07秒压强的分布图，结果如图2所示，其中三角形实线代表没有应用本发明方法的计算结果，矩形虚线代表应用原始针特殊状态方程下的高阶方案的计算结果，菱形虚线代表应用本发明方法的计算结果，可见，使用本发明的方法处理一般状态方程下的带对称性源项多介质耦合问题时，能达到与原始高阶方案相同的效果，能有效地消除界面两侧的压强差，能够保持界面处压强的连续性。
[0113]
实施例2
[0114]
本算例为球形活塞对气体的冲击膨胀问题，如图3所示，是一个典型的球对称系统下，给定速度、加速度的运动界面问题。外部流场为完全气体，流场状态为(u0,p0,ρ0)＝(0,1.0e+4,1.0)，满足状态方程p＝(γ-1)ρe，γ＝7；中间为膨胀球形活塞，初始时间为t0＝3.0e-5，初始位置为r
b0
＝30，在其上附着有一个激波，即激波初始半径为s0＝30，已知其界面位置rb，速度ub满足
[0115]
为了说明本发明的方法的效果，对t＝0.003889秒密度的分布图，结果如图4所示，其中黑色粗实线代表该问题的精确解，三角形虚线代表没有应用高阶方案的计算结果，矩形虚线代表应用原始特殊状态方程下的高阶方案的计算结果，菱形虚线代表应用本方案的计算结果，可见，使用本发明的方法处理一般状态方程下的带对称性源项界面问题时，能有效地消除激波错位现象，能够保持界面处物理性质的连续性。
[0116]
针对多介质的情况，即每相邻的两个介质界面，都采用本发明的方法进行计算。
[0117]
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
技术特征：1.一种多介质界面两侧物理量一阶空间导数预测值计算方法，其特征在于，通过求解具有线性初值分布的一般状态方程下的带对称性源项多介质广义riemann问题，获得界面两侧物理量一阶空间导数的预测值；包括以下步骤：1)建立关于介质一的关系式：其中，u
a
为界面处速度随时间的变化率，p
a
为界面处压强随时间的变化率，ρ
1*
为介质一在界面处的密度的预测值，c
1*
为介质一在界面处的声速的预测值，z
l
为仅和介质一初值相关的参数，且：其中ρ
l
为介质一在界面处的密度极限状态，p
′
l
为介质一在界面处的压强的一阶空间导数，c
l
为介质一在界面处的声速极限状态，u
′
l
为介质一在界面处的速度的一阶空间导数，m＝1,2,3是坐标维数，x
cd
为界面位置坐标，u
l
为介质一在界面处的速度极限状态；2)建立关于介质二的关系式：其中，ρ
2*
为介质二在界面处的密度的预测值，c
2*
为介质二在界面处的声速的预测值，z
r
为仅和介质二初值相关的参数，且其中ρ
r
为介质二在界面处的密度极限状态，p
′
r
为介质二在界面处的压强的一阶空间导数，c
r
为介质二在界面处的声速极限状态，u
′
r
为介质二在界面处的速度的一阶空间导数，u
r
为介质二在界面处的速度极限状态；3)联立步骤1)和步骤2)的关系式，获取物质界面处速度随时间的变化率u
a
和压强随时间的变化率p
a
：：4)获取界面处的介质一的速度的一阶空间导数预测值u
′
1*
、压强的一阶空间导数预测值p
′
1*
，界面处的介质二的速度的一阶空间导数预测值u
′
2*
、压强的一阶空间导数预测值p
′
2*
：p
′
1*
＝-ρ
1*
u
a
p
′
2*
＝-ρ
2*
u
a
其中，u
1*
为介质一在界面处的速度的预测值，u
2*
为介质二在界面处的速度的预测值；5)获取界面处的介质一的密度的一阶空间导数预测值ρ
′
1*
，界面处的介质二的密度的一阶空间导数预测值ρ
′
2*
：建立介质一关于密度空间导数预测值ρ
′
1*
和介质一在界面处的密度的一阶空间导数ρ
′
l
的关系式：其中和是介质一内能e
l
分别对压强p
l
和密度ρ
l
的偏导数，和是介质一预测内能e
1*
分别对预测压强p
1*
和预测密度ρ
1*
的偏导数，求得：建立介质二关于密度空间导数预测值ρ
′
2*
和介质二在界面处的密度的一阶空间导数ρ
′
r
的关系式：其中和是介质二内能e
r
分别对压强p
r
和密度ρ
r
的偏导数，和是介质二预测内能e
2*
分别对预测压强p
2*
和预测密度ρ
2*
的偏导数，求得：6)获取界面两侧介质一的物理量的一阶空间导数预测值u
′
1*
和介质二的物理量的一阶空间导数预测值u
′
2*
：u
′
1*
:＝(ρ
′
1*
,u
′
1*
,u
′
1*
）u
′
2*
:＝(ρ
′
2*
,u
′
2*
,u
′
2*
)。2.根据权利要求1所述的一种多介质界面两侧物理量一阶空间导数预测值计算方法，其特征在于，将所获得的介质一的物理量的一阶空间导数预测值和介质二的物理量的一阶空间导数预测值用于对一般状态方程下的带对称性源项多介质界面高阶处理方法中。3.根据权利要求2所述的一种多介质界面两侧物理量一阶空间导数预测值计算方法，其特征在于，所述的一般状态方程下的带对称性源项多介质界面高阶处理方法包括：s1：读取流场数据，获取多介质界面处两侧的物理量的状态值及一阶空间导数，建立界面处的多介质riemann问题以及具有线性初值分布的一般状态方程下的带对称性源项多介质广义riemann问题；s2：使用已有多介质riemann问题求解器求解步骤s1中建立的界面处的多介质riemann问题，获得界面两侧物理量的预测值；s3：求解步骤s1中建立的具有线性初值分布的一般状态方程下的带对称性源项多介质广义riemann问题，获得界面两侧物理量一阶空间导数的预测值；
s4：在界面两侧建立两个虚拟区域分别作为两种介质的边界区域，使用步骤s2中获得的界面两侧物理量的预测值和步骤s3中获得的界面两侧物理量一阶空间导数的预测值对两个虚拟区域分别进行线性赋值，将原多介质问题解耦为两个单介质问题，采用时间离散方法在时间上推进，分别计算，得到流场的仿真结果。
技术总结本发明公开了一种多介质界面两侧物理量一阶空间导数预测值计算方法，包括：分别建立界面两侧关于建立关于介质一、二的关系式并联立得到物质界面处速度随时间的变化率和压强随时间的变化率，获取界面处的介质一、二的速度的一阶空间导数预测值、压强的一阶空间导数预测值，获取界面处的介质一、二的密度的一阶空间导数预测值，并最终获取界面两侧介质一、二的一阶空间导数预测值。采用本发明所得到的界面两侧物理量一阶空间导数预测值可用于一般状态方程下的带对称性源项多介质界面高阶处理，且不再局限于指定的特殊状态方程，可以广泛应用于处理一般状态方程下的多介质耦合问题，具有很高的普适性。具有很高的普适性。具有很高的普适性。
技术研发人员：刘铁钢 张孝涛 冯成亮 张斌
受保护的技术使用者：北京航空航天大学
技术研发日：2022.07.26
技术公布日：2022/11/1
