计算不可压缩流体- NS方程求解算法
在上一篇介紹差分格式可以看到求解不可壓縮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
總結
以上是生活随笔為你收集整理的计算不可压缩流体- NS方程求解算法的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: DetachedCriteria和Cri
- 下一篇: 高版本号chrome安装flashpla