Hybrid fluid-solid simulation
名工大:尾形修司
流体-弾性体連成問題: 古典的問題設定
設定
◉ 流体領域 :
!
"f
r ˙ v
f=# $
t %
f
◉ 固体領域 :
!
" S
r ˙ v S =# $
t %
s
◉ 境界 :
!
r v f
=r v s 速度一致(粘着条件)
!
t "
f# ˆ n
f+
t "
s# ˆ n
s= 0 応力のつりあい
!
"
!
"s
!
"f
!
"
→ Navier-Stokes方程式
→ 粗視化した弾性体
現実は様々な可能性あり(接触角,疎水性/親水性)
現実は塑性変形もありうる
!
"s!
"f
流体-弾性体連成問題: 数値解法例
直接数値シミュレーション(DNS)法
◉ 流体領域でのみ,Navier-Stokes方程式をメッシュ上で解く
◉ 流体-固体境界の移動・変形の取り扱いは,一般に困難である.
例えばArbitrary Lagrangian-Eulerian (ALE) methodにより弾性体の移動・変形に伴い,メッシュを動的に再構成する:あるときは留まって(Euler的),またあるときは動いて(Lagrange的),観測する
流体-弾性体連成問題: 数値解法例
Immersed boundary method (IBM)
◉ 流体は,全領域にNavier-Stokes方程式を適用し,Euler的に解く
C.S. Peskin, J. Comp. Phys. 25, 220 (1977)
◉ 境界上の粘着条件は,Navier-Stokes方程式の 外力項を通じてフィードバックループにより 反復的に実現する
◉ 境界上の固体点で流体に働く力を 流体格子点上で流体に対して働く力は
!
r F (
r X ,t)
!
r f (
r x ,t) =
r F (
r X (s,t),t)"(
r x #
r X (s,t))$ ds
として
実際の数値計算ではデルタ関数を離散化する
!
"(x,x ') =1# | x # x' |, for | x # x ' |< 1
0, otherwise
$ % &
!
"#r u
#t+
r u $ %
r u
&
' (
)
* + +%" = µ,
r u +
r f
…境界上の固体点の移動
!
"r X
"t=
r u =
r u (
r X (s,t), t# )$(
r x %
r X (s,t))ds
…外力を含むNS方程式◉
!
r X
!
r x
Boltzmann方程式からの流体方程式導出 (1)
数密度n,温度T, 平均速度uの条件での,速度vのMaxwell-Boltzmann分布
!
f"0
dr v = n …(1)
!
f0(v"# $ u" )d
r v = 0
!
f0(v"# $ u" )(v% $ u% )d
r v = nkBT&" ,%
!
f0(v"# $ u" )(v% $ u% )(v& $ u& )d
r v = 0
!
f0(v"# $ u" )(v% $ u% )(
r v $
r u )
2dr v = 5n(kBT)
2&" ,%
GaussIntegrals:
…(2)
…(3)…(4)
…(5)
!
f0
=n
(2"kBT)3 / 2exp #
(r v #
r u )
2
2kBT
$
% &
'
( )
Boltzmann方程式 (元は希薄ガスを想定)
!
"t +r v # "r
x +
r F # " r
v ( ) f = d $ v 1% d $ v
2dv
2( $ f
1$ f 2& f
1f
2)P
12' $ 1 $ 2
Bhatnagar-Gross-Krook近似 (単一緩和時間近似)
!
"t f +r v # " r
x f +
r F # "r
v f =
1
$( f
0% f )
!
" f # f0$ %(&t f
0+
r v ' & r
x f
0+
r F ' &r
v f
0) …(6)
↑外力:
!
r F
は, と のずれの程度
!
"
!
f
!
f0
↑2体散乱
を表していると見なせる
程度での長波長極限の挙動をこれから考える
!
"1
Boltzmann方程式からの流体方程式導出 (2)
質量保存
!
"t dr v # f + "$ d
r v # fv$ + F$ d
r v "v$# f =
1
%dr v # ( f
0 & f )
!
" #tn + #$ (nu$ ) = 0
運動量保存
BGK-Boltzmann方程式の両辺を速度積分すると
BGK-Boltzmann方程式の両辺に をかけて積分すると
!
v"
!
"t dr v fv#$ + "% d
r v fv#$ v% + F% d
r v "v%$ fv# =
1
&dr v ( f
0 ' f )v#$
!
" #t (nu$ ) + #% dr v fv$v%& ' nF$ = 0
…(7)
…(8)
式(8)の第2項は,式(6)を適用して部分積分すると
!
dr v fv"# v$ % d
r v f
0# v"v$ & ' ((t dr v f
0# v"v$ + () dr v # f
0v"v$v) + nF"u$ + nu"F$ )
!
dr v " f
0v#v$ = d
r v " f
0u#u$ + nkBT%# ,$式(2),(3)より なので
!
" 0の項を集めると
!
"t(nu# ) +"$ (nu#u$ ) = %"# (nkBT) + nF#
↑Navier-Stokes方程式において粘性項をゼロとした,Euler方程式が得られた.
↑連続の式が得られる
…(10)
さらに式(7)を用いて,
!
"tu# + u$"$u# = %
1
n"# (nkBT) + F# …(11)
…(9)
式(1)と(2)を使うと
式(8)の
↓部分積分
(同一文字はsummation)
(微分は直後にのみかかる)
Boltzmann方程式からの流体方程式導出 (3)
!
"t ( dr v f
0v#$ v% ) = "t (nu#u% + nkBT&# ,% )
= "t (nu# )u% + nu#"tu% + "tnkBT&#,% + n"t (kBT)&# ,%
= '"( (nu#u( )u% '"# (nkBT)u% ' nF#u%
' nu#u("(u% ' u#"% (nkBT) ' nu#F%
'"( (nu( )kBT&#,% ' nu("( (kBT&#,% ) ' 2
3"(u(nkBT
= '"( (nu#u%u)) '"% (nkBT)u# '"# (nkBT)u%
' n(F#u% + u#F% ) '"( (nkBTu))&# ,% '2
3nkBT"(u(
!
"# dv$ f0v%v&v# = "& (nkBTu% ) + "% (nkBTu& ) + "# (nkBTu# )'%,& + "# (nu%u&u# )
↓式変形には,式(4)を使った
!
"t dr v f
0v#$ v% + "& d
r v f
0$ v#v%v& = nkBT("#u% + "%u# ) '2
3nkBT"&u&よって
!
n"tu# + nu$"$u# = %"# (nkBT) + nF# + "$ [&("$u# + "#u$ %
2
3"'u'(#,$ )]
!
with " # nkBT$ % consider as viscocity ↑有名なNavier-Stokes方程式(Newtonian)!
項中:
!
"1式(8)の
式変形には,式(10)と(11)を使い,さらにエネルギー保存から導ける
!
"t(k
BT) + u#"# (kBT) $ %
2
3"#u#kBT
を使った
←
式(8)の
!
"1項をまとめると
!
" 0項と
また
熱輸送方程式もBoltzmann方程式に を掛け積分して得られる
!
(r v "
r u )
2
粘性率,熱伝導度,音速は,現実とは関係無く互いに で結びついている!!
!
"
Lattice Boltzmann method (LBM) (1)
◉ コーディングが楽な流体計算法(非圧縮流体に適する)
◉ 常に格子点上のみに存在し,離散化速度(Δtで1 or 0格子点移動)をもつ, 多数の仮想粒子
◉ 空間に規則的な格子点を設定(LB方程式の粗視化に相当する)
◉ 仮想粒子は,移動 & 散乱で質量保存則,運動量保存則に従う
!
r e
i
!
fi(r x +
r c i"t, t + "t) # f i
*(r x ,t) = f i(
r x ,t) +$i({ f i(t)}) + Fi
!
fi(r x ,t)時刻t で離散化速度 を持つ粒子数
!
r c
i
衝突項
!
"i
外力項
!
Fi
左例(D2Q9)では速度0も入れて i={0,…,8}
参考: 蔦原等「格子気体法・格子ボルツマン法」 コロナ社(1999).
!
r x
◉ bounce-back rule等により,複雑境界へも,簡単に適用できる方法
!
"x(lattice spacing)
Lattice Boltzmann method (LBM) (2)
◉ 巨視的スケールでは,BGK-Bolzmann方程式と同様に, LBMも連続の式,Navier-Stokes方程式を満たす
等温流体でのBGK衝突モデル
!
"i = #( fi # f ieq) /$
!
fi
eq= f j
j
"#
$ % %
&
' ( ( wi 1+
r c i )
r v
cs
2+(r c i )
r v )
2
2cs
4*
v2
2cs
2
+
, -
.
/ 0
温度:
!
T
仮想粒子質量:
!
m
!
"r v = m fi
r c i
i
#局所流速 :
!
r v
!
" = m fii
#局所密度: !
wi
重み:
!
wi=1
i
"#
$ %
&
' (
Maxwell-Boltzmann分布を,流速vの2次まで展開して得る平衡分布
緩和時間:
!
"#t
!
"x
"t
2# $ 0.5
3◉ 動粘性率:
!
cs
2
="x
2
3"t2
音速cs:
単一時間緩和係数
…単一時間緩和係数に関係する
!
fieq
= ni
"
!
fieq(vi," # u" ) = 0
i
$
!
fieq(vi," # u" )(vi,$ # u$ ) = nkBT%" ,$
i
&
等を満たすように決める
LBMでの複雑境界の取り扱い法
Bounce-back rule
擬似摩擦力
◉ 仮想粒子の移動後に,固体内部での流れを打ち消す速度分布を追加する
◉ 壁の移動速度 (UB)を反映させる
e.g., Buxton et al, Phys. Rev. E 71, 56707 (2005)
Ahlrichs and Dunwerg,J. Chem. Phys. 111, 8225 (1999)
◉ 流体点と固体点の速度差に比例して,摩擦力を互いに働かせる
!
r F s = "# (
r v s "
r v f )
!
r F f
= "r F s
…固体への摩擦力
…流体への摩擦力
◉ 固体の内部に流体は侵入しない
◉ 注意:固体の内部に流体が侵入する
!
fi(r x i,t) = f i(
r x i,t) + "f i(
r x i,t,
r U B)
v UB
−(v −UB)+UB= −v+2UB
UB
!
m << MB
(反発係数1)
Preliminary problem
固定壁に挟まれたポアズイユ流内に,rodを置いておく.
Rod=粗視化粒子(CGP)法で作られた粒子群
LBM–CGP interaction in 2D (1)
◉ Immersed boundary methodの考え方を援用する
… 3rd-order Lagrange interpolation
◉ LBM流体速度分布の固体壁での反射 (局所的な固体面方向の情報を用いない簡単な手法)
… bounce back with CG velocity (assuming “heavy” solid)
◉ LBM流体からCG固体への運動量移行 (at X)
LBM–CGP interaction in 2D (2)
!
r F (
r x ,t) =
r g (
r X (s,t), t)"(
r x #
r X (s,t))ds
$S%
!
fi(r X ,t) = "(
r X ,
r x #,$
#,$
% ) f i(r x # ,$ ,t)
… 固体から 流体への力
!
f j (r X ,t + "t) = f i
*(r X ,t) # 2 fk
*
k$( )wi
r c i %
r U CG(r X )
cs
2
!
2mfi
*(X,t) " 2m( fk
*
k# )wi
r c i $
r U CG(
r X )
cs
2
%
& '
(
) *
i
collision
#r c i
X(s,t)g(X,t)
x
!
redistribute : f i(r X ,t + "t)# f i(
r x ,t + "t) … to LBM points
!
"S
[方向jは,方向iの逆向き]
Conservation of total energy
10 CGP
20 CGP
= Ar (atom)× 1024 21 CG particle
Dynamics of CG particles alone
collision
streaming
CG-dynamics
apply force ufluid +F Δt/ρ
boundary condition on the wall
fluid force on CG particles
get macro fields
CG-evolution
velocity and position of CG particles
Lagrange interp.
scaling to CG world
rescaling to LBM world
Time advancing algorithm
!
r U
CG[LBM unit] =
vatomic unit
unit
3cs
t M
"1r P
CG[atomic unit]
!
FCG,ext
=3"c
s
2#x
fatomic unit
unitf [LBM unit]
Scaling of velocity:
Scaling of force:
400
= Ar (atom)× 102421 CG particle
1 dX[CGP] = (0.8 or 1.2)dx[LBM]
dX
dx
20 (10) CGP
50 CGP
H=100Umax
Characteristics of fluid-solid system
!
Re =U
max
µH = 400 or 500
!
Umax = (0.1 or 0.15) 3cs
Re = 400
Hybrid LBM-CGP simulation
Hard rod
Internalstressof rod
Re = 500
Hybrid LBM-CGP simulation
Soft rod
Internalstressof rod
Re=500Re=400
Energy in the solid body
Discussion: Scaling of time
MD計算: Ar atom (Lennard-Jones potential)
CGP計算
LBM計算
TCG
τ
理論値 : 400 ps v.s. 1ps (Air) Too far!採用値: Divide 1 LBM step into 5 CGP steps
Time scale!
c" = dx(LBM) = 3cs"
!
" =dx
3cs
~ 400ps for air, 100ps for liquid
~ 400#tCG for air, 100#tCG for liquid
今回は,現実の気体より1-2桁程度早い音速に設定した.[つまり,現実の固体より特性振動時間が1-2桁程度長い,柔らかい固体を設定したことに相当するとも言える]
d=空間次元