CN104699909A - 一种保持强稳定性的变步长多步法时间离散算法 - Google Patents
一种保持强稳定性的变步长多步法时间离散算法 Download PDFInfo
- Publication number
- CN104699909A CN104699909A CN201510132095.7A CN201510132095A CN104699909A CN 104699909 A CN104699909 A CN 104699909A CN 201510132095 A CN201510132095 A CN 201510132095A CN 104699909 A CN104699909 A CN 104699909A
- Authority
- CN
- China
- Prior art keywords
- lambda
- discrete
- algorithm
- step size
- variable step
- 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
Links
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Complex Calculations (AREA)
Abstract
本发明提出了一种保持强稳定性的变步长多步法时间离散算法,其包括半离散Hermite本质无振荡有限体积算法的构造,初始步使用了三阶保强稳定的Runge-Kutta算法;本发明给出了两类时间离散的接口处理,特别是初始步的时间步长的选取。本发明的优点在于变时间步长和保强稳定性。通过实例展示了这个新算法在实例中的计算机模拟的高精度、高分辨率以及时间高效性等优势。这种技术还可以与多种空间算法组合,适用于可压缩流体模拟,最终能够推广到实际复杂流场的模拟等大规模科学计算问题。
Description
技术领域
本发明涉及计算流体力学数值方法领域,特别是一种保持强稳定性的变步长多步法时间离散算法。
背景技术
在可压缩流体计算中,欲模拟的物理量通常快速变化,在狭窄区域内以大梯度变化,为准确计算机模拟流场的细节,许多高分辨技术被不断地提出和完善。
目前,求解的双曲守恒律的主流的高精度技术为间断Galerkin有限元方法,有限体积或者有限差分加权本质无振荡方法等技术,通常的这些技术在求解欧拉方程组时需要与保强稳定性Runge-Kutta方法,或者Lax-Wendroff方法相结合。最近,有人给出Lax-Wendroff方法与间断Galerkin有限元方法等高精度方法相结合的高精度算法设计,这些算法在计算时间效率上优于保强稳定性Runge-Kutta方法,然而算法实现比保强稳定性Runge-Kutta方法复杂。
又有人提出易于算法实现的总变差稳定多步法离散技术,与保强稳定性Runge-Kutta方法相比,总变差稳定多步法离散技术在计算效率上可以设计出比保强稳定性Runge-Kutta方法更高效率的方法,然而现有的保强稳定性的多步法时间离散技术基于等步长,不能用于实际问题计算。
发明内容
本发明的主要目的在于克服现有技术中的上述缺陷,提出一种应用于计算流体力学的实际问题的模拟的保持强稳定性的变步长多步法时间离散算法。
本发明采用如下技术方案:
一种保持强稳定性的变步长多步法时间离散算法,令现有的半离散有限体积Hermite加权本质无振荡算法表示为L(Q),并得到其对应的半离散方程其特征在于:令tn为时间节点,n为时间步,Δtn=tn-tn-1;给出对于所述半离散方程的以下离散算法
优选的,还包括另一离散算法
h1=Δtn-3,h2=Δtn-2,h3=Δtn-1,h4=Δtn,h5=Δtn+1,
优选的,还包括另一离散算法
其中h1=Δtn-4,,h2=Δtn-3,h3=Δtn-2,h4=Δtn-1,h5=Δtn,h6=Δtn+1;
优选的,所述离散算法的初始步使用经典保强稳定的Runge-Kutta时间离散技术,即:
初始步的时间步长可以选取正常计算的时间步长,其中n为整数。
由上述对本发明的描述可知,与现有技术相比,本发明具有如下有益效果:
本发明算法的时间步长可以根据发展方程的情况来自动调整时间步长;其中提供三种离散算法以供选择,而方法二和方法三的CFL数比方法一大,数值计算中CFL数越大,计算效率越高,因此这两种技术还具有更高的计算效率。另外,本发明所提出的算法可以与其它主流的守恒律的高精度技术结合。
附图说明
图1是30°斜坡示意图;
图2是本发明算法的整体流程示意图;
图3是本发明算法模拟双马赫问题的数值结果(网格2400*600等分的数值结果);
图4是图3的局部放大图。
具体实施方式
以下通过具体实施方式对本发明作进一步的描述。
本发明对于现有的半离散有限体积Hermite加权本质无振荡算法,我们用记号L(Q)表示这部分算法,记号解释L是一种算法,而Q是算法L所要处理的变量。
一、我们叙述L的算法过程:
算法L用于模拟无粘绝热空气动力学,为了清楚描述所模拟的空气动力学物理过程,我们以二维气体动力学问题为例,我们知道问题必然满足三大物理基本定律,其数学模型为
Ut+F(U)x+G(U)y=0;
其中
上述方程组中第一方程为质量守恒定律,第二个和第三个方程为动量守恒定律,第四个方程为能量守恒定律。相应的变量符号说明:单位体积的总能ρ,(u,v),P分别为流体的密度,速度,压强;γ=1.4为理想气体状态系数。
为了方便地描述模拟上述物理过程的算法,将其方程组用下面一般的符号描述双曲守恒律
符号说明:div(ab)=ax+by。
对上面的守恒律方程组作一个梯度grad(我们用在直角坐标下说明这个符号, 得到
其中的含义为
而Hermite有限体积算法的控制方程由上面两个方程(1)和(2)组成
Qt+div H(Q)=0,
其中Q=(U,grad U)T,
对上面的控制方程组在小区域Ωj(物理表面记为)上积分,
再通过Gauss公式,把积分转化为带方向的积分可以得到半离散有限体积方法
其中|Ωj|是区域Ωj的体积,ds可以理解为对某一物体表面细分后的面,对应的这个面的向外的方向为n。
另由S个组成,即那么对上面方程右端的积分就可以用在各个上的离散点数为q的高斯积分离散
其中:Gsl和ωl是相应的高斯积分点和权重。(q点高斯积分离散指的是 )
对于在某一高斯点的数值流通量H(Q(Gsl,t))·n可由近似或者精确Riemann解代替(这是由不同之间按照某一公式计算得到,这里算法不是本发明所要考虑的,不再详细叙述)。
到此为止,这部分关于算法L(L就是公式(3)中)叙述完了,这是因为上述过程中其实是基于的公式过程。因此这里我们记为简记为L(Q)。叙述完算法L,得到一个半离散方程
二、对于本发明提出的一种保持强稳定性的变步长多步法时间离散算法,给出了这类方法中的3种技术:
记号说明:tn是时间节点,n为时间步,Δtn=tn-tn-1。
方法一、保强稳定性的变步长多步法时间离散技术:
其中h1=Δtn-2,h2=Δtn-1,h3=Δtn,h4=Δtn+1和可以看到上面的时间步长h为变时间,如果Δtn=Δt,那么这一个方法变为等时间步长多步法
这个方法的CFL数为1/3,因此本发明的多步法技术的CFL数也取为1/3。
CFL数是算法稳定相关的参数,CFL数越大,则计算时间效率越好。
方法二、保强稳定性的变步长多步法时间离散技术:
其中h1=Δtn-3,h2=Δtn-2,h3=Δtn-1,h4=Δtn,h5=Δtn+1和这种技术的CFL数可取为1/2,同样选取相对应等步长方法的CFL数。
方法三、保强稳定性的变步长多步法时间离散技术:
其中h1=Δtn-4,h2=Δtn-3,h3=Δtn-2,h4=Δtn-1,h5=Δtn,h6=Δtn+1和这种技术的CFL数可取为0.567,同样选取相对应等步长方法的CFL数。
如上面所述的方法一公式中得到Qn+1时需要Qn和Qn-3。我们计算从n=0开始,因此无法直接使用多步法,需要使用以下的保强稳定性Runge-Kutta方法开得到Q0,Q1,Q2,Q3。
初始步的时间步长可以选取正常计算的时间步长,符号说明本说明书中n作为上下标表示为整数。
以下使用本发明算法模拟计算流体力学中的双马赫问题:
1.问题描述:
这是二维气体动力学间断解问题,其数学模型为
Ut+F(U)x+G(U)y=0,
其中
而单位体积的总能其中ρ,(u,v),P分别为流体的密度,流体速度,压强,理想气体状态系数γ=1.4。
初始时刻,Mach 10的正激波撞向一个30°斜坡,如图1。那么激波与斜坡之间成60°夹角。使用所发明的离散技术计算流场的演化,模拟出tprint=0.2时刻流场的密度分布。
2.计算机模拟参数及边界条件:
斜坡平面作为x轴,建立平面直角坐标系,计算区域取为[0,4]×[0,1],反射壁面处于计算区域的内部,一个Mach数为10的斜强激波放置在处,并与x轴成60°夹角。在之前的底边壁面采用准确的激波波后条件,其他壁面采用反射边界条件。左侧波后区域的边界条件为入流边界条件,右侧为出流边界条件。上边界为精确边界条件即给出斜激波。
3.算法L过程:
如图2整体流程图的第一步,初值条件的赋值,时间层nnt=0,tnum=0:
1)对计算区域划分:
把计算区域[0,4]×[0,1]作在x轴和y轴上分别做2400等分和600等分,对划分单元编号,x轴和y轴的编号分别为i和j,则这些划分后的单元可用Iij表示,单元关于坐标轴的长高分别为Δx,Δy。
2)赋初始条件,即赋计算量当y>(x-1/6)tan(π/3)时,赋波后条件,否则赋波前条件。
每一步时间计算的边界条件:
3)赋边界条件:边界条件给定见上面2的描述。例如上边界条件,当时,赋波后条件,否则赋波前条件。
赋初始条件和边界条件后,求近似或者精确Riemann解,即说明中的公式(3);到此完成了算法L(L(Q))。
4.多步法时间离散算法过程,这里叙述方法一的过程,方法二和方法三的过程与方法一一样,不一样之处在于他们的计算公式不同:
若采取上文中方法一,初始步数nnt为4,即需要保强稳定性Runge-Kutta方法(4)求出Q1,Q2,Q3,如图2,当步数nnt<初始步数,保强稳定性Runge-Kutta方法执行初始步。nnt为时间计算步,tnum为计算时间,dt为计算步长度,tprint为最后的计算时间。
注意每执行完保强稳定性Runge-Kutta方法的中间步后要赋3)中的边界条件,时间更新tnum=tnum+dt,其中dt为Δt。
执行完初始步后,进行方法一
其中h1=Δtn-2,h2=Δtn-1,h3=Δtn,h4=Δtn+1和
反复通过上面公式计算从Qn到Qn+1,直到达到指定时间tprint=0.2,注意每执行完公式(7)就要赋3)中的边界条件,时间更新tnum=tnum+dt,其中dt为Δtn+1。
结果性能分析:
图3、图4为网格2400*600等分的数值结果,从图3可以看出发明的技术很好的捕捉了激波,从图4可以看出发明可以模拟出详细的涡结构。
本发明的优点在于变时间步长和保强稳定性。通过实例展示了这个新算法在实例中的计算机模拟的高精度、高分辨率以及时间高效性等优势。这种技术还可以与多种空间算法组合,适用于可压缩流体模拟,最终能够推广到实际复杂流场的模拟等大规模科学计算问题。
上述仅为本发明的具体实施方式,但本发明的设计构思并不局限于此,凡利用此构思对本发明进行非实质性的改动,均应属于侵犯本发明保护范围的行为。
Claims (4)
1.一种保持强稳定性的变步长多步法时间离散算法,令现有的半离散有限体积Hermite加权本质无振荡算法表示为L(Q),并得到其对应的半离散方程其特征在于:令tn为时间节点,n为时间步,Δtn=tn-tn-1;给出对于所述半离散方程的以下离散算法
h1=Δtn-2,h2=Δtn-1,h3=Δtn,h4=Δtn+1;CFL数可取为1/3。
2.如权利要求1所述的一种保持强稳定性的变步长多步法时间离散算法,其特征在于:还包括另一离散算法
h1=Δtn-3,h2=Δtn-2,h3=Δtn-1,h4=Δtn,h5=Δtn+1,
3.如权利要求1所述的一种保持强稳定性的变步长多步法时间离散算法,其特征在于:还包括另一离散算法
其中h1=Δtn-4,h2=Δtn-3,h3=Δtn-2,h4=Δtn-1,h5=Δtn,h6=Δtn+1;
4.如权利要求1所述的一种保持强稳定性的变步长多步法时间离散算法,其特征在于:所述离散算法的初始步使用经典保强稳定的Runge-Kutta时间离散技术,即:
Q(1)=Qn+ΔtL(Qn)
初始步的时间步长可以选取正常计算的时间步长,其中n为整数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510132095.7A CN104699909B (zh) | 2015-03-25 | 2015-03-25 | 一种保持强稳定性的变步长多步法时间离散方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510132095.7A CN104699909B (zh) | 2015-03-25 | 2015-03-25 | 一种保持强稳定性的变步长多步法时间离散方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104699909A true CN104699909A (zh) | 2015-06-10 |
CN104699909B CN104699909B (zh) | 2018-11-02 |
Family
ID=53347025
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510132095.7A Expired - Fee Related CN104699909B (zh) | 2015-03-25 | 2015-03-25 | 一种保持强稳定性的变步长多步法时间离散方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104699909B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107992696A (zh) * | 2017-12-13 | 2018-05-04 | 电子科技大学 | 一种复杂色散媒质中的改进的指数时间积分构造方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20050061908A1 (en) * | 2000-05-31 | 2005-03-24 | Kevin Kremeyer | Shock wave modification method and system |
US20070011647A1 (en) * | 2003-04-06 | 2007-01-11 | Daniel Abrams | Optimized photomasks for photolithography |
CN104268322A (zh) * | 2014-06-27 | 2015-01-07 | 北京航空航天大学 | Weno差分方法的一种边界处理技术 |
-
2015
- 2015-03-25 CN CN201510132095.7A patent/CN104699909B/zh not_active Expired - Fee Related
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20050061908A1 (en) * | 2000-05-31 | 2005-03-24 | Kevin Kremeyer | Shock wave modification method and system |
US20070011647A1 (en) * | 2003-04-06 | 2007-01-11 | Daniel Abrams | Optimized photomasks for photolithography |
CN104268322A (zh) * | 2014-06-27 | 2015-01-07 | 北京航空航天大学 | Weno差分方法的一种边界处理技术 |
Non-Patent Citations (3)
Title |
---|
JUN ZHU,ETC: "Finite volume Hermite WENO schemes for solving the Hamilton-Jacobi equations II:Unstructured meshes", 《COMPUTERS AND MATHEMATICS WITH APPLICATIONS》 * |
JUNZHU,ETC: "Finite Volume Hermite WENO schemes for solving the Hamilton-Jacobi Equation", 《COMMUNICATION COMPUTER PHYSICS》 * |
QIUJIANXIAN,ETC: "Hermite WENO schemes for Hamilton-Jacobi equations on unstructured meshes", 《JOURNAL OF COMPUTATIONAL PHYSICS》 * |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107992696A (zh) * | 2017-12-13 | 2018-05-04 | 电子科技大学 | 一种复杂色散媒质中的改进的指数时间积分构造方法 |
CN107992696B (zh) * | 2017-12-13 | 2020-12-29 | 电子科技大学 | 一种复杂色散媒质中的改进的指数时间积分构造方法 |
Also Published As
Publication number | Publication date |
---|---|
CN104699909B (zh) | 2018-11-02 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Greve et al. | Initial results of the SeaRISE numerical experiments with the models SICOPOLIS and IcIES for the Greenland ice sheet | |
Zalesny et al. | Numerical model of the circulation of the Black Sea and the Sea of Azov | |
CN108763683B (zh) | 一种三角函数框架下新weno格式构造方法 | |
Abedian et al. | A high-order weighted essentially non-oscillatory (WENO) finite difference scheme for nonlinear degenerate parabolic equations | |
CN110457806A (zh) | 基于交错网格的中心五阶weno格式的全流场模拟方法 | |
Castonguay et al. | Application of high-order energy stable flux reconstruction schemes to the Euler equations | |
CN104504189B (zh) | 随机激励下大规模结构设计方法 | |
CN102495932A (zh) | 一种基于响应面建模和改进粒子群算法的有限元模型修正方法 | |
Zhang et al. | A class of DG/FV hybrid schemes for conservation law IV: 2D viscous flows and implicit algorithm for steady cases | |
CN107066720A (zh) | 一种基于piv技术的可压缩流体压力场的计算方法和装置 | |
CN102141064A (zh) | 空间过滤法建立湍流模型的构建方法 | |
Katta et al. | High-order finite-volume transport on the cubed sphere: Comparison between 1D and 2D reconstruction schemes | |
CN104091065A (zh) | 一种求解浅水问题模拟间断水流数值的方法 | |
CN107526105A (zh) | 一种波场模拟交错网格有限差分方法 | |
CN107657075A (zh) | 模拟地下水介质交界面处达西速度的区域分解有限元法 | |
CN102163263B (zh) | 风机叶片振动位移及其威布尔分布拟合方法 | |
CN104699909A (zh) | 一种保持强稳定性的变步长多步法时间离散算法 | |
CN103970610A (zh) | 一种供水管网中节点流量的监控方法 | |
CN105335552A (zh) | 不可伸展带状物体的几何性质描述模型及动力学模拟方法 | |
Turner et al. | MPAS-Seaice: a new variable resolution sea-ice model | |
Mai-Duy et al. | A compact 9 point stencil based on integrated RBFs for the convection–diffusion equation | |
CN105973319A (zh) | 一种控制棒驱动机构排污系统水力特性计算方法 | |
CN106021711A (zh) | 一种面向密集频率结构振动特征值的随机摄动方法 | |
CN106066912A (zh) | 一种多分区结构化网格的生成方法 | |
Zhao et al. | The discontinuous Petrov—Galerkin method for one-dimensional compressible Euler equations in the Lagrangian coordinate |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C53 | Correction of patent of invention or patent application | ||
CB03 | Change of inventor or designer information |
Inventor after: Qiu Jianxian Inventor after: Cai Xiaofeng Inventor before: Qiu Jianxian Inventor before: Zhang Fengyan Inventor before: Sun Chunpeng |
|
COR | Change of bibliographic data |
Free format text: CORRECT: INVENTOR; FROM: QIU JIANXIAN ZHANG FENGYAN SUN CHUNPENG TO: QIU JIANXIAN CAI XIAOFENG |
|
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20181102 Termination date: 20190325 |