计算不可压缩流体- NS方程求解算法
2014-03-20 02:11
253 查看
在上一篇介绍差分格式可以看到求解不可压缩NS方程,主要是求解速度的动量方程和压力的波松方程。本文介绍一些历史中的算法,下一文介绍目前常用的算法。
一 Marker and Celler (MAC) 算法
在自然差分格式中给出PPE如下: laplace ( p) = laplace ( body_force) - rho * div ( u * grad_u) (1) 。 该方程不保证 div_ u = 0 条件。
MAC 将PPE方程如下定义: laplace(p) = miu * laplace( Div) - (du^2/dxx + duv/dxy + dvu/dyx + dv^2/dyy) - dD/dt (2), 即Div 项不直接扔掉,特别是 dD/dt项。
离散(2)式并令 (n+1) 时步的 Div = 0 (保证div _u = 0 条件),得到 (n+1)时步的压力场;将该压力值带入上一篇介绍的半隐式离散动量方程, 采用forward Euler时间离散求解(n+1)时步的速度场。
边界处理同自然差分格式。
二 SOLA算法
MAC每步都要计算PPE,代价大。考虑 压力-速度耦合,且 速度满足不可压缩条件(div_ u = 0), 那么存在隐函数 Div(P) = 0. 构造不动点迭代(牛顿迭代)如下:
p^(m+1) = p^(m) - D/ D'(p) % p^(m) 表示第m迭代步的压力值。
其中 D‘(p) = dD/dp 构造怎么构造呢? 还是从 半隐式离散动量方程入手: 对方程取散度,并令 div_ u ^ (n) = 0 %即n时步的速度满足不可压缩条件, 在对p求导。
du^(n+1) / dp = 2k/ h/ h; dv^(n+1) / dp = 2k / h/ h; (k = dt, h = grid spacing)
故 D' = 4k/ h/ h。 亦即压力迭代格式 p^(m+1) = p ^(m) - 4k D_ij /h /h (3)
由于动量方程显式离散,所以迭代过程压力变化多少,对应速度变化多少。即 n+1 时步
u_(i,j) ^(m+1) = u_(i,j) ^(m) + k/h * dp
v_(i,j) ^(m+1) =v_(i,j) ^(m) + k/h * dp
为满足单元边界通量守恒,同时更新n+1时步的 u_(i-1,j ) , v(i, j-1):
u_(i -1 ,j) ^(m+1) = u_(i-1,j) ^(m) - k/h * dp
v_(i,j -1 ) ^(m+1) = v_(i,j-1) ^(m) - k/h * dp
注意这里变化是负号了。
带入上面四式到FV离散的散度项,即 D^(m) = D^(m+1) - 4k * dp/ h/ h
令Div^(m+1) = 0 --> dp = - 4k D_Ij/ h / h (4)
(3) = (4)
三 人工可压缩性
在div_u = 0 方程中加上 drho/dt, 求解至稳态,得到 rho, 定义 drho/ dp = delta % delta压缩系数 --> 压力值
应用较广,诸如SPH等粒子算法中。但是,这个稳态密度对应的压力,其实还是没有物理意义的,因为它跟动量方程里面的压力没半点关系。或者说 delta不是一个动量方程相关的量。所以,经验值了。
下一篇介绍后两者: 投影法 SIMPLE
一 Marker and Celler (MAC) 算法
在自然差分格式中给出PPE如下: laplace ( p) = laplace ( body_force) - rho * div ( u * grad_u) (1) 。 该方程不保证 div_ u = 0 条件。
MAC 将PPE方程如下定义: laplace(p) = miu * laplace( Div) - (du^2/dxx + duv/dxy + dvu/dyx + dv^2/dyy) - dD/dt (2), 即Div 项不直接扔掉,特别是 dD/dt项。
离散(2)式并令 (n+1) 时步的 Div = 0 (保证div _u = 0 条件),得到 (n+1)时步的压力场;将该压力值带入上一篇介绍的半隐式离散动量方程, 采用forward Euler时间离散求解(n+1)时步的速度场。
边界处理同自然差分格式。
二 SOLA算法
MAC每步都要计算PPE,代价大。考虑 压力-速度耦合,且 速度满足不可压缩条件(div_ u = 0), 那么存在隐函数 Div(P) = 0. 构造不动点迭代(牛顿迭代)如下:
p^(m+1) = p^(m) - D/ D'(p) % p^(m) 表示第m迭代步的压力值。
其中 D‘(p) = dD/dp 构造怎么构造呢? 还是从 半隐式离散动量方程入手: 对方程取散度,并令 div_ u ^ (n) = 0 %即n时步的速度满足不可压缩条件, 在对p求导。
du^(n+1) / dp = 2k/ h/ h; dv^(n+1) / dp = 2k / h/ h; (k = dt, h = grid spacing)
故 D' = 4k/ h/ h。 亦即压力迭代格式 p^(m+1) = p ^(m) - 4k D_ij /h /h (3)
由于动量方程显式离散,所以迭代过程压力变化多少,对应速度变化多少。即 n+1 时步
u_(i,j) ^(m+1) = u_(i,j) ^(m) + k/h * dp
v_(i,j) ^(m+1) =v_(i,j) ^(m) + k/h * dp
为满足单元边界通量守恒,同时更新n+1时步的 u_(i-1,j ) , v(i, j-1):
u_(i -1 ,j) ^(m+1) = u_(i-1,j) ^(m) - k/h * dp
v_(i,j -1 ) ^(m+1) = v_(i,j-1) ^(m) - k/h * dp
注意这里变化是负号了。
带入上面四式到FV离散的散度项,即 D^(m) = D^(m+1) - 4k * dp/ h/ h
令Div^(m+1) = 0 --> dp = - 4k D_Ij/ h / h (4)
(3) = (4)
三 人工可压缩性
在div_u = 0 方程中加上 drho/dt, 求解至稳态,得到 rho, 定义 drho/ dp = delta % delta压缩系数 --> 压力值
应用较广,诸如SPH等粒子算法中。但是,这个稳态密度对应的压力,其实还是没有物理意义的,因为它跟动量方程里面的压力没半点关系。或者说 delta不是一个动量方程相关的量。所以,经验值了。
下一篇介绍后两者: 投影法 SIMPLE
相关文章推荐
- 计算不可压缩流体--差分格式的两个问题
- 计算不可压缩流体 -- 差分/有限体积离散格式
- 计算不可压缩流体 -- 数学基础
- [物理学与PDEs]第2章习题8 一维定常粘性不可压缩流体的求解
- 人工压缩算法--定常原始变量不可压缩N-S方程
- GPU求解粘性不可压流体
- LA3485二分+求解积分方程+辛普森算法计算积分
- 稀疏线性系统求解算法 之 存储结构(MCRF) 强于二维数组、三元组、行压缩、修正行压缩等
- 赛码网基础算法——简单计算(求解加减乘除等复杂算式)
- [物理学与PDEs]第2章习题7 一维不可压理想流体的求解
- mysql计算生日算法
- Hadoop计算能力调度器算法解析
- 字符串计算 算法开源
- YYMMDD转换成4位字符压缩表示算法java实现——应用各位不等进制的思想
- USACO 蔡勒公式 (计算任何一日属一星期中哪一日的算法)
- 我就给一个PHP逆波兰表达式的算法吧---工资计算专用
- 我就给一个PHP逆波兰表达式的算法吧---工资计算专用
- 计算n!的最后一位非零数字的算法
- 计算几何常用算法概览[z]
- GZIP压缩原理分析(15)——第五章 Deflate算法详解(五06) 预备知识(05) 预备知识总结