格子Boltzmann

格子Boltzmann 方法模拟C/C 复合材料

颗粒沉积过程

罗思璇

()

Particle Deposition Process Simulation in C/C Composites by Lattice-Boltzmann Method

Luo Sixuan

()

Abstract: Lattice Boltzmann method is used here to study the particle deposition process on C/C composites surface. This method considered the boudary condition change during particle deposition. Finally, the deposition pattern is obtained. Keywords: LB Method; flow-particle coupling; C/C composites; deposition

摘要：本文使用格子Boltzmann 方法研究了固体火箭发动机中C/C 复合材料表面上颗粒的沉积模态。该方法考虑了沉积过程中边界形貌的变化对流场的影响，最终得到了颗粒在碳纤维表面的沉积形态。 关键词：LB 方法；流固耦合；C/C 复合材料；沉积

0 引言

C/C 复合材料是目前新材料领域重点研究和开发的一种新型超高温热结构材料，具有密度小，比强度大、热膨胀系数低、热导率高等特点，是理想的航空航天高温材料[1, 2]。

C/C 复合材料在工作过程中其表面流过的工质为高温燃气。高温燃气中通常带有燃烧产生的固体颗粒，如选用较高比冲的含铝推进剂时会产生一定量的凝聚相(Al2O3颗粒)。固体颗粒在C/C 复合材料表面的沉积、冲刷及烧蚀会造成材料内型面的破坏，甚至影响气动性能。

本文使用格子Boltzmann 方法模拟C/C 复合材料中碳纤维上颗粒沉积过程及形态。

1模拟流场的格子Boltzmann 模型

格子Boltzmann 方法是近二十年来刚发展起来的，一种以“半晶格分离法”为处理方式的新型热量逐级传递数值方法,最初是在研究电磁场中的流动现象时被提出的，并且该方法可以确定流体域、固体域和温度场在边界处的连续性,十分适合针对复杂几何形状流固耦合传热问题的数值分析。与传统的经典CFD 方法相比，格子波尔兹曼算法具有很多优点。因而近年来受到国内外学者的广泛关注,并迅速在气固两相流和传热等研究领域得到应用。

格子Boltzmann 方法将流体抽象为微观的虚拟颗粒，通过这些颗粒在规则的网格点上进行碰撞和迁移来达到模拟流场的目的。分布函数f i (x ,t )表示t 时刻，x 网格点上，速度为c i 流体颗粒的概率密度，流场的宏观量通过对分布函数进行统计而得到。本文使用D3Q15模型模拟流场，流体宏观密度ρ和动量ρu 计算如下：

10

Q i i f ρ-==∑，1

Q i i i f ρ-==∑u c

(1)

本文使用BGK 碰撞算子[3]，流场演化方程为：

eq (,)(,)[(,)(,)]i i i i i f x t t t f x t f x t f x t τ+⋅∆+∆-=-c

(2)

其中∆t 为时间步长，τ为无量纲松弛时间，eq i f 为平衡态分布函数，在D2Q9模型中如下计算：

2

eq

2222s s s

1[1()]22i i i

i f c c c ρα⋅⋅=++-c u c u u

(3)

其中u 为流体粒子速度，ρ为流体密度；αi 是与模型有关的权系数，D2Q9模型中，α0=4/9，αi =1/9（i=1,3,5,7），αi =1/36（i=2,4,6,8）；c s 是当地声速，3s c c ，c =∆x/∆t 。流体粘度和压力计算公式为：

2

2

(21)

2s s c P c ντρ⎧=-⎪⎨

⎪=⎩

(4)

2描述颗粒运动的CA 概率模型

Dupuis 和Chopard 最早提出了LB-LGA 模型[4, 5]，并将其运用到雪的沉积以及水底管道下方沙层的冲

刷等定性描述；Przekop 等[6, 7]则利用这种模型定性研究单纤维捕集颗粒过程、纳米和微米尺度混合纤维捕集颗粒过程。但是这些工作均认为颗粒速度取决于流体微团速度，两者之比为某经验常数，无法定量考虑流体对颗粒的曳力以及其他外力对颗粒运动的影响，无法正确描述颗粒与流体微团之间的轨迹滑移、颗粒与流体之间的相互作用，因此只能定性地描述气固两相流动。本文使用随机布朗力[8, 9]的方法来描述颗粒的随机布朗扩散运动，因此颗粒的受力方程如下(本文只考虑流体曳力与随机布朗力)：

p p

D B p

d dt

τ-=+=

+u u u F F (5) p p dx u dt

=

(6)

其中u p 为颗粒速度，τp 为颗粒弛豫时间尺度，τp =ρp d p 2/(18μ) ，μ为气体动力粘度，F B 为随机布朗力，ς为

平均值0、方差为1的高斯随机数，d p 为颗粒直径，k B 为Boltzmann 常数，T 为温度。

通过对方程(5)进行二次积分计算,可以依次得到颗粒的速度和位移如下

1p p f B p p

p

exp()()(1exp())n n t

t

u u u F τττ+∆∆=⋅-

++⋅⋅--

(7) 1p p p f f p B p p p

()(1exp())((1exp()))n n n t t

x x u u u t t F ττττ+∆∆=+---+⋅∆+∆+--⋅⋅⋅

(8)

上标n 表示当前时刻，n+1表示下一时刻。

由此可得到颗粒在Δt 内的实际位移Δx (=1p p n n

x x +-)，把Δx 在某网格方向的投影与该方向网格长度的比值作为颗粒沿该网格线运动到相邻格点的概率：

)7,5,3,1(),,

0max(=⋅∆=i dx

x p i

i e (9)

其中dx 为网格步长。最终可以确定颗粒的最终格点位置:

1p p 11335577n n x x μμμμ+=++++e e e e

(10)

其中μi 为一个布尔量，取1的概率为p i 。

以图 1为例，p 1>0, p 3>0, p 5=0, p 7=0(因为e i =-e i +4，p i ×p i +4=0)。此时通过两个均与分布于[0,1]区间的随机数r 1和r 2，可确定颗粒最终位置，

*

1123p p

*

1123p p 1*

1123p p 3*1123p p 2,,,,,,,,r p r p r p r p t

r p r p t r p r p t

⎧>>=⎪

⎪<>=+∆⎪⎨><=+∆⎪⎪<<=+∆⎪⎩x x x x e x x e x x e 如果如果如果如果 (11)