有限体积离散

通常,各种数值方法(包含有限体积方法)用于将数学模型转换成代数方程组。此转换涉及在空间和时间上对控制方程进行离散化。然后使用代数多重网格求解器求解生成的线性方程。

对于非稳态问题,要分析的物理时间间隔被细分成任意数量的子间隔(称作时间步)。

通用传输方程

将适当的本构关系引入守恒方程即可获得闭合的方程组。所有守恒方程都可根据通用传输方程编写。通过将通用传输方程对控制体积 V 进行积分并应用高斯的散度定理,可获得传输方程的以下积分形式:

1. EQUATION_DISPLAY
ddtV ρϕdVTransient Term+A ρvϕdaConvective Flux=AΓϕdaDiffusive Flux+VSϕdVSource Term
(879)

其中, 表示标量属性的传输, 为控制体积的表面积, 表示表面矢量。ϕAda例如,通过设置 等于ϕ 1,u,v,w,E,orH,Yi 并为扩散系数 和源项选择适当的值,可获得特殊形式的质量、动量、能量和组分守恒的偏微分方程。Γ

Eqn. (879) 具有四个截然不同的项:
  • 瞬态项,表示控制体积内流体属性 ϕ 的时间变化率。
  • 对流通量,表示对流导致的、流体属性 ϕ 通过整个控制体积边界的净减小率。
  • 扩散通量,对应于扩散导致的、流体属性 ϕ 通过整个控制体积边界的净增加率。
  • 源项,表示控制体积内流体属性 ϕ 的生成/破坏。
下图显示了描述通用传输方程离散化的两种多面体计算网格单元:

表面积分
使用求积近似在两个级别上计算 Eqn. (879) 中的表面积分:
  • 积分用网格单元面上同一位置处的变量值表示。

    Simcenter STAR-CCM+ 采用二阶中点规则,即,按照网格单元面中心处的值与网格单元面网格面积的乘积计算积分:

    AJϕdafJfϕaf

    其中,Jϕ 为流体属性 ϕ 的对流或扩散通量,af 为网格单元的面 f 的表面积矢量,f 为网格单元的所有网格单元面的总和。

  • 网格单元面中心处的值未知,并通过根据网格单元中心处的值进行插值来近似。这些插值格式(即离散格式)会在对流通量扩散通量下介绍。
体积积分
Eqn. (879) 中的源项需要对网格单元的体积进行积分。Simcenter STAR-CCM+ 通过计算网格单元中心处的源项的平均值与该网格单元体积的乘积得到近似体积积分。此近似为二阶精确。
2. EQUATION_DISPLAY
VSϕVSϕ0V0
(880)

要获得这些近似的二阶精度,网格单元面中心需要为面积加权中心并且网格单元中心需要为体积中心。

将该积分近似应用于 Eqn. (879) 会生成以下半离散传输方程:

3. EQUATION_DISPLAY
ddt(ρϕV)0+f [ρϕ(va )]f=f (Γϕa)f+(SϕV)0
(881)

对流通量

面处的离散化对流项可重新排列如下:

4. EQUATION_DISPLAY
(ϕρva )f=(m˙ϕ)f=m˙fϕf
(882)

其中,m˙f 为面处的质量流率。

Eqn. (882) 需要面处的流体属性 ϕf 的值。根据网格单元值计算流体属性面值 ϕf 所采用的方式对数值格式的稳定性和精度具有重大影响。

标准化变量图
标准化变量图 (NVD) 用于分析对流离散格式的有界性属性。

下图中的部分 (a) 显示了网格单元面 f 附近的三个网格单元,跨这三个网格单元的速度 vf 已知。节点变量值用 αDαCαU 标记,表示彼此相对的顺风位置、中间位置和迎风位置。

对流有界性标准在 NVD 图中显示为图的部分 (b):



面 f 附近的标准化变量 ξ(r,t) 定义为:

5. EQUATION_DISPLAY
ξ(r,t)=α(r,t)-αUαD-αU
(883)

标准化面值:

6. EQUATION_DISPLAY
ξf=αf-αUαD-αU
(884)

通过仅使用位于点 U、C 和 D 处的 α 节点值的任何差分格式计算,可编写为以下形式:

ξf=ξf(ξC)

其中:

7. EQUATION_DISPLAY
ξC=αC-αUαD-αU
(885)

为防止求解中的不实际振荡,αC(以及后续的 αf)的局部范围必须介于 αUαD 之间,这意味着:

8. EQUATION_DISPLAY
αUαCαD
(886)
9. EQUATION_DISPLAY
αDαCαU
(887)

如果域中每个点都满足此标准,则整个求解中不会存在不实际振荡。

对流差分格式的有界性标准可以显示在 NVD 图中,将 ξf 显示为 ξC 的函数,即上图的部分 (b) 中的阴影区域(包括线 ξf=ξC)。它还可以通过以下条件表示:

  • 0ξC1 时,有界区域位于线 ξf=ξC 之上以及 ξf=1 之下。
  • ξC<0ξC>1 时,ξf 等于 ξC

NVD 仅与对流传输有关。如果存在源或汇,则可能违反 Eqn. (886)Eqn. (887) 中给出的条件。当变量具有物理边界时,有界性标准的重要性尤其明确;例如,相体积分数不得变为负数或大于 1。

一阶迎风

一阶迎风差分格式 (FOU) 以阶跃函数的形式近似网格单元面值 ϕf。该近似为一阶精确。对流通量计算如下,具体取决于流向:

10. EQUATION_DISPLAY
(m˙ϕ)f={m˙fϕ0       for m˙f0m˙fϕ1        for m˙f<0
(888)

其中,正向质量通量正在离开网格单元,负向质量通量正在进入网格单元。此离散格式基于流体的传输属性:由于对流的作用,ϕ 仅向下游传输。如果流体的流线与网格线对齐,则生成的近似准确。但如果流线与网格线形成角度,则 UDS 不准确。这些误差称作人工或数值扩散。UDS 的优点在于没有条件限制,从而可稳定和帮助求解器实现稳健收敛。很常见的做法是使用一阶对流获得初始求解,然后切换到二阶对流以获得更精确且最终收敛的求解。

二阶迎风

对于二阶迎风 (SOU) 格式,对流通量计算如下:

11. EQUATION_DISPLAY
(m˙ϕ)f={m˙fϕf,0       for m˙f0m˙fϕf,1        for m˙f<0
(889)

其中,面值 ϕf,0 ϕf,1 在面的任一侧上从网格单元中心值进行线性插值:

ϕf,0=ϕ0+s0(ϕ)r,0ϕf,1=ϕ1+s1(ϕ)r,1

其中:

s0=xf-x0s1=xf-x1

网格单元 0 和 1 的梯度 (ϕ)r,0 (ϕ)r,1 根据 Eqn. (911) 计算为受限重构梯度。

中心差分

中心差分格式 (CDS) 通过在两个最邻近网格单元中心值之间进行线性插值,来近似网格单元面中心值。

对于中心差分格式,对流通量计算如下:

12. EQUATION_DISPLAY
(m˙ϕ)f=m˙f[fϕ0+(1-f)ϕ1]
(890)

其中,线性插值因子定义为:

f = V 1 V 0 + V 1

f 与网格延伸有关。对于均匀网格, f 的值为 0.5。

中心差分是正式二阶精确的。但是,在大多数稳态情况下,它容易出现色散误差和稳定性问题。色散误差使得难以离散化不允许出现超调的正定量(例如温度或湍动能)。

二阶迎风的中心差分的一个重要优点是,当用于离散速度(非正定量)时,它会保留湍动能。因此,它是大涡模拟 (LES) 中的一种非常有用的格式,而迎风格式会导致湍动能异常快速地衰减。

有界中心差分

对于 [232][258] 中介绍的有界中心差分格式,对流通量计算如下:

13. EQUATION_DISPLAY
(m˙ϕ)f={m˙ϕFOUfor ξ<0  or  1<ξm˙(σϕCD+(1-σ)ϕSOU)for 0ξ1
(891)
其中:
ϕFOU 通过一阶迎风插值获得的网格单元面中心值
ϕSOU 通过二阶迎风插值获得的网格单元面中心值
ϕCD 通过中心差分插值获得的网格单元面中心值
ξ 基于局部条件计算的标准化变量图 (NVD) 值

ξ 的平滑单调函数如下:

14. EQUATION_DISPLAY
σ=σ(ξ)
(892)

ξubfξ 时,它满足 σ(0)=0 σ(ξ)=1 ξubf 称作迎风混合因子,其值可确保格式在精度和稳定性之间的适当平衡。因此, ξubf 的值越小,精度越高,而值越大,格式的稳定性越高。

如果无法满足对流有界性标准(例如当 ξ<0  or  1<ξ 时),为了保持有界性,有界中心差分格式会转换成一阶迎风格式。而中心差分格式(正式二阶精确格式)不会如此。因此,有界中心差分格式可能比中心差分格式更加耗散,尤其是在网格较粗糙的情况下。

总之,有界中心差分格式可有效平衡精度(与二阶迎风格式相比明显提高)和稳定性(由于其有界性有保证)。建议将该格式用于复杂湍流的 LES。
混合 MUSCL 三阶/中心差分

该格式对于稳态和非稳态模拟都适用,并且具有一个模型参数 σMUSCL3,用于控制格式中的数值耗散。对于有界中心差分格式,在非平滑流区域中,此格式使用标准化变量图 (NVD) 值 ξ,通过切换到一阶格式来确保格式的有界性。

检测到平滑局部流条件时,该格式将构建为 MUSCL 三阶迎风和三阶中心差分重构格式之间的混合。

对流通量计算如下:

15. EQUATION_DISPLAY
(m˙ϕ)f={m˙ϕFOUforξ<0orξ>1m˙(σMUSCL3ϕMUSCL3+(1σMUSCL3)ϕCD3)for0ξ1
(893)

混合因子 σMUSCL3 由用户控制,且必须根据实际问题或模型确定。

在高速下,需要谨慎限制 MUSCL 三阶迎风重构值 ϕMUSCL3,才能不影响精度的正式阶数,同时防止发生伪振荡。高速限制基于加权基本无振荡 (WENO) 原理,并且仅应用于面值重构的二次部分。基于 WENO 的平均根据三个不同的模板计算。

在靠近强冲击的区域中,ϕMUSCL3 精度降为二阶。

需要根据基于马赫数的通量限制器在无限和 WENO 限制 ϕMUSCL3 之间切换。对于不可压缩模拟,ϕMUSCL3不受限的(通量限制器为 1),因此,只需使用 NVD 值即可确保格式的有界性。总之,与二阶和有界中心差分 (BCD) 格式相比,混合 MUSCL 三阶/CD 格式可提高(减小)耗散。同时,它很稳定(由于其有界性)且能够模拟从不可压缩流态到高速可压缩流态的稳态流和非稳态流。

混合二阶迎风/中心差分

对于稳态计算,二阶迎风是最精确的格式,而中心差分格式的动能保留属性则使它更适用于大涡模拟 (LES)。在此格式中,这两种方法都按以下方程混合:

16. EQUATION_DISPLAY
(m˙ϕ)f={m˙f{σHUϕf,0+(1-σHU)[fϕf,0+(1-f)ϕf,1]}    for m˙f0m˙f{σHUϕf,1+(1-σHU)[fϕf,0+(1-f)ϕf,1]}      for m˙f<0
(894)

其中, σ H U 为根据流态选择的适当混合因子(请参见 Eqn. (1450)Eqn. (1306))。对于二阶迎风格式,面值 ϕ f , 0 ϕ f , 1 是使用重构梯度从面的任一侧上的网格单元值线性插值得到的。

混合二阶迎风/有界 - 中心

以下混合二阶迎风/有界中心格式是复杂湍流的分离涡模拟 (DES) 和尺度解析混合 (SRH) 模拟的默认格式:

17. EQUATION_DISPLAY
(m˙ϕ)f=m˙(σHUϕSOU+(1-σHU)ϕBCD)
(895)

其中,在面 处,对于二阶迎风格式和有界中心差分格式, 的值分别为 和 。 f ϕ ϕ S O U ϕ B C D σ H U 为根据流态选择的适当混合因子(对于 DES,请参见 Eqn. (1450);对于 SRH,请参见 Eqn. (1306))。

扩散通量

Eqn. (879) 中通过网格单元的内部网格单元面的扩散通量离散化如下:

18. EQUATION_DISPLAY
Df=(Γϕa)f
(896)

其中, Γ 为面扩散率, ϕ 为流体属性 ϕ 的梯度, a 为表面积矢量。

要获得隐式包含网格单元值 ϕ0 ϕ 1 的内部面梯度的精确二阶表达式,需要使用以下分解:

19. EQUATION_DISPLAY
ϕf=(ϕ1-ϕ0)α+ϕ¯(ϕ¯ds)α
(897)

其中:

α=aadsds=x1-x0ϕ¯=(ϕ0+ϕ1)2

然后可将内部面处的扩散通量写为:

20. EQUATION_DISPLAY
Df=Γfϕfa=Γf[(ϕ1-ϕ0)αa+ϕ¯a-(ϕ¯ds)αa]
(898)

其中, Γ f 为网格单元值的谐波平均值。

在边界面上使用类似的分解:

21. EQUATION_DISPLAY
Df=Γfϕfa=Γf[(ϕf-ϕ0)αa+ϕ0a-(ϕ0ds)αa]
(899)

其中:

22. EQUATION_DISPLAY
α=aads
(900)
23. EQUATION_DISPLAY
ds = x f - x 0
(901)

Eqn. (898) 中的第二和第三项表示二阶梯度(或交叉扩散)贡献。它们对于保持非正交网格的精度至关重要。

上面给出的公式假设网格单元 0 和 1 的形心位于面的相对两侧。此外,假设其位置与面网格面积矢量指向网格单元 0 外部的约定一致。要防止产生非物理解,谨慎使用 a ds 之间的角度不超过 90 度的网格。

为确保使用的网格有效,Simcenter STAR-CCM+ 中提供了一个诊断工具,允许此角度(称为偏斜角)以度为单位计算并存储在相邻网格单元中。存储在每个网格单元中的值表示每个网格单元面的最大偏斜角。

由于可以采用进一步的诊断工具检查二阶梯度贡献是否失稳,因此可选择忽略这些项,在这种情况下,Eqn. (898) 将简化为:

24. EQUATION_DISPLAY
DfΓf[(ϕ1-ϕ0)αa]
(902)

与内部面类似,Eqn. (899) 中的第二和第三项表示边界二阶梯度(或交叉扩散)贡献。

与内部一样,由于有进一步的诊断工具,因此可以选择忽略这些边界二阶梯度,这样 Eqn. (899) 将简化为:

25. EQUATION_DISPLAY
DfΓ[(ϕf-ϕ0)αa]
(903)

瞬态项

对于瞬态模拟,时间为额外坐标。通常,除了空间离散化之外,数值计算过程还要求对此坐标进行离散化,即,将总时间间隔细分为多个时间步。需要在不同的时间级别下获取控制方程的求解,时间 t 处的求解需要之前时间级别的求解。时间积分格式需要区分用于积分的时间级别数,以及在哪个时间级别上对通量和源积分。

隐式时间积分
欧拉隐式格式(一种一阶临时格式)使用当前时间级别 n + 1 和上一时间级别 n 中的求解来近似 Eqn. (879) 中的瞬态项:
26. EQUATION_DISPLAY
ddt(ρϕV)0=(ρϕV)0n+1-(ρϕV)0nΔt
(904)

在向后微分公式 BDF2 中,非稳态项的基本二阶时间离散使用当前时间级别 n + 1 求解以及前两个时间级别 n n - 1 中的求解。

BDF2 为:

27. EQUATION_DISPLAY
t(ρχϕV)=(32(ρχϕV)n+12(ρχϕV)n+12(ρχϕV)n1)1Δt
(905)

在二阶时间模拟的第一个时间步中,由于只有两个时间级别可用,因此使用了一阶离散。

使用更多时间级别的二阶格式为 BDF 的线性组合。

用于四个时间级别的格式 BDF2Opt(4) 为:

28. EQUATION_DISPLAY
BDF2Opt(4)=12(BDF3+BDF2)
(906)

其中,BDF3 为:

t(ρχϕV)=(116(ρχϕV)n+13(ρχϕV)n+32(ρχϕV)n113(ρχϕV)n2)1Δt

该格式使用前三个时间级别以及当前时间级别。在时间级别可用时,它在第一步使用一阶离散,并在第二步使用二阶离散。

用于五个时间级别的格式 BDF2Opt(5) 为:

29. EQUATION_DISPLAY
BDF2Opt(5)=(112)BDF4(5222)BDF3+(12+5222)BDF2
(907)

其中,BDF4 为:

t(ρχϕV)=(2512(ρχϕV)n+14(ρχϕV)n+3(ρχϕV)n143(ρχϕV)n2+14(ρχϕV)n3)1Δt

在时间级别可用时,该格式使用前四个时间级别并且在前三步使用前三个离散。

上述公式在整个过程中使用恒定时间步 Δt。但是,Simcenter STAR-CCM+ 实际上使用可变时间步来适应间隔不均匀的时间级别。

梯度

除了变量值外,还需要在网格单元中心处和网格单元面中心处的变量梯度执行以下操作:

  • 构造网格单元面处的变量值
  • 计算扩散项的二阶梯度
  • 计算压力-速度耦合的压力梯度
  • 计算湍流模型的应变率和旋转速率
最小二乘法

网格单元-0 的梯度计算如下:

30. EQUATION_DISPLAY
ϕ = ( β ) LSQ_grad ( ϕ )
(908)

其中:

31. EQUATION_DISPLAY
LSQ_grad ( ϕ ) = f M 1 ( w f ds f ) δ ϕ f = f q f δ ϕ f
(909)

其中:

  • ds f = x c 1 x c 0 为两个网格单元中心之间的矢量。
  • δ ϕ f = ϕ c 1 ϕ c 0 为标量 的网格单元中心值之差。 ϕ
  • 最小二乘张量 定义为: M
    M = f w f d s f d s f
    重量:
    w f = | A f | | d s f | V
    其中:
    • | A f | 是面面积矢量的幅值。
    • | d s f | 为网格单元中心之间的长度。
    • V 为网格单元体积。

标量 计算如下: β

32. EQUATION_DISPLAY
β = min ( 1 , C max user max f , g ( x f 0 q g ) )
(910)

其中:

  • C max user 为用户指定的最大重构系数(默认值 1)。
  • max f , g 针对范围 1 到 内所有网格单元面的所有配对 。 N faces f , g
  • x f 0 = x f x c 0 为面-网格单元矢量。
  • q f Eqn. (909) 中定义。

梯度限制器

如果使用不受限重构梯度重构,重构的面值可能位于相邻网格单元(通过面连接)中的网格单元值范围之外。为此,Simcenter STAR-CCM+ 会查找相邻网格单元值的下限和上限,然后使用这些界限来限制重构梯度。

可使用以下梯度限制器:

  • Venkatakrishnan
  • Min-Mod
  • 已修正 Venkatakrishnan

通过位于任意面形心处的网格单元 0 值重构的面值 ( ϕf,0 ) 由以下公式给出:

33. EQUATION_DISPLAY
ϕ f , 0 = ϕ 0 + s 0 ( ϕ ) 0
(911)

其中, s 0 = x f - x 0 x f x 0 分别为面形心和网格单元形心, ϕ r , 0 为重构梯度。

每个网格单元 0 都需要受限重构梯度,这样重构的面值才不会超过相邻网格单元形心值的最大值和最小值,包括网格单元 0 中的值。比例因子 α 定义为受限值和不受限值的比率,即:

34. EQUATION_DISPLAY
( ϕ ) r , 0 = α ϕ
(912)

对于每个网格单元 0,这两个数量定义为:

ϕ0max=max (ϕ0,ϕneighbors)ϕ0min=min (ϕ0,ϕneighbors)

其中, ϕneighbors 表示与网格单元 0 具有公共面的每个相邻网格单元中的网格单元值。

这两个数量还可以定义为:

Δmax=ϕ0max-ϕ0Δmin=ϕ0min-ϕ0

对于网格单元 0 的每个面 f,Δf 定义为:

35. EQUATION_DISPLAY
Δf=ϕf,0-ϕ0=s0(ϕ)r,0u
(913)

现在定义为:

36. EQUATION_DISPLAY
rf={ΔfΔmax   for Δf>0ΔfΔmin   for Δf0
(914)
Venkatakrishnan

Venkatakrishnan [292] 限制器针对面给出:

37. EQUATION_DISPLAY
αf=2rf+1rf(2rf+1)+1
(915)
Min-Mod

Min-Mod 限制器针对面给出:

38. EQUATION_DISPLAY
αf=min(1,1/rf)
(916)
已修正 Venkatakrishnan

已修正 Venkatakrishnan 限制器针对面给出:

39. EQUATION_DISPLAY
αf=rf(2rf+1)+1rf(rf(2rf+1)+1)+1
(917)

Venkatakrishnan 限制器为默认选项。某些情况下,与默认选项相比,Min-Mod 和已修正 Venkatakrishnan 限制器的耗散更小,尤其在气动声学的非稳态模拟中或使用 LES/DES 进行湍流建模时。

α 的网格单元值由以下公式给出:

40. EQUATION_DISPLAY
α=min (αf)
(918)
总变化边界梯度

梯度过度受限时,收敛可能缓慢甚至停止。此问题最常发生于近常数场的梯度中,或者当限制器在 0 和 1 之间循环而没有停在收敛值上时。

要避免出现此问题,可以应用总变化边界梯度限制。重构值的变化源于额外的差值

41. EQUATION_DISPLAY
δ=ψmax(Δmax-Δmin)
(919)

其中, ψ 为介于 0 和 1 之间的分数, max(Δmax-Δmin) 为模拟中任意位置处的局部最大值与局部最小值之间的最大差值。

ψ 的默认值为 0.05。 ψ 的值越大,意味着梯度的受限程度越低。这可以提高精度,尤其是对于平滑流场梯度重构。但是,值过大可能会在流场不连续处产生较大的超调或欠调,从而导致求解器稳定性降低以及精度降低。

TVB 梯度限制适用于使用大涡模拟 (LES) 或分离涡模拟 (DES) 的高保真度模拟,如气动声学和湍流模拟。

边界条件

用于计算上述对流和扩散通量的前述表达式适用于所有内部网格单元面。在与计算域的边界重合的面上,应用边界条件。对边界表面的积分表示为已知边界数据和内部网格单元值的函数。

对于狄利克雷边界条件:
  • 计算对流通量的方法为,将 Eqn. (882) 中的 替换为边界值ϕf ϕb
  • 计算扩散通量的方法为,在 Eqn. (898) 中,将 替换为边界值 并将 替换为 ϕ1ϕbϕ¯ϕ0

对于诺伊曼边界条件,扩散通量直接使用边界处的变量值计算,此类变量值通过离散化梯度近似值获得。

代数方程组

Δϕp=ϕpk+1-ϕpk 时,已传输变量 ϕ 的代数系统以增量形式编写:

42. EQUATION_DISPLAY
apωΔϕp+n anΔϕn=r
(920)

其中, ϕ 是迭代 k+1 (当前迭代)或 k (上一迭代)中的已传输变量。将计算网格单元 p 的所有相邻网格单元 n 的总和。右侧 r 为残差。系数 和 直接从离散项求得。 a p a n ω 为亚松弛因子。

右侧:

43. EQUATION_DISPLAY
r=[ddt(ρϕV)+f(ρϕva)ff(Γϕa)fSϕV]=0
(921)

称为残差,表示原始方程 (Eqn. (881)) 在迭代 k 处的离散形式。之后,基于定义,离散方程完全满足时,残差将变为零。

由于问题的非线性性质,因此需要进行迭代求解。有两种级别的迭代;控制求解更新的外部迭代环,以及控制线性系统的迭代求解的内环。由于外部迭代重复多次,因此仅在每个外部迭代中就足以近似求解线性系统。线性系统的迭代求解使用代数多重网格 (AMG) 完成。