CN107273625A - 一种三维不可压缩非定常n‑s方程有限元数值求解方法 - Google Patents

一种三维不可压缩非定常n‑s方程有限元数值求解方法 Download PDF

Info

Publication number
CN107273625A
CN107273625A CN201710481774.4A CN201710481774A CN107273625A CN 107273625 A CN107273625 A CN 107273625A CN 201710481774 A CN201710481774 A CN 201710481774A CN 107273625 A CN107273625 A CN 107273625A
Authority
CN
China
Prior art keywords
mrow
msub
mfrac
msubsup
moment
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.)
Pending
Application number
CN201710481774.4A
Other languages
English (en)
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.)
Southwest University of Science and Technology
Original Assignee
Southwest 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 Southwest University of Science and Technology filed Critical Southwest University of Science and Technology
Priority to CN201710481774.4A priority Critical patent/CN107273625A/zh
Publication of CN107273625A publication Critical patent/CN107273625A/zh
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

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)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

本发明公开了一种三维不可压缩非定常N‑S方程有限元数值求解方法,首先应用数学分裂的方法将N‑S方程分裂成扩散方程、对流方程以及压力修正方程,再根据n时刻值的初值和边界条件,由扩散方程计算得到n+1时刻的第一速度过渡值,然后由对流方程算出n+1时刻第二速度过渡值,根据压力泊松方程计算得到n+1时刻的压力值,根据速度修正方程计算得到n+1时刻的速度值,采用支反力公式计算流体对边界的作用力,最后判断当前时间步是否达到模拟时长,若达到则停止计算,否则进入下一步时间循环计算直到达到模拟时长。本发明中速度‑压力插值函数无需满足LBB条件,有效避免了有限元修正权函数的困难,同时计算效率大大提高。

Description

一种三维不可压缩非定常N-S方程有限元数值求解方法
技术领域
本发明属于流体力学技术领域,具体涉及一种三维不可压缩非定常N-S方程有限元数值求解方法的设计。
背景技术
流体流动是工程实际中普遍存在的现象,比如在航空、机械、建筑、动力、水力、化工、能源、环境、生物等工程领域,存在着大量的与流体流动相关的问题。在某种意义上讲,正是在流体力学问题的研究不断取得新成果的前提下,才促进了这些工程技术领域的大力发展。
湍流运动又是流体运动的常见现象,Navier-Stokes(N-S)方程是描述湍流运动的基本控制方程,它能够从本质上描述流体运动的规律,如大气运动、海洋流动、轴承润滑、透平机械内部流动等。但是标准Calerkin有限元法用于非定常不可压N-S方程求解时,会出现数值振荡。主要原因有两点:一是由于压力变量不出现在连续方程中,离散过程中速度场和压力场有限元插值函数的组合不恰当,使得LBB条件下不能满足,从而引起压力场的数值振荡;二是动量方程中存在非线性的对流项,当雷诺数较高对流占优时,标准Galerkin有限元求解法将导致数值振荡。
为避免速度-压力插值函数不满足LBB条件而导致的压力场的数值振荡,Chorin提出了投影法,速度-压力的插值函数无需满足LBB条件并且能够高效地计算非定常问题。在针对对流占优导致的数值解振荡,学者们发展了多种稳定化有限元解决该问题。目前应用较为广泛的有SUPG法、Taylor-Calerkin法、特征-Galerkin法等。Wang等将Taylor展开引入到特征线-Galerkin法,结合算子分裂法的优点提出特征线算子分裂(CBOS)有限元。对流项借鉴CBS算法,避免了P-G法等其他有限元法修正权函数的困难。但是由于压力插值函数比速度插值函数低一阶,使得不能随意取速度-压力插值函数,对流方程采用显示求解,同时时间步长要求较为苛刻,计算成本高。
发明内容
本发明的目的是针对现有技术中的上述不足,提出一种三维不可压缩非定常N-S方程有限元数值求解方法,通过引入一个无需满足散度为零约束条件的中间变量对压力和速度解耦求解,使得计算速度更快;此时速度与压力插值函数可以采用同阶插值,避免速度与压力插值函数不满足LBB条件而导致的压力场的数值振荡。
本发明的技术方案为:一种三维不可压缩非定常N-S方程有限元数值求解方法,包括以下步骤:
S1、建立工程对象数字模型,并对工程对象数字模型进行计算网格划分;
S2、设置模拟计算时的材料参数、边界条件、模拟时长及时间步长;
S3、将N-S方程应用于计算网格划分后的工程对象数字模型,应用数学分裂的方法将N-S方程:
▽·ui=0 (1)
分裂成扩散方程:
对流方程:
压力修正方程:
式(1)-(2)中,ui表示流体在x,y,z三个方向上的速度(u,v,w),下标i=1,2,3,p为压力,t为时间,Re为雷诺数,fi表示流体在x,y,z三个方向上的质量力(fx,fy,fz),▽为哈密顿算子,Δ为拉普拉斯算子;式(3)中,表示n+1时刻的第一速度过渡值;式(4)中,表示n+1时刻的第二速度过渡值;式(5)-(6)中,表示n+1时刻的速度值,pn+1表示n+1时刻的压力值。
S4、根据材料参数、边界条件以及扩散方程计算得到n+1时刻的第一速度过渡值;
S5、根据第一速度过渡值和对流方程计算得到n+1时刻的第二速度过渡值;
S6、根据压力修正方程得到压力泊松方程,结合第二速度过渡值计算得到n+1时刻的压力值;
S7、根据压力修正方程得到速度修正方程,结合第二速度过渡值及n+1时刻的压力值计算得到n+1时刻的速度值;
S8、采用支反力公式计算流体对边界的作用力;
S9、判断当前时间步是否达到模拟时长,若达到则停止计算,否则返回步骤S2。
本发明的有益效果是:本发明首次将N-S方程应用于三维模型,同时借用投影法特征线算子分裂有限元的思路,将三维不可压缩非定常N-S方程的压力、速度进行解耦计算,速度与压力插值函数可以采用同阶插值,使得速度-压力插值函数无需满足LBB条件,有效的避免了P-G法等有限元修正权函数的困难,同时由于压力和速度的解耦计算,使得计算效率大大提高。
进一步地,步骤S4具体为:
对于扩散方程(3)采用三步有限元格式求解n+1时刻的第一速度过渡值
式(7)-(9)中,Δt为时间步长,表示扩散方程在n+Δt/3时刻的解,示扩散方程在n+Δt/2时刻的解,xj表示方向,下标j=1,2,3,分别对应xj取x方向、y方向、z方向。
上述进一步方案的有益效果为:扩散方程采用三步有限元格式,使得由n时刻出发计算n+1时刻函数值时具有三阶精度,提高了整体计算精度。
进一步地,步骤S5具体为:
对于对流方程(4)采用相容方程替代偏微分方程,得到:
对对流方程的相容方程(10)采用多步有限元格式,得到:
式(10)-(11)中,h为每个整体时间步内的分步数,下标i=1,2,3,j=1,2,3,k=1,2,3,子时间步长若l=1则若l=h则将第一速度过渡值代入式(10)-(11),计算得到n+1时刻的第二速度过渡值
上述进一步方案的有益效果为:对流方程采用多步格式计算,降低了对整体时间步长的要求,提高了算法的稳定性。
进一步地,步骤S8具体为:
采用支反力公式计算流体对边界的作用力:
式(14)中,n表示n时刻,Γ为边界面积分项,Ω为体积分项,δui表示速度虚位移,δp表示压力虚位移。
上述进一步方案的有益效果为:借鉴固体力学支反力概念来计算流体对固体边界的作用力,并将该算法首次应用与三维模型中。传统求解支反力的方法是根据积分公式出发,其中涉及角度的正、余弦值,给计算复杂模型中流体对固体边界作用力带来了一定的困难;而借鉴固体力学支反力概念来计算流体对固体边界的作用力是根据N-S方程的扩散项的弱形式出发,使得求解复杂模型中流体对固体边界的作用力变得简单而不受模型形状的制约。
进一步地,步骤S2中的边界条件包括速度边界条件和出口对流边界条件,出口对流边界条件表示为:
式中uc为入口来流的最大速度。
上述进一步方案的有益效果为:引入了出口边界条件来加快计算速度,同时使得自由出口边界条件的稳定度大大提高。
附图说明
图1所示为本发明实施例提供的一种三维不可压缩非定常N-S方程有限元数值求解方法流程图。
图2所示为本发明实施例提供的方腔流计算模型示意图。
图3所示为本发明实施例提供的雷诺数为1000时三维顶板驱动方腔流等值面图。
图4所示为本发明实施例提供的三维方腔内中垂面上x速度分量与现有数值模拟对比示意图。
具体实施方式
现在将参考附图来详细描述本发明的示例性实施方式。应当理解,附图中示出和描述的实施方式仅仅是示例性的,意在阐释本发明的原理和精神,而并非限制本发明的范围。
本发明实施例提供了一种三维不可压缩非定常N-S方程有限元数值求解方法,如图1所示,包括以下步骤S1-S9:
S1、建立工程对象数字模型,并对工程对象数字模型进行计算网格划分。本发明实施例中,工程对象数字模型是采用数值分析通用图形用户界面进行创建的,如图2所示,数字模型为长宽高都为1的方腔,采用无量纲形式,模型中的网格被划分为六面体八节点单元。
S2、设置模拟时的相关参数及条件,具体包括:
(1)材料参数。
本发明实施例中,材料参数包括材料密度、雷诺数Re、质量力fi,参数具体数值设置如下表:
(2)边界条件。
本发明实施例中,边界条件包括速度边界条件和出口对流边界条件,出口对流边界条件表示为:
式中uc为入口来流的最大速度。
速度边界条件设置为:顶板设置无因次速度u=1,v=0,w=0,四周壁为不可滑移边界条件(u=0,v=0,w=0)。在设置边界条件的同时应当设置相应的压力参考点。
(3)模拟时长及时间步长。
本发明实施例中,模拟时长设置为Tmax=50,时间步长设置为Δt=0.01。
S3、将N-S方程应用于计算网格划分后的工程对象数字模型,应用数学分裂的方法将N-S方程:
▽·ui=0 (1)
分裂成扩散方程:
对流方程:
压力修正方程:
式(1)-(2)中,ui表示流体在x,y,z三个方向上的速度(u,v,w),下标i=1,2,3,p为压力,t为时间,Re为雷诺数,fi表示流体在x,y,z三个方向上的质量力(fx,fy,fz),▽为哈密顿算子,Δ为拉普拉斯算子;式(3)中,表示n+1时刻的第一速度过渡值;式(4)中,表示n+1时刻的第二速度过渡值;式(5)-(6)中,表示n+1时刻的速度值,pn+1表示n+1时刻的压力值。
S4、根据材料参数、边界条件以及扩散方程计算得到n+1时刻的第一速度过渡值。
本发明实施例中,对于扩散方程(3)采用三步有限元格式求解:
式(7)-(9)中即在每一整体时间步内将扩散方程分为3步进行求解,其中Δt为时间步长,表示扩散方程在n+Δt/3时刻的解,示扩散方程在n+Δt/2时刻的解,xj表示方向,下标j=1,2,3,分别对应xj取x方向、y方向、z方向。
代入n时刻值的速度初值和步骤S2设定的边界条件即可求得n+1时刻的第一速度过渡值
S5、根据第一速度过渡值和对流方程计算得到n+1时刻的第二速度过渡值。
本发明实施例中,对于对流方程(4)采用相容方程替代偏微分方程,得到:
对对流方程的相容方程(10)采用多步有限元格式,得到:
式(10)-(11)中,在每一个整体时间步内将对流方程分为h步进行求解,h为每个整体时间步内的分步数,下标i=1,2,3,j=1,2,3,k=1,2,3,子时间步长若l=1则若l=h则将第一速度过渡值代入式(10)-(11),计算得到n+1时刻的第二速度过渡值
S6、对式(5)两边同时取散度并联立式(6),得到压力泊松方程:
将第二速度过渡值代入式(12),计算得到n+1时刻的压力值pn+1
S7、整理式(5)得到速度修正方程:
将第二速度过渡值及n+1时刻的压力值pn+1代入式(13),计算得到n+1时刻的速度值
S8、采用支反力公式计算流体对边界的作用力:
式(14)中,n表示n时刻,Γ为边界面积分项,Ω为体积分项,δui表示速度虚位移,δp表示压力虚位移。
S9、判断当前时间步是否达到模拟时长Tmax,若达到则停止计算,否则返回步骤S2进入下一步时间循环计算,直到达到模拟时长Tmax
如图3所示,从左至右分别为三维顶板驱动方腔流雷诺数为1000时,速度值为0.18的等值面图、压力等值面图、涡量等值面图。根据速度值为0.18的等值面图,由于侧面壁的阻滞作用,流体在侧面壁附近堆积并发生翻卷,从而导致方腔流角涡的形成;从压力等值面图和涡量等值面图也可以看出由于侧壁的影响而存在的三维效应。
如图4所示为三维方腔内中垂面上x速度分量与现有数值模拟对比示意图。通过将本发明算法计算得出的三维方腔内中垂面上x速度分量与KL.Wong等学者对中、低雷诺数下的方腔流进行的数值模拟结果进行对比可知,二者曲线基本吻合,说明本发明有较高的计算精度及稳定性。
本领域的普通技术人员将会意识到,这里所述的实施例是为了帮助读者理解本发明的原理,应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合,这些变形和组合仍然在本发明的保护范围内。

Claims (9)

1.一种三维不可压缩非定常N-S方程有限元数值求解方法,其特征在于,包括以下步骤:
S1、建立工程对象数字模型,并对工程对象数字模型进行计算网格划分;
S2、设置模拟计算时的材料参数、边界条件、模拟时长及时间步长;
S3、将N-S方程应用于计算网格划分后的工程对象数字模型,并应用数学分裂的方法将所述N-S方程分裂成扩散方程、对流方程以及压力修正方程;
S4、根据所述材料参数、边界条件以及扩散方程计算得到n+1时刻的第一速度过渡值;
S5、根据所述第一速度过渡值和对流方程计算得到n+1时刻的第二速度过渡值;
S6、根据所述压力修正方程得到压力泊松方程,结合所述第二速度过渡值计算得到n+1时刻的压力值;
S7、根据所述压力修正方程得到速度修正方程,结合所述第二速度过渡值及n+1时刻的压力值计算得到n+1时刻的速度值;
S8、采用支反力公式计算流体对边界的作用力;
S9、判断当前时间步是否达到模拟时长,若达到则停止计算,否则返回步骤S2。
2.根据权利要求1所述的三维不可压缩非定常N-S方程有限元数值求解方法,其特征在于,所述步骤S2中的材料参数包括材料密度、雷诺数、质量力。
3.根据权利要求2所述的三维不可压缩非定常N-S方程有限元数值求解方法,其特征在于,所述步骤S3具体为:
应用数学分裂的方法将N-S方程:
<mrow> <mo>&amp;dtri;</mo> <mo>&amp;CenterDot;</mo> <msub> <mi>u</mi> <mi>i</mi> </msub> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>u</mi> <mi>i</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mrow> <mo>(</mo> <msub> <mi>u</mi> <mi>i</mi> </msub> <mo>&amp;CenterDot;</mo> <mo>&amp;dtri;</mo> <mo>)</mo> </mrow> <msub> <mi>u</mi> <mi>i</mi> </msub> <mo>=</mo> <mo>-</mo> <mo>&amp;dtri;</mo> <mo>&amp;CenterDot;</mo> <mi>p</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mi>Re</mi> </mfrac> <mi>&amp;Delta;</mi> <mo>&amp;CenterDot;</mo> <msub> <mi>u</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>f</mi> <mi>i</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
分裂成扩散方程:
<mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>-</mo> <mfrac> <mn>1</mn> <mi>Re</mi> </mfrac> <mi>&amp;Delta;</mi> <mo>&amp;CenterDot;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> <mo>=</mo> <msub> <mi>f</mi> <mi>i</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
对流方程:
<mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>2</mn> </msub> </mrow> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>2</mn> </msub> </mrow> </msubsup> <mo>&amp;CenterDot;</mo> <mo>&amp;dtri;</mo> <mo>&amp;CenterDot;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>2</mn> </msub> </mrow> </msubsup> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
压力修正方程:
<mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mo>-</mo> <mo>&amp;dtri;</mo> <mo>&amp;CenterDot;</mo> <msup> <mi>p</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <mo>&amp;dtri;</mo> <mo>&amp;CenterDot;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
式(1)-(2)中,ui表示流体在x,y,z三个方向上的速度(u,v,w),下标i=1,2,3,p为压力,t为时间,Re为雷诺数,fi表示流体在x,y,z三个方向上的质量力(fx,fy,fz),为哈密顿算子,Δ为拉普拉斯算子;式(3)中,表示n+1时刻的第一速度过渡值;式(4)中,表示n+1时刻的第二速度过渡值;式(5)-(6)中,表示n+1时刻的速度值,pn+1表示n+1时刻的压力值。
4.根据权利要求3所述的三维不可压缩非定常N-S方程有限元数值求解方法,其特征在于,所述步骤S4具体为:
对于扩散方程(3)采用三步有限元格式求解n+1时刻的第一速度过渡值
<mrow> <mfrac> <mrow> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mfrac> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> <mn>3</mn> </mfrac> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mi>n</mi> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>/</mo> <mn>3</mn> </mrow> </mfrac> <mo>=</mo> <mfrac> <mn>1</mn> <mi>Re</mi> </mfrac> <mfrac> <mo>&amp;part;</mo> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> </mfrac> <mrow> <mo>(</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mi>n</mi> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <mfrac> <mrow> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mfrac> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> <mn>2</mn> </mfrac> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mi>n</mi> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>/</mo> <mn>2</mn> </mrow> </mfrac> <mo>=</mo> <mfrac> <mn>1</mn> <mi>Re</mi> </mfrac> <mfrac> <mo>&amp;part;</mo> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> </mfrac> <mrow> <mo>(</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mfrac> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> <mn>2</mn> </mfrac> </mrow> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <mfrac> <mrow> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mi>n</mi> </msubsup> </mrow> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mn>1</mn> <mi>Re</mi> </mfrac> <mfrac> <mo>&amp;part;</mo> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> </mfrac> <mrow> <mo>(</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mfrac> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> <mn>2</mn> </mfrac> </mrow> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>
式(7)-(9)中,Δt为时间步长,表示扩散方程在n+Δt/3时刻的解,示扩散方程在n+Δt/2时刻的解,xj表示方向,下标j=1,2,3,分别对应xj取x方向、y方向、z方向。
5.根据权利要求4所述的三维不可压缩非定常N-S方程有限元数值求解方法,其特征在于,所述步骤S5具体为:
对于对流方程(4)采用相容方程替代偏微分方程,得到:
<mrow> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>2</mn> </msub> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> <mo>=</mo> <mo>-</mo> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>&amp;CenterDot;</mo> <msubsup> <mi>u</mi> <mi>j</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <msup> <mi>&amp;Delta;t</mi> <mn>2</mn> </msup> </mrow> <mn>2</mn> </mfrac> <msubsup> <mi>u</mi> <mi>k</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> <mfrac> <mo>&amp;part;</mo> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>k</mi> </msub> </mrow> </mfrac> <mrow> <mo>(</mo> <msubsup> <mi>u</mi> <mi>j</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>
对对流方程的相容方程(10)采用多步有限元格式,得到:
<mrow> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mi>h</mi> </mfrac> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mfrac> <mrow> <mi>l</mi> <mo>-</mo> <mn>1</mn> </mrow> <mi>h</mi> </mfrac> </mrow> </msubsup> <mo>=</mo> <mo>-</mo> <mi>&amp;delta;</mi> <mi>t</mi> <mo>&amp;CenterDot;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mfrac> <mrow> <mi>l</mi> <mo>-</mo> <mn>1</mn> </mrow> <mi>h</mi> </mfrac> </mrow> </msubsup> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mfrac> <mrow> <mi>l</mi> <mo>-</mo> <mn>1</mn> </mrow> <mi>h</mi> </mfrac> </mrow> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <msup> <mi>&amp;delta;t</mi> <mn>2</mn> </msup> </mrow> <mn>2</mn> </mfrac> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mfrac> <mrow> <mi>l</mi> <mo>-</mo> <mn>1</mn> </mrow> <mi>h</mi> </mfrac> </mrow> </msubsup> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mfrac> <mrow> <mi>l</mi> <mo>-</mo> <mn>1</mn> </mrow> <mi>h</mi> </mfrac> </mrow> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>j</mi> </msub> </mrow> </mfrac> <msubsup> <mi>u</mi> <mi>k</mi> <mrow> <mi>n</mi> <mo>+</mo> <mfrac> <mrow> <mi>l</mi> <mo>-</mo> <mn>1</mn> </mrow> <mi>h</mi> </mfrac> </mrow> </msubsup> <mfrac> <mrow> <msub> <mi>&amp;delta;u</mi> <mi>i</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>x</mi> <mi>k</mi> </msub> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>
式(10)-(11)中,h为每个整体时间步内的分步数,下标i=1,2,3,j=1,2,3,k=1,2,3,子时间步长l=1,2,...,h,若l=1则若l=h则将第一速度过渡值代入式(10)-(11),计算得到n+1时刻的第二速度过渡值
6.根据权利要求5所述的三维不可压缩非定常N-S方程有限元数值求解方法,其特征在于,所述步骤S6具体为:
对式(5)两边同时取散度并联立式(6),得到压力泊松方程:
<mrow> <mi>&amp;Delta;</mi> <mo>&amp;CenterDot;</mo> <msup> <mi>p</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>&amp;Delta;</mi> <mi>t</mi> </mrow> </mfrac> <mo>&amp;dtri;</mo> <mo>&amp;CenterDot;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>2</mn> </msub> </mrow> </msubsup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow>
将第二速度过渡值代入式(12),计算得到n+1时刻的压力值pn+1
7.根据权利要求6所述的三维不可压缩非定常N-S方程有限元数值求解方法,其特征在于,所述步骤S7具体为:
整理式(5)得到速度修正方程:
<mrow> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>2</mn> </msub> </mrow> </msubsup> <mo>=</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mi>&amp;Delta;</mi> </mfrac> <mo>&amp;dtri;</mo> <mo>&amp;CenterDot;</mo> <msup> <mi>p</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow> 2
将第二速度过渡值及n+1时刻的压力值pn+1代入式(13),计算得到n+1时刻的速度值
8.根据权利要求3所述的三维不可压缩非定常N-S方程有限元数值求解方法,其特征在于,所述步骤S8具体为:
采用支反力公式计算流体对边界的作用力:
<mrow> <mtable> <mtr> <mtd> <mrow> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <msub> <mo>&amp;Integral;</mo> <mi>&amp;Gamma;</mi> </msub> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <mi>n</mi> </mrow> </mfrac> <msub> <mi>&amp;delta;u</mi> <mi>i</mi> </msub> <mo>+</mo> <mi>p</mi> <mo>&amp;CenterDot;</mo> <mi>n</mi> <mi>d</mi> <mi>&amp;Gamma;</mi> <mo>=</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <msub> <mo>&amp;Integral;</mo> <mi>&amp;Omega;</mi> </msub> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <msub> <mi>&amp;delta;u</mi> <mi>i</mi> </msub> <mi>d</mi> <mi>&amp;Omega;</mi> <mo>-</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <msub> <mo>&amp;Integral;</mo> <mi>&amp;Omega;</mi> </msub> <mrow> <mo>(</mo> <mi>p</mi> <mo>&amp;dtri;</mo> <mo>&amp;CenterDot;</mo> <msub> <mi>&amp;delta;u</mi> <mi>i</mi> </msub> <mo>+</mo> <mo>&amp;dtri;</mo> <mo>&amp;CenterDot;</mo> <msubsup> <mi>&amp;delta;u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> <mi>&amp;delta;</mi> <mi>p</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>&amp;Omega;</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>-</mo> <mfrac> <mn>1</mn> <mi>Re</mi> </mfrac> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <msub> <mo>&amp;Integral;</mo> <mi>&amp;Omega;</mi> </msub> <mo>&amp;dtri;</mo> <mo>&amp;CenterDot;</mo> <msubsup> <mi>&amp;delta;u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> <mo>&amp;dtri;</mo> <mo>&amp;CenterDot;</mo> <msub> <mi>&amp;delta;u</mi> <mi>i</mi> </msub> <mi>d</mi> <mi>&amp;Omega;</mi> <mo>-</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <msub> <mo>&amp;Integral;</mo> <mi>&amp;Omega;</mi> </msub> <msub> <mi>f</mi> <mi>i</mi> </msub> <msub> <mi>&amp;delta;u</mi> <mi>i</mi> </msub> <mi>d</mi> <mi>&amp;Omega;</mi> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>
式(14)中,n表示n时刻,Γ为边界面积分项,Ω为体积分项,δui表示速度虚位移,δp表示压力虚位移。
9.根据权利要求3所述的三维不可压缩非定常N-S方程有限元数值求解方法,其特征在于,所述步骤S2中的边界条件包括速度边界条件和出口对流边界条件,所述出口对流边界条件表示为:
<mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> </mrow> <mrow> <mo>&amp;part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <msub> <mi>u</mi> <mi>c</mi> </msub> <mo>&amp;dtri;</mo> <mo>&amp;CenterDot;</mo> <msubsup> <mi>u</mi> <mi>i</mi> <mrow> <mi>n</mi> <mo>+</mo> <msub> <mi>&amp;theta;</mi> <mn>1</mn> </msub> </mrow> </msubsup> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> </mrow>
式中uc为入口来流的最大速度。
CN201710481774.4A 2017-06-22 2017-06-22 一种三维不可压缩非定常n‑s方程有限元数值求解方法 Pending CN107273625A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710481774.4A CN107273625A (zh) 2017-06-22 2017-06-22 一种三维不可压缩非定常n‑s方程有限元数值求解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710481774.4A CN107273625A (zh) 2017-06-22 2017-06-22 一种三维不可压缩非定常n‑s方程有限元数值求解方法

Publications (1)

Publication Number Publication Date
CN107273625A true CN107273625A (zh) 2017-10-20

Family

ID=60068103

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710481774.4A Pending CN107273625A (zh) 2017-06-22 2017-06-22 一种三维不可压缩非定常n‑s方程有限元数值求解方法

Country Status (1)

Country Link
CN (1) CN107273625A (zh)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108460188A (zh) * 2018-02-05 2018-08-28 电子科技大学 一种应用于pic静电模型的电荷分配有限元fem求解算法
CN108509741A (zh) * 2018-04-10 2018-09-07 西南科技大学 一种泥石流方程的有限元数值求解方法
CN108962392A (zh) * 2018-05-29 2018-12-07 杭州晟视科技有限公司 一种速度场的修正方法、装置及存储介质
CN108984874A (zh) * 2018-07-02 2018-12-11 中国航空发动机研究院 获取不可压缩流动的流场的数值模拟方法
CN109726433A (zh) * 2018-11-30 2019-05-07 电子科技大学 基于曲面边界条件的三维无粘低速绕流的数值模拟方法
CN113505551A (zh) * 2021-09-09 2021-10-15 中国空气动力研究与发展中心计算空气动力研究所 来流变化诱发非定常的模拟方法、系统、存储介质和终端

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010210281A (ja) * 2009-03-06 2010-09-24 Nikon Corp 流れ算出方法、流れ算出装置およびプログラム
CN102339326A (zh) * 2010-07-16 2012-02-01 中国石油化工股份有限公司 一种分析模拟缝洞型油藏流体流动的方法
CN104331579A (zh) * 2014-11-19 2015-02-04 中国石油大学(华东) 一种低渗透储层原油边界层的模拟方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010210281A (ja) * 2009-03-06 2010-09-24 Nikon Corp 流れ算出方法、流れ算出装置およびプログラム
CN102339326A (zh) * 2010-07-16 2012-02-01 中国石油化工股份有限公司 一种分析模拟缝洞型油藏流体流动的方法
CN104331579A (zh) * 2014-11-19 2015-02-04 中国石油大学(华东) 一种低渗透储层原油边界层的模拟方法

Non-Patent Citations (10)

* Cited by examiner, † Cited by third party
Title
M.BREUER ET AL.: "Accurate computations of the laminar flow past a square cylinder based on two different methods: lattice-Boltzmann and finite-volume", 《INTERNATIONAL JOURNAL OF HEAT AND FLUID FLOW》 *
WANG DAGUO ET AL.: "Characteristic-based operator-splitting finite element method for Navier-Stokes equations", 《SCIENCE CHINA-TECHNOLOGICAL SCIENCES》 *
冯立伟: "非定常不可压Navier-Stokes方程的有限元解法", 《万方数据库》 *
水庆象等: "N-S方程基于投影法的特征线算子分裂有限元求解", 《力学学报》 *
江春波等: "求解不可压缩流动的分步有限元格式", 《清华大学学报(自然科学版)》 *
深圳大学复变函数与场论教研组: "《复变函数与场论简明教程》", 30 November 2012, 西安电子科技大学出版社 *
王伟等: "线性剪切来流对方柱绕流特性影响的数值分析", 《水利水电科技进展》 *
王晖等: "基于流体动力学的放空管应力分析", 《特种设备安全技术》 *
胡彬等: "特征线算子分裂有限元的圆柱绕流大涡模拟", 《水利水电科技进展》 *
齐金苑等: "《水利工程建设百科全书 勘测设计•施工技术•质量管理卷 第1册》", 31 August 2003, 当代中国音像出版 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108460188A (zh) * 2018-02-05 2018-08-28 电子科技大学 一种应用于pic静电模型的电荷分配有限元fem求解算法
CN108509741A (zh) * 2018-04-10 2018-09-07 西南科技大学 一种泥石流方程的有限元数值求解方法
CN108509741B (zh) * 2018-04-10 2022-07-29 西南科技大学 一种泥石流方程的有限元数值求解方法
CN108962392A (zh) * 2018-05-29 2018-12-07 杭州晟视科技有限公司 一种速度场的修正方法、装置及存储介质
CN108984874A (zh) * 2018-07-02 2018-12-11 中国航空发动机研究院 获取不可压缩流动的流场的数值模拟方法
CN109726433A (zh) * 2018-11-30 2019-05-07 电子科技大学 基于曲面边界条件的三维无粘低速绕流的数值模拟方法
CN113505551A (zh) * 2021-09-09 2021-10-15 中国空气动力研究与发展中心计算空气动力研究所 来流变化诱发非定常的模拟方法、系统、存储介质和终端

Similar Documents

Publication Publication Date Title
CN107273625A (zh) 一种三维不可压缩非定常n‑s方程有限元数值求解方法
Liang et al. Phase-field-based lattice Boltzmann modeling of large-density-ratio two-phase flows
Sun et al. Numerical investigation on the unsteady cavitation shedding dynamics over a hydrofoil in thermo-sensitive fluid
CN104143027B (zh) 一种基于sph算法的流体热运动仿真系统
Stuhmiller The influence of interfacial pressure forces on the character of two-phase flow model equations
Chen et al. Lattice Boltzmann model for incompressible axisymmetric flows
CN114662425B (zh) 一种水轮机启停工况流场仿真预测方法及系统
CN105975700B (zh) 一种模拟超声空泡动力学行为的数值方法
Rogers et al. On the accuracy of the pseudocompressibility method in solving the incompressible Navier-Stokes equations
CN101833605A (zh) 基于流体体积模型模具微细流道磨料流精密加工控制方法
Hua et al. Level set, phase-field, and immersed boundary methods for two-phase fluid flows
Li et al. Large eddy simulation of unsteady shedding behavior in cavitating flows with time-average validation
Rezvaya et al. Solving the hydrodynamical tasks using CFD programs
CN109101752A (zh) 一种复杂水工建筑物局部结构自振频率计算方法
Kidron et al. Robust Cartesian grid flow solver for high-Reynolds-number turbulent flow simulations
CN106773782A (zh) 一种气动伺服弹性混合建模方法
Li et al. A numerical investigation of the flow between rotating conical cylinders of two different configurations
CN108509741A (zh) 一种泥石流方程的有限元数值求解方法
Liu et al. Simulation of incompressible multiphase flows with complex geometry using etching multiblock method
Yu et al. A coupled lattice Boltzmann and particle level set method for free-surface flows
Dubrovskyi et al. The application of the extended cells method to simulate the flow of combustion gases in the LPRE chamber
Cui et al. 3-D CFD validation of an axisymmetric jet in cross-flow (NASA Langley Workshop Validation: Case 2)
Li et al. Simulating the frontal instability of lock-exchange density currents with dissipative particle dynamics
Mateescu et al. Analysis of unsteady three-dimensional confined flows with oscillating walls and variable inflow velocity
Chiu Computational fluid dynamics simulations of hydraulic energy absorber

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
RJ01 Rejection of invention patent application after publication

Application publication date: 20171020

RJ01 Rejection of invention patent application after publication