CFD有限体积7-瞬态项离散

最新推荐文章于 2024-10-21 10:48:54 发布

觅°

最新推荐文章于 2024-10-21 10:48:54 发布

阅读量913

点赞数

分类专栏： CFD 笔记文章标签：学习

本文链接：https://blog.csdn.net/qq_41817886/article/details/129383127

版权

笔记同时被 2 个专栏收录

11 篇文章

订阅专栏

CFD

8 篇文章

订阅专栏

参考资料：

F. Moukalled, L. Mangani, M. Darwish. The finite volume method in computational fluid dynamics. An advanced introduction with Openfoam and Matlab[M]. Chapter13：Temporal discretization: the transient term

简介：本章介绍了一种瞬态项有限积分离散方法（书上介绍了两种，本文省略了有限差分法），基本思想是在一个拟时间单元上采用与对流项离散类似的方法对瞬态项离散。

1. 瞬态项离散

瞬态问题的控制方程为

$\frac{\partial \rho \phi}{\partial t}+\iota (\phi)=0$

经空间积分和离散后，上式可写为

$\frac{\partial \rho_C\phi_C}{\partial t}V_C+L(\phi_C^{t})=0$

其中 $L(\phi_C^{t})$ 是在一个参考时刻t时的非瞬态项空间离散结果，可表达为

$L(\phi_C^{t})=a_C\phi_C^{t}+\sum_{F\sim NB(C)}a_C\phi_F^{t}-b_C$

当t趋向于无穷大时，瞬态离散方程趋近于定常方程，这一情况与利用时间推进法求解定常问题一致。现介绍瞬态项离散方法。

瞬态项的有限体积离散方法与对流项的离散方法十分相近，除了积分是对时间积分而不是对空间积分。对 $\frac{\partial \rho_C\phi_C}{\partial t}V_C+L(\phi_C^{t})=0$ 从 $[t-\Delta t/2, t+\Delta t/2]$ 进行积分，可得

$\int_{t-\Delta t/2}^{t+\Delta t/2}\frac{\partial \rho_C\phi_C}{\partial t}V_Cdt+\int_{t-\Delta t/2}^{t+\Delta t/2}L(\phi_C)dt=0$

将VC视为常数，上式结果为

$V_C(\rho_C\phi_C)^{t+\Delta t/2}-V_C(\rho_C\phi_C)^{t-\Delta t/2}+L(\phi_C)\Delta t=0$

等价于

$\frac{V_C(\rho_C\phi_C)^{t+\Delta t/2}-V_C(\rho_C\phi_C)^{t-\Delta t/2}}{\Delta t}+L(\phi_C)=0$

现需要根据 $t-\Delta t$ 和t处（单元上）参数值确定 $t-\Delta t/2$ 和 $t+\Delta t/2$ 处（界面上）参数值（也就是确定插值型线），最后离散方程为

$FluxT=FluxC\phi_C+FluxC^{k-1}\phi_C^{k-1}+FluxV$

上标k-1表示上一时间步计算结果，最后更新代数方程组系数：

$\\a_C=a_C+FluxC\\b_C=b_C-FluxC^{k-1}\phi_C^{k-1}-FluxV$

1）一阶瞬态格式

一阶隐式和显式欧拉格式分别基于一阶迎风格式和背风格式。

a. 一阶隐式欧拉格式

首先介绍基于一阶迎风格式的隐式欧拉格式。在隐式欧拉格式下，界面上参数值取上一时间步的参数值：

$\\(\rho_C\phi_C)^{t+\Delta t/2}=(\rho_C\phi_C)^{t}\\(\rho_C\phi_C)^{t-\Delta t/2}=(\rho_C\phi_C)^{t-\Delta t}$

因此 $\frac{V_C(\rho_C\phi_C)^{t+\Delta t/2}-V_C(\rho_C\phi_C)^{t-\Delta t/2}}{\Delta t}+L(\phi_C)=0$ 取为

$\frac{V_C(\rho_C\phi_C)^{t}-V_C(\rho_C\phi_C)^{t-\Delta t}}{\Delta t}+L(\phi_C)=0$

由于该格式为一阶格式，会产生数值耗散（上标圆圈与k-1意义一致），由泰勒展开式得

整理为

代入离散方程 $\frac{V_C(\rho_C\phi_C)^{t}-V_C(\rho_C\phi_C)^{t-\Delta t}}{\Delta t}+L(\phi_C)=0$ 中，有

该格式即使在时间步长很大的情况下，依旧是稳定的。

b. 一阶显式欧拉格式

一阶显式欧拉格式采用的是背风插值型线，即界面 $t-\Delta t/2$ ， $t+\Delta t/2$ 上参数值取下一时间步参数值：

$\\(\rho_C \phi_C)^{t+\Delta t/2}=(\rho_C\phi_C)^{t+\Delta t}\\(\rho_C\phi_C)^{t-\Delta t/2}=(\rho_C\phi_C)^{t}$

因此 $\frac{V_C(\rho_C\phi_C)^{t+\Delta t/2}-V_C(\rho_C\phi_C)^{t-\Delta t/2}}{\Delta t}+L(\phi_C)=0$ 取为

$\frac{V_C(\rho_C\phi_C)^{t+\Delta t}-V_C(\rho_C\phi_C)^{t}}{\Delta t}+L(\phi_C^{t})=0$

该格式为显式格式，其包含假设： $\rho\phi$ 的对流在整个时间步始终占主导。该格式会产生数值反耗散，t+dt对t做泰勒展开，为

整理为

代入离散方程 $\frac{V_C(\rho_C\phi_C)^{t+\Delta t/2}-V_C(\rho_C\phi_C)^{t-\Delta t/2}}{\Delta t}+L(\phi_C)=0$ 中，有

与数值反耗散相关的是数值不稳定性，且随着时间步dt的增大，数值不稳定性增强。

2）二阶瞬态欧拉格式

二阶瞬态格式同样可以采用线性插值型线构造，i.e., 可以利用对称线性型线（中心差分）构造出Crank-Nicolson (CN)格式或利用考虑远上游单元的线性插值型线构造出Adams-Moulton格式，该格式也被称为二阶迎风欧拉格式（SOUE）。

a. CN格式（中心差分型线）

当时间步均匀时，界面上参数值可以表达为

$(\rho \phi)^{t+\Delta t/2}=\frac{(\rho \phi)^{t+\Delta t}+(\rho \phi)^{t}}{2}$

$(\rho \phi)^{t-\Delta t/2}=\frac{(\rho \phi)^t+(\rho \phi)^{t-\Delta t}}{2}$

代入离散方程 $\frac{V_C(\rho_C\phi_C)^{t+\Delta t/2}-V_C(\rho_C\phi_C)^{t-\Delta t/2}}{\Delta t}+L(\phi_C)=0$ 中，有

该格式由t-dt和t时间步参数值求解t+dt时间步参数值，为显式方法，需要CFL数限制时间步长以保证格式的稳定性。

现分析CN格式的数值精度，将t+dt和t-dt参数对t参数进行泰勒展开

代入离散方程中，有

$\frac{\partial (\rho \phi)}{\partial t} |_t+\frac{\partial ^3(\rho \phi)}{\partial t^3}|_t\frac{\Delta t^2}{6}+o(\Delta t^4)+\frac{L(\phi_C^t)}{V_C}=0$

$\frac{\partial (\rho \phi)}{\partial t} |_t+\frac{L(\phi_C^t)}{V_C}=-\frac{\partial ^3(\rho \phi)}{\partial t^3}|_t\frac{\Delta t^2}{6}-o(\Delta t^4)$

可见该格式为二阶精度（书上截断误差给的是dt^3，有问题？）。

b. 二阶迎风欧拉格式（SOUE）

二阶迎风格式：

$\phi_f=\frac{3}{2}\phi_C-\frac{1}{2}\phi_U$

其中C为界面上游，U为界面远上游单元。因此

$(\rho_C \phi_C)^{t+\Delta t/2}=\frac{3}{2}(\rho_C\phi_C)^t-\frac{1}{2}(\rho_C\phi_C)^{t-\Delta t}$

$(\rho_C \phi_C)^{t-\Delta t/2}=\frac{3}{2}(\rho_C\phi_C)^{t-\Delta t}-\frac{1}{2}(\rho_C\phi_C)^{t-2\Delta t}$

代入离散方程 $\frac{V_C(\rho_C\phi_C)^{t+\Delta t/2}-V_C(\rho_C\phi_C)^{t-\Delta t/2}}{\Delta t}+L(\phi_C)=0$ 中，有

已知t-2dt和t-dt参数，通过上函数式求解t时刻参数，为隐式格式

现分析SOUE格式的数值精度。将t-dt和t-2dt参数对t参数做泰勒展开：

代入离散方程中

因此SOUE格式为二阶精度格式。

2. 初始条件

采用以下时间离散“网格”

一阶迎风格式决定了

$\\(\rho_C \phi_C)^{t_{initial}+\Delta t/2}=(\rho_C \phi_C)^{t_{initial}}\\(\rho_C \phi_C)^{t_{initial}+\frac{3}{2}\Delta t}=(\rho_C \phi_C)^{t_{initial}+\Delta t}$

3. 非均匀时间步长

在实际应用中常采用变时间步长，在大的时间步长不能违反CFL条件。CFL条件：Courant等人最先发展了CFL条件，根据符号相反原则，要求

$a_C+a_C^{o}\leq 0$

例如对于一维对流问题，离散方程系数为

$\\a_C=\rho_C^{o}u_C^{o}\Delta y_C\\a_C^{o}=-\frac{\rho_C^{o}V_C}{\Delta t}\\a_C+a_C^{o}\leq 0$

等价于

$\Delta t\leq \frac{\Delta x_C}{u_C^{o}}$

对流问题CFL为

$CFL^{conv}=\frac{|v_C^{o}|\Delta t}{\Delta x_C}$

因此一维对流问题CFL应小于等于1。

有限体积法的非均匀时间步：时间单元定义为 $\Delta t$ ，两个相邻瞬态单元的时间步长为 $\delta t$ ，对于均匀时间步长有 $\Delta t=\delta t$ ，对于非均匀时间步长， $\delta t=(\Delta t+\Delta t^{o})/2$ 。

以CN格式和SOUE格式为例说明非均匀时间步长的应用方法。

1）CN格式

均匀网格时，CN格式为

$(\rho \phi)^{t+\Delta t/2}=\frac{(\rho \phi)^{t+\Delta t}+(\rho \phi)^{t}}{2}$

$(\rho \phi)^{t-\Delta t/2}=\frac{(\rho \phi)^t+(\rho \phi)^{t-\Delta t}}{2}$

非均匀网格时，上一时间步和当前时间步的参数值都会影响格式插值型线，如

$(\rho_C\phi_C)^{t-\Delta t/2}=\frac{\Delta t^{o}}{\Delta t+\Delta t^{o}}(\rho_C\phi_C)^t+\frac{\Delta t}{\Delta t+\Delta t^{o}}(\rho_C\phi_C)^{t-(\Delta t^{o}+\Delta t)/2}$

$(\rho_C \phi_C)^{t-\Delta t/2-\Delta t^{o}}=\frac{\Delta t^{oo}}{\Delta t^{o}+\Delta t^{oo}}(\rho_C \phi_C)^{t-(\Delta t^{o}+\Delta t)/2}+\frac{\Delta t^{o}}{\Delta t^{o}+\Delta t^{oo}}(\rho_C\phi_C)^{t-\Delta t^{o}-(\Delta t+\Delta t^{oo})/2}$

代入离散方程

$\frac{(\rho_C\phi_C)^{t-\Delta t/2}-(\rho_C\phi_C)^{t-\Delta t/2-\Delta t^{o}}}{\Delta t}V_C+L(\phi_C^{o})=0$

得到

$\frac{\Delta t^{o}}{\Delta t+\Delta t^{o}}\frac{V_C}{\Delta t}(\rho_C \phi_C)+[\frac{\Delta t}{\Delta t+\Delta t^{o}}-\frac{\Delta t^{oo}}{\Delta t^{o}+\Delta t^{oo}}]\frac{V_C}{\Delta t}(\rho_C \phi_C)^{o}-\frac{\Delta t^{o}}{\Delta t^{o}+\Delta t^{oo}}\frac{V_C}{\Delta t}(\rho_C \phi_C)^{oo}+L(\phi_C^{o})=0$

2）SOUE格式