一种基于多维黎曼解的任意拉格朗日欧拉方法
【专利摘要】本发明公开了一种基于二维黎曼解的任意拉格朗日欧拉方法，用以求解二维可压缩流体力学方程组及以之为基础的相关问题。整个算法的流程包括：网格生成与物理量的初始分布；确定下一时刻网格的策略和方式；确定网格单元边界流量的二维黎曼解法器算法；网格单元边界流量的表达形式；获得下一时刻的计算结果的数值方法。新的发明提供了一种确定网格单元边界流量的新的二维黎曼解法器算法，能够克服传统的一维黎曼解法导致的数值激波不稳定、网格易扭曲等现象，改正了现有二维黎曼解方法过于复杂而难以实施的缺点。本发明是一种简洁、健壮和精确的数值算法组件，适合多介质大变形问题的数值模拟，适合现在流行的有限差分、有限体积和有限元方法的形式。
【专利说明】一种基于多维黎曼解的任意拉格朗日欧拉方法
【技术领域】
[0001]本发明属于流体力学数值模拟【技术领域】，涉及一种基于多维黎曼解的任意拉格朗日欧拉方法，具体地说，涉及一种利用二维黎曼解模拟多介质可压缩流动的任意拉格朗日欧拉方法。
【背景技术】
[0002]多介质大变形问题的数值模拟，是流体力学问题数值模拟中具有挑战性的课题。它需要求解涉及多种介质相互作用、并发生剧烈运动的问题(如，介质撞击和作用的弹塑性问题，爆轰问题，高温高压的辐射流体力学问题)。在这类问题中，介质常常具有复杂的物态性质，状态量(密度、压强)会发生突然的改变(激波、爆轰波出现)，作用区域有可能会发生几百倍的变化。这种问题的计算方法包含了复杂的流程，如网格运动、介质界面的追踪和可压缩流体力学方程组的数值解法。
[0003]在工程中目前应用的主要方法有欧拉方法、拉格朗日方法和任意拉格朗日欧拉方法。欧拉方法使用固定网将，实施简便。主要困难一个是难以精确地计算区域变动巨大的问题，另一个是难以正确地描述不同介质处于同一个网格内的状态(特另是当介质的物态性质相差较大、或者存在辐射等多物理作用过程时)。在弹塑性问题和辐射流体力学问题中，人们较少使用这一方法，更多的是采用拉格朗日方法和任意拉格朗日欧拉方法。拉格朗日方法使用的网格随着介质的运动而运动，可以清晰地分辨介质界面，不需要计算介质在人为混合后的状态。但这种方法非常容易造成网格扭曲，导致计算中断。任意拉格朗日欧拉方法可以在计算过程中调节网格的运动方式，例如仅在介质界面处让网格随着介质运动，而在其它地方则调整网格，使之保持较好的几何品质。该方法在理论上同时兼具欧拉与拉格朗日方法的优点，是一种较理想的解决方案。但在实际的计算过程中，计算结果还是会出现严重的错误，如虚假的涡度误差和数值激波不稳定现象。这些错误扭曲了真实的物理图像，并常常导致计算过程的中断。计算出现错误的主要原因是流体力学方程组的解法器出现了问题。
[0004]流体力学方程组的解法器是一种利用当前时刻的物理量，获得下一时刻的物理量的算法实施组件。其关键之处是设计当流体通过网格边界时数值通量(流量)的算法。目前在工程应用中广泛使用的流体力学方程组解法器，都是沿网格边界法线方向，求解一维的黎曼问题，并以获得的结果作为流量的近似值使用的。这种方法对一般流体力学问题是适用的，但在多介质大变形等极端条件下，由于缺少多维信息，对多维效应引起的涡度误差和数值激波不稳定性等难点问题难以解决。依靠加密计算网格，减少时间步长，高精度技术等传统方法并不能改正这些缺陷。有时候加密网格还导致计算终止得更快。采用包含多维信息的多维黎曼解法器，是解决这类问题的可行方案。
[0005]遗憾的是，针对二维黎曼问题的解法器还没有在任何这类问题的应用软件中实施。这是由于二维黎曼问题的精确求解理论上还没有完全解决，即使使用近似解法器方法，实施过程也非常复杂。现有的研究成果，基本都停留在学术层面，大部分仅仅应用在四边形网格的欧拉方法中。法国的Despres和Maire等人，分别设计了一种可以在非结构网格(六边形或者三角形)上使用的二维近似黎曼解法器，但仅仅适用于拉格朗日方法。现有的任意拉格朗日欧拉方法在网格边界的通量计算中，还没有实施简单、适用于复杂网格体系、健壮有效的二维黎曼解法器。
【发明内容】
[0006]为了克服现有技术中的缺陷，本发明提供了一种利用二维黎曼解模拟多介质可压缩流动的任意拉格朗日欧拉方法，针对复杂的流体力学计算问题，研制的方法能适应复杂的网格体系，具备分辨多介质问题的能力，并取得良好的计算效果。与传统方法使用局部一维的黎曼解方法不同的是，本发明设计了一个两维的黎曼问题解法器，利用多维信息改正传统方法的不足。本发明是一种健壮、简洁的数值算法，可以明显地减少传统方法的误差。其技术方案为，
[0007]一种基于多维黎曼解的任意拉格朗日欧拉方法，包括以下步骤:
[0008]I)网格生成与物理量的初始分布
[0009]所有的物理量(密度、速度、压强、能量)均定义在网格的中心。
[0010]2)下一时刻网格的确定和生成
[0011]针对不同的网格策略(欧拉方法，自适应网格方法和拉格朗日方法)，求解特定的网格方程，并形成软件程序，主要包括椭圆型网格生成器，自适应网格生成器等。
[0012]3) 二维黎曼解法器组件
[0013]传统方法缺陷:工程应用中的一维黎曼解算法，不具备二维黎曼解的性质。文献中已有的二维黎曼解算法需要利用前一个时刻的物理量，求解复杂的波系相互作用情况。由于理论不完备，这造成算法的实施非常困难。
[0014]新的发明想法:如果二维黎曼解法器仅仅提供经过网格边界的流量算法，而不计算复杂波系的相互作用，就可以大大减少问题的难度。
[0015]新的发明遇到的困难:如何从前一个时刻网格中心的物理量，获得直到下一个时刻的网格节点和网格边界的物理量。如果采用直接的插值等数值技术，由于物理量存在间断，将会完全失效。
[0016]新的发明主要内容包括:节点运动速度的确定方法和网格边界通量的计算方法。
[0017]I)节点运动速度的确定方法
[0018]在某一网格V。的所有边界上(看图1(a))，求解沿着边界qq+法方向的局部一维黎曼问题，获得网格边界上的运动速度。当围绕所要求解的节点< (看图1(b))的所有边界上的边运动速度全部求出时，利用这些边速度的法方向分量，采用最小二乘方法求出网格节点的运动速度 <。其具体规则是节点的运动速度在网格边的投影与这些边速度的法方
向分量的差，在最小二乘意义下最小。计算过程使用的权函数，需要特殊的设计。
[0019]2)网格边界通量的计算方法
[0020]利用节点处的二维流体运动速度< ;将此速度在网格边界qq+处的投影，更新为
相应边界的新的流动速度，再以两边的状态和此流动速度，设计沿qq+法向的新的一维黎曼解。此时由于网格节点处的运动速度，进入了边界通量的计算中，传统一维黎曼解中的间断跳跃条件不再成立。一种单边的新的跳跃条件被构造出来，一个边界上具有了两个流量(看图2(b))，分别对应不同的网格单元。这种方法能保证网格的运动方式与边界通量完全相容。沿边界的切向动量，也由于引进了节点速度的平均作用，导致计算更加稳定。
[0021]3)边界通量表达形式
[0022]给定了两层网格节点的移动速度后，根据黎曼解的结构，设计出边界通量的表达形式。4)获得下一时刻的计算结果
[0023]如果计算时间已经到了需要的时刻，则计算终止，否则回到第2)步。
[0024]本发明的有益效果:
[0025]本发明适合复杂的网格，特别是无结构网格及结构和非结构混合网格；清晰地分辨介质界面；明显地改善现有流体力学算法的健壮性，提高计算精度。
【专利附图】
【附图说明】
[0026]图1是网格及其节点的几何示意图，其中图1(a)是网格单元V。及相邻单元，图1(b)是网格节点及相邻单元；
[0027]图2是构造边界流量的算法示意图，其中图2(a)是传统算法，图2(b)是新提出的边界流量算法；
[0028]图3是新发明的局部一维黎曼解的波系分解图，其中图3(a)是传统算法，图3 (b)是新提出的边界通量；
[0029]图4是激波衍射问`题密度轮廓图，其中图4 (a)是采用传统的HLLC —维近似黎曼解算法，图4(b)是采用新的二维黎曼解算法；
[0030]图5是初始网格一致倾斜时(Saltzmann网格)，其中图5 (a)是采用传统的HLLC一维近似黎曼解算法，图5(b)是采用新的二维黎曼解算法，
[0031]图6是初始网格为三角形网格的双马赫问题，其中图6(a)是采用传统的HLLC —维近似黎曼解算法，图6(b)是采用新的二维黎曼解算法；
[0032]图7是初始网格为三角形网格的双马赫问题，其中图7(a)是为密度轮廓，图7(b)是为网格图；
[0033]图8是初始网格为三角形网格的多介质激波管问题，图8(a)为初始网格，图8 (b)为计算过程中的拉格朗日网格，图8 (c)为计算过程中的自适应移动网格；图8 (d)是沿着横向的密度分布，图8(e)是沿着横向的压强分布。
【具体实施方式】
[0034]下面结合附图和【具体实施方式】对本发明的技术方案作进一步详细地说明。
[0035]使用一种基于多维黎曼解的任意拉格朗日欧拉方法，求解如下的流体力学方程组
謂 U) f)G(U)
[0036]— + -十」+ —=Hs
ο? axay
[0037]其中H表示多物理过程的源项，U，F，G分别是
【权利要求】
1.一种基于多维黎曼解的任意拉格朗日欧拉方法，其特征在于，包括以下步骤: 1)网格生成与物理量的初始分布； 2)下一时刻网格的确定和生成:需要求解与求解特定的网格方程，并形成软件程序，主要包括椭圆型网格生成器，自适应网格生成器； 3)黎曼解法器组件 构造一个全新的两维近似黎曼解法，该解法器包括节点速度的确定方法和网格边界通量的计算方法，在任意网格下，网格节点的速度算法:已知网格单元中心的物理量，需要求解2X2的代数方程组来获得节点处的流动速度，点处的流动速度与网格边界通过一些特定方法获得的法向速度差在最小二乘意义下最小； 构造边界流量的方法:利用了节点处的二维流体速度，以此速度在网格边界qq+处的投影，作为相应边界的流动速度，再以两边的状态和此流动速度；构造沿qq+法向的一维黎曼解，由于网格节点处的运动速度，进入了边界通量的计算中，这种方法能保证网格的运动方式与边界通量完全相容，沿边界的切向动量，也引进了节点速度的作用； 4)边界通量表达形式:给定了两层网格节点的移动速度后，根据黎曼解的结构，写出边界通量的表达形式； 5)获得下一时刻的计算结果 如果计算时间已经到了需要的时刻，则计算终止，否则回到第2)步。
【文档编号】G06F17/50GK103823916SQ201310498537
【公开日】2014年5月28日 申请日期:2013年10月23日 优先权日:2013年10月23日
【发明者】沈智军, 闫伟, 袁光伟 申请人:沈智军, 闫伟, 袁光伟