格子Boltzmann数值计算的性能优化
偏微分方程求解的一种新颖方法——格子Boltzmann模型

7 6 邻 节 点
大 学 数 学
第2 7卷
) 撞 , 一 个 节 点 上 从 相 邻 节 点 运 动 来 的 粒 子 发 生 碰 撞 , 据 质 量 、 量 和 能 量 守 恒 规 则 改 碰 在 根 动
其 中 r 松弛 时间 尺度 , 制达 到平衡 的速 度 ( 是 控 可根据 需 要 进行 设 置 ) 由于稳 定 性 的原 因 , 过 实 际测 , 经 算 r必须 大于 1e /.
事 实 上 不 同 的 网 格 剖 分 有 着 不 同 的平 衡 分 布 函数 , B 建 立 模 型 的 核 心 问 题 就 是 根 据 不 同 的 网 格 L M
[ 键 词 ] 格 子 B l ma n方 法 ; 衡 态 分 布 函 数 ; Q 关 ot n z 平 D2 9模 型 ; a i — tk s 程 ; 流一 扩 散 方 程 N ve So e 方 r 对 [ 图 分 类 号 ] O2 1 8 中 4 .2 [ 献标识码]A 文 [ 章 编 号 ] 17 —4 4 2 1 ) 30 7 —8 文 6 21 5 (0 1 0 —0 50
在 低 Mah 马赫 ) 的假 设下 ( l c)其 中粒子平 衡态 分布 函数 为 c( 数 I , U《
~ P
[ + 一 ] +
且
C =c 4 /  ̄,
。 /, 一4 9
1 U 一 3 一(2 一 4 1 9, ∞ 一 6 7 8 1 3 , — / 5 一∞ 一∞ — / 6
第2 7卷 第 3期
格子波兹曼方法

格子波兹曼方法
格子波兹曼方法(Lattice Boltzmann Method, LBM)是一种广泛应用于计算流
体力学领域的数值方法。
它基于分子动力学模型,通过离散化空间网格和时间步长来模拟复杂的流体流动问题。
格子波兹曼方法通过将流体宏观物理量离散化到网格上的节点,使用分布函数
描述流体粒子的运动。
流体粒子在相邻节点之间以一种特定的方式进行碰撞和传播,模拟流体的宏观行为。
格子波兹曼方法相对于传统的Navier-Stokes方程求解方法具有多个优势。
首先,它因其并行化的能力而广泛应用于高性能计算中。
其次,LBM的离散化框架使得
它在处理具有复杂边界条件和多相流问题时更加灵活。
此外,LBM对于非连续和
非均匀流体介质的模拟效果也相对较好。
格子波兹曼方法在各个领域都有广泛的应用。
在流体力学领域,LBM被用于
模拟自由表面流动、湍流现象和多孔介质中的流动行为。
在微观领域,LBM也被
用于模拟微观流体力学现象,例如微管流动和纳米颗粒悬浮体的输运行为。
除了流体力学领域,格子波兹曼方法还被应用于其他科学领域。
例如,它被用
于模拟热传导、传质过程、相变以及复杂物质的输运现象。
此外,LBM还被用于
模拟生物流体力学、地下水流动、大气动力学和地震波传播等问题。
综上所述,格子波兹曼方法是一个高效且灵活的数值方法,用于模拟复杂的流
体流动问题。
它在计算流体力学领域以及其他科学领域都有广泛的应用前景。
这种方法的进一步发展和应用将有助于我们更好地理解和预测流体行为,并解决相关领域的实际问题。
颗粒群绕流传质特性的格子Boltzmann模拟

颗粒群绕流传质特性的格子Boltzmann模拟王利民;吴承优【摘要】为考察气固两相流中颗粒团聚对气固传质过程的影响,采用耦合传质的格子Boltzmann模型对在空间均匀分布的颗粒结构和非均匀分布的颗粒团聚物结构条件下的绕流过程进行了直接数值模拟,得到了气体绕颗粒流动的速度分布和伴随流动的浓度分布,发现颗粒团聚物对绕流过程的速度分布和浓度分布都具有明显的影响.通过对颗粒绕流的气固传质过程进行定量分析发现,在两种结构条件下计算得到的传质舍伍德数均随着雷诺数的增大而呈指数函数形式增大,并且在均匀结构条件下的传质舍伍德数一般为非均匀结构条件下的3~5倍,颗粒团聚物的存在将严重影响颗粒绕流过程的气固传质效率.通过格子Boltzmann方法建立的气固相间传质模型,可以为研究气固两相流介尺度结构的传递特性提供理论依据.【期刊名称】《化学反应工程与工艺》【年(卷),期】2016(032)004【总页数】8页(P289-296)【关键词】颗粒群绕流;颗粒团聚物;格子玻尔兹曼方法;传质;数值模拟【作者】王利民;吴承优【作者单位】中国科学院过程工程研究所多相复杂系统国家重点实验室,北京100190;中国科学院过程工程研究所多相复杂系统国家重点实验室,北京100190;江西理工大学资源与环境工程学院,江西赣州341000【正文语种】中文【中图分类】TQ021.4气固流动中,团聚物对气固传递特性的影响很大[1]。
从20世纪80年代,李静海等[2]开始研究气固反应器层次的介尺度问题,从颗粒聚团现象入手,认为聚团的形成来源于气体和颗粒的各自运动趋势在竞争中的协调,从而建立了能量最小多尺度(Energy-Minimization Multi-Scale,EMMS)模型,通过EMMS模型计算发现稀相和密相的曳力系数存在几个数量级的差异[3],而双流体模型对微元网格内部的平均化处理却忽略了该非均匀结构的影响。
杨宁等[4]利用EMMS 模型考虑了宏观非均匀结构对计算网格的曳力修正,更准确地预测了快速流化床的颗粒夹带量。
格子玻尔兹曼方法

格子玻尔兹曼方法格子玻尔兹曼方法(Lattice Boltzmann Method,简称LBM)是一种基于微观粒子动力学的计算流体力学方法,它是由Lattice Gas Automata(LGA)经过演化和发展而来的。
LBM是一种离散的方法,它通过在空间网格上模拟分子碰撞和传输过程来描述流体的宏观运动。
与传统的有限差分法、有限体积法相比,LBM具有计算效率高、并行性好、适应复杂边界条件等优点,因此在流体力学领域得到了广泛的应用。
LBM的基本思想是将流体系统离散化,将连续的流体宏观运动转化为离散的微观碰撞和传输过程。
在LBM中,流体被看作是由大量微观粒子组成的,这些微观粒子在空间网格上按照一定的规则进行碰撞和传输。
通过对微观粒子的运动状态进行统计,可以得到流体的宏观性质,如密度、速度等。
LBM的核心是格子玻尔兹曼方程(Lattice Boltzmann Equation,简称LBE),它描述了微观粒子在空间网格上的运动规律。
在LBM中,流体的宏观性质由分布函数来描述,分布函数是表示在某一时刻某一空间点上流体微观粒子的分布情况。
在每个时间步内,分布函数按照一定的规则进行碰撞和传输,通过迭代计算可以得到流体在空间网格上的演化过程。
LBM的计算过程可以并行化,因此在计算效率上具有明显的优势。
LBM的另一个优点是它对复杂边界条件的处理能力强。
由于LBM是基于离散网格的方法,因此可以比较容易地处理复杂的边界条件,如曲面边界、移动边界等。
这使得LBM在模拟复杂流体系统时具有一定的优势。
除此之外,LBM还有一些其他的优点,如对多相流、多孔介质流动等复杂流体现象的模拟能力强,对于非稳态流动和湍流流动的模拟也有一定的优势。
总之,格子玻尔兹曼方法作为一种新兴的计算流体力学方法,具有诸多优点,逐渐得到了流体力学领域的广泛关注和应用。
随着计算机硬件性能的不断提升,LBM的应用前景将更加广阔,相信它会在流体力学领域发挥越来越重要的作用。
KdV方程的格子Boltzmann模型求解

第46卷 第1期华北理工大学学报(自然科学版)V o l .46 N o .12024年01月J o u r n a l o fN o r t hC h i n aU n i v e r s i t y o f S c i e n c e a n dT e c h n o l o g y (N a t u r a l S c i e n c eE d i t i o n )J a n .2024收稿日期:2023-03-17 修回日期:2023-12-18基金项目:国家自然科学基金项目(32070669)㊂ 第一作者:陈梦涵,女,硕士研究生㊂研究方向:应用数理统计㊂E -m a i l :2094953965@q q .c o m. 通讯作者:王希胤,男,博士,教授㊂研究方向:生物信息学㊂E -m a i l :w a n g x i y i n @v i p.s i n a .c o m. D O I :10.3969/j.i s s n .2095-2716.2024.01.013文章编号:2095-2716(2024)01-0103-08K d V 方程的格子B o l t z m a n n 模型求解陈梦涵,王希胤,李金(华北理工大学理学院,河北唐山063210)关键词:K d V 方程;D 1Q 5模型;格子B o l t z m a n n 方法摘 要:浅水波模型被广泛地用于模拟水波传播的动力学行为㊂很多问题,如强非线性问题㊁非平衡问题㊁实际应用中发生的问题等,使得传统的理论研究手段通常无能为力㊂文章首先给出了格子B o l t z m a n n 方法(L B M )的基本理论,然后利用经典的一维五速度(D 1Q 5)的离散速度模型,给出K o r t e w e g -d eV r i e s (K d V )方程中含有修正项的格子B o l t z m a n n (L B )模型推导公式,最后进行数值模拟,将K d V 方程的精确解和模拟解进行比较,然后验证修正模型的精确性㊂实验结果表明,用格子B o l t z m a n n 方法对K d V 方程进行求解,其模拟解和精确解吻合度较高㊂中图分类号:O 212.1 文献标识码:A波动现象是自然界最普遍的现象之一,对于水波的研究也一直是科学和工程研究领域的重要课题㊂尽管波动问题所涉及的领域不同,但是描述波动现象的方程却是相同的㊂由于水波千姿百态,用肉眼就可以观察到,因此很早就引起了人们的注意,可以说是人们最为熟悉的一种波㊂波动是物质运动的重要形式,广泛存在于自然界㊂波动中被传递的物理量的扰动和振动有多种形式,例如,弦线中的波㊁空气或固体中的声波㊁水波㊁电磁波,等等㊂为了更加具体地研究各种波动,就产生了各种形式的波动方程,因此,浅水波方程也成为重要的研究对象之一㊂而K o r t e w e g -d eV r i e s (K d V )方程是1895年由荷兰数学家科特韦格(K o r t e w e g )和德弗里斯(d eV r i e s )在研究浅水中小振幅长波运动时共同发现的一种单向运动浅水波偏微分方程㊂在求解偏微分方程地过程中,我们经常用到的数值计算方法有:有限元法(F i n i t eE l e m e n tM e t h o d s),有限差分法(F i n i t eD i f e r e n c eM e t h o d s ),有限体积法(F i n i t eV o l u m e M e t h o d s )和格子B o l z m a n n (L B M )方法等㊂其中,有限差分法虽然相对其他三种方法而言简便易行,而且有丰富多样的离散方法,但是它在求解问题时对求解区域的适应性比较差㊂有限元法虽然采用的网格剖分更加灵活,从一定程度上讲对求解区域具有更强的适应性,但是它在求解间断问题时会受到很大的限制,达不到有限差分法的效果㊂而有限体积法可以被视为是上述两种方法的结合,虽然能够充分利用有限元网格灵活性和克服差分法对网格适应性差的缺陷,但是数值实验较难进行[1]㊂作为一种新兴的数值模拟方法,L B M 基于B o l t z m a n n 方程的离散,是一种自下而上的求解方法㊂它描述了微观粒子的碰撞和迁移,利用分布函数(一种概率密度分布函数)来确定粒子的分布,即分布函数描述了流体的宏观运动㊂近年来,由于L B M 具有计算简便㊁良好的并行性㊁处理不规则的复杂边界容易且对于源项的考虑简单等诸多优势,已经自然而然地发展成为了求解浅水波方程的一种新方法[2]㊂在以往的研究中,英国利物浦大学的教授Z h o u [3]较为全面地阐述了浅水波方程的L B M 理论,包括外力的不同处理格式㊁湍流模型的构造㊁多种边界条件的处理方法以及对于许多经典浅水波问题的验证㊂中山大学环境科学与工程学院的L i 和H u a n g [4]进行了对流-扩散方程与浅水波方程耦合的研究,并采用L B 的多松弛模型和双松弛模型分别对流场和污染物场进行了模拟㊂文献[5]提出了一类粘性浅水方程的晶格B o l t z m a n n (L B )模型,该模型采用源项的二阶矩来恢复控制方程中的粘性,并消除C h a p m a n -E n s k o g 分析过程中产生的附加误差㊂文献[6]建立了一种适用于浅水方程的晶格玻尔兹曼模型(L A B S W E ),它用源项如床面坡度,床面摩擦力来求解方程㊂通过求解定常和非定常流动问题,验证了该模型的有效性㊂鉴于以上背景,文章首先给出了格子B o l t z m a n n 方法(L B M )的基本理论,然后利用经典的一维五速度(D 1Q 5)的离散速度模型,给出K d V 方程中含有修正项的格子B o l t z m a n n (L B )模型推导公式,最后进行数值模拟,将K d V 方程的精确解和模拟解进行比较,然后验证修正模型的精确性㊂实验结果表明,用格子B o l t z m a n n 方法对K d V 方程进行求解,其模拟解和精确解吻合度较高㊂1方法格子玻尔兹曼方法(L a t t i c eB o l t z m a n n M e t h o d)是一种基于微观介观的流体力学计算方法,适用于二维或三维流体流动问题的模拟[7]㊂图1是格子玻尔兹曼方法的流程图:其主要思想是离散化流体的分布函数,通过对分布函数的演化来模拟流体的运动㊂近年来,L B M 由于计算简单㊁并行性好㊁易于处理复杂不规则的边界及能简单方便地考虑源项等优势,已经发展成为求解浅水波方程(S W E s )的一种新方法㊂下面是格子玻尔兹曼方法的求解步骤:图1 L B M 方法流程图(1)确定格子和速度模型:首先需要确定流场的离散化格点和速度模型㊂通常情况下,将流体分成若干个小区域,每个小区域都对应一个格点,格点上有一组离散的速度向量㊂(2)定义分布函数:为了描述格点上流体的状态,引入一个分布函数g ,用来表示在每个格点上,每个速度方向上的粒子数密度㊂它是时间和位置的函数,通常用离散的速度和离散的时间步长表示㊂(3)离散B o l t z m a n n 方程:基于B o l t z m a n n 方程,对分布函数进行离散化,得到离散化的B o l t z m a n n 方程,它描述了分布函数的演化过程㊂在格子玻尔兹曼方法中,B o l t z m a n n 方程可以看成是一个简单的微分方程,其左侧是分布函数的时间导数,右侧是一个碰撞项和一个弛豫项,用于描述粒子之间的相互作用和粒子与流体之间的相互作用㊂401 华北理工大学学报(自然科学版) 第46卷(4)离散碰撞项和弛豫项:将碰撞项和弛豫项进行离散化,得到离散化的碰撞算子和弛豫算子㊂碰撞算子用于描述粒子之间的相互作用,而弛豫算子用于描述粒子与流体之间的相互作用㊂(5)迭代求解:通过迭代求解离散化的B o l t z m a n n 方程,计算出每个格点上的分布函数,从而得到流场的速度场和密度场㊂(6)计算宏观量:根据格点上的分布函数,可以计算出宏观量,如速度㊁密度㊁压力等㊂(7)处理边界条件:对于边界处的格点,需要根据具体的物理问题设置边界条件㊂(8)模拟结束:当达到预设的模拟时间或达到收敛条件时,模拟结束㊂2模型简介2.1 离散速度模型L B M 中一维离散速度模型最常见的是D 1Q 3模型和D 1Q 5模型[8],具体如下:(1)对于D 1Q 3模型(见图2),模型参数如下: c ң=c 01-1[], c s =c3, ωi =2/3,c i2=01/6,c i2=c 2{(1)图2 D 1Q 3离散速度模型其中,ωi 为权重系数,c =Δx /Δt 为粒子迁移速度,c s 是与当地声速相关的量㊂(2)对于D 1Q 5模型(见图3),模型参数如下: c ң=c 01-12-2[],c s =c (2) ω0=12,ω1=ω2=16,ω3=ω4=112(3)图3 D 1Q 5离散速度模型其中,ωi 为权重系数,c =Δx /Δt 为粒子迁移速度,c s 是与当地声速相关的量㊂2.2 K d V 方程将L B 模型应用于K d V 方程中,需要将K d V 方程离散化成网格上的方程组,然后通过L B 模型求解这个方程组㊂具体来说,L B 模型中的速度分布函数被定义为格点上的波高,通过计算速度分布函数在不同时间和空间的演化来模拟K d V 方程的行为㊂与传统的有限差分法和有限元法相比,L B 模型具有计算效率高㊁适合并行计算等优点,因此在模拟非线性波等问题时得到了广泛应用㊂考虑非线性偏微分方程一般形式[9]: ∂u ∂t +αu ∂u ∂x +βu n ∂2u ∂x 2-γ∂2u ∂x 2+δ∂3u ∂x3=0(4)其中,u =u x ,t ()是物质在空间x 处和时刻t 时的密度,α,β,γ,δ为参数㊂当β=0,γ=0时,方程(4)化为K d V 方程:∂u ∂t +αu ∂u ∂x +δ∂3u∂x3=0(5)501 第1期 陈梦涵,等:K d V 方程的格子B o l t z m a n n 模型求解3模型推导采用D 1Q 5模型给出K d V 方程含有修正项的L B 模型推导,给出的演化方程为: g i x ң+c ңi Δt ,t +Δt ()-g i x ң,t ()=-1τg i x ң,t ()-g e q i x ң,t ()()+Δt h i x ң,t ()(6)其中,τ表示松弛时间,c ңi 为沿着流动方向的单位速度矢量,gi x ң,t ()表示在t 时刻,位于x ң处沿着离散速度方向c ңi 运动的粒子分布函数,Δt h i x ң,t ()为修正项㊂根据参考文献[10-12],将修正函数h i x ң,t ()和平衡态分布函数g e qi x ң,t ()分别定义为: h i x ң,t ()=λ1,i u 2+λ2,iu n +1,i =0~4(7) g e qix ң,t ()=ρ1u (8)其中,λ1,i ,λ2,i 和p i 为调整参数㊂宏观变量u x ,t ()定义为[13-15]: u =ðigi (9)为了得到稳定的宏观变量u ,假设分布函数g i 也处于平衡状态,且有: u =ðig e qi =ði g 0()i =ðigi (10)由(10)可得, ði gn ()i=0,n ⩾1(11)为了能够恢复到宏观方程,平衡态分布函数g e q i 和修正函数需要满足下面几个约束条件:ðic i g e qi()=0,ðic 2i g e q i ()=c 2λu ,ðic3i g e qi ()=c 3ηu (12) ði hi=0,ði c i h i ()=c λ1u 2+c λ2u n +1,ðic2i h i ()=0(13)使用多尺度分析将方程恢复到宏观方程,引入1个离散的时间尺度和3个连续的时间尺度,其具体表达形式为:t 0,t 1=εt ,t 2=ε2t ,t 3=ε3t (14)对分布函数g i 和时间导数进行Ch a p m a n -E n s k o g 展开,可得: g i =g 0()i +εg1()i +ε2g 2()i +ε3g 3()i +O ε4()(15)∂∂t =∂∂t 0+ε∂∂t 1+ε2∂∂t 2+ε3∂∂t 3+O ε4()(16)其中,ε表示任意小的参数,在宏观方程的推导过程中,不妨假设ε=Δt ,将(6)式的左边对时间和空间进行泰勒展开,并保留Οε4()项,可得 ε∂g i ∂t +ε∂∂x c i g i ()+12ε2∂2g i ∂t 2+ε2∂2∂t ∂x c i g i ()+12ε2∂2∂x 2c 2i g i ()+16ε3∂∂t +c i ∂∂x æèçöø÷3g i =-1τg i -g e qi ()+εh i +Οε4()(17)将(15)和(16)代入(17)式中得,并对比左右两边可得ε的同阶项:可以得出,O ε()系数: g e qi =g0()i (18)O ε1()系数:601 华北理工大学学报(自然科学版) 第46卷∂∂x c i g e qi ()=-1τg 1()i +h i (19)O ε2()系数: ∂g e q i ∂t 1+∂∂x c i g 1()i ()+12∂2∂x 2c 2i g e qi ()=-1τg 2()i (20)O ε3()系数: ∂g e q i ∂t 2+∂g 1()i ∂t 1+∂∂x c i g 2()i ()+∂2∂x ∂t 1c i g e qi ()+12∂2∂x 2c 2i g 1()i ()+16∂3∂x 3c 3i g 1()i ()=-1τg 3()i (21)结合约束条件(12)和(13),将方程(19)两边分别乘以c i 和c 2i 后并对i 求和得:ði c i g 1()i =τði c i h i -τ∂∂x ðic 2i ge q ()i ()=c τλ1u 2+λ2u n +1()-c 2τλ∂u ∂x (22) ði c 2i g 1()i =τðic 2i h i -τ∂∂x ði c 3i g e q ()i ()=-c 3τη∂u ∂x (23)同理,结合约束条件(12)和(13),将方程(20)式两边分别乘以c i 后并对i 求和,得出:ðic i g 2()i=τ∂∂t 1ðic 2i g e q i ()+∂∂x ði c 2i g 1()i ()+12∂2∂x 2ði c 3i g e q ()i ()æèçöø÷ =c 3τ2η∂2u ∂x 2-12c 3τη∂2u ∂x 2=c 3ητ2-12τæèçöø÷∂2u ∂x2(24)结合(7)(8)(9)和(19),将方程(16)两边对i 求和,得出: ∂u ∂t 1+2c τλ1u ∂u ∂x +n +1()c τλ2u n ∂u ∂x +c 2λ12-τæèçöø÷∂2u ∂x 2=0(25)同理,结合(10)(11)(12)和(22)~(24),将方程(20)两边对i 求和,得出: ∂u ∂t 2+c 3ητ2-τ-16æèçöø÷∂3u ∂x 3=0(26)将3.19()ˑε+3.20()ˑε2,可得: ∂u ∂t 1+2c τλ1εu ∂u ∂x +n +1()c τλ2εu n ∂u ∂x +c 2λε12-τæèçöø÷∂2u ∂x 2+c 3ηε2τ2-τ-16æèçöø÷∂3u ∂x 3=0(27)将方程(27)和(4)对比可得: α=2c τλ1ε,β=n +1()c τλ2ε(28) γ=c 2λ12-τæèçöø÷ε,δ=c 3ητ2-τ-16æèçöø÷ε2(29)其中,τ=12+112+δε2ηc 3,λ=γεc 2τ-1/2(),λ1=α2τεc ,λ2=βn +1()τεc 可以得出,方程(27)就是一维K D V 方程的L B 模型㊂由方程(7)和(13)式,得出修正函数h i 为: h 0=-12λ1u 2-12λ2u n +1h 1=h 2=13λ1u 2+13λ2u n +1h 3=16λ1u 2+16λ2u n +1h 4=-13λ1u 2-13λ2u n +1ìîíïïïïïïïïïï(30)701 第1期 陈梦涵,等:K d V 方程的格子B o l t z m a n n 模型求解联立方程(10)和(12)式,可得平衡态分布函数f e q i 为:g e q 0=η2+6ρ4-λ+1æèçöø÷u =ρ0u g e q 1=λ-η2-4ρ4æèçöø÷u =ρ1u g e q2=λ2-η6-4ρ4æèçöø÷u =ρ2u g e q 3=η6+ρ4æèçöø÷u =ρ3u g e q 4=ρ4u ìîíïïïïïïïïïïïï(31)4数值模拟将L B 模型应用于K d V 方程中,需要将K d V 方程离散化成网格上的方程组,然后通过L B 模型求解这个方程组㊂具体来说,L B 模型中的速度分布函数被定义为格点上的波高,通过计算速度分布函数在不同时间和空间的演化来模拟K d V 方程的行为㊂与传统的有限差分法和有限元法相比,L B 模型具有计算效率高㊁适合并行计算等优点,因此在模拟非线性波等问题时得到了广泛应用㊂(1)考虑如下的K d V 方程:∂u ∂t +6u ∂u ∂x +∂3u∂x3=0设置参数如下:Δx =0.001,Δt =0.00001,边界条件u 0,t ()=u 255128π,t æèçöø÷=0,t >0,初始条件u x ,0()=s i n x ,x ɪ0,255128πéëêêùûúú,该方程的精确解为: u x ,t ()=c 212s e c h 2c 212x -c 21t ()+l n c 2-l n c 3éëêêùûúú取c 1=2c ,c 2=c 3=1,得出该方程的精确解为: u x ,t ()=2c 2s e c h 2c x -4c 2t ()[],x ɪ0,255128πéëêêùûúú图4 t =0.01,L B 模拟解和精确解对比 图5 t =0.25,L B 模拟解和精确解对比模拟结果如图4㊁图5所示㊂图4和图5分别给出了t =0.01和t =0.25时刻的L B 模拟解和解析解的对比图,从图中可以看出:在t =0.25之前,模拟解和解析解吻合的程度较高,但是随之时间的推移,模拟解与解析解存在一定的偏离,这主要是原因有扰动项O ε4(),它在一定程度会对孤子高度㊁速度以及形状有影801 华北理工大学学报(自然科学版) 第46卷响,且当t >0.25时,方程的模拟解和精确解差别较大㊂(2)考虑如下的K d V 方程:∂u ∂t +αu ∂u ∂x +δ∂3u∂x3=0其中,u =u x ,t (),u 为波动地振幅,x 为波横向传播的位移,t 为时间㊂初始条件为:u x ,0()=3A s e c h 2B x +C (),x ɪ0,2[]边界条件为:u 0,t ()=u 2,t ()=0,t >0解析解为:u x ,t ()=3A s e c h 2B x -D t +C (),x ɪ0,2[]其中,B =12αAδ,D =αA B .设置参数如下:Δx =0.001,Δt =5ˑ10-4,α=1,δ=4.84ˑ10-4,模拟结果如图(6)和图(7)所示. 图6 t =0.0005时,模拟解与解析解对比 图7 t =0.0020时,模拟解与解析解对比图6和图7分别给出了t =0.0005和t =0.002时刻的L B 模拟解和解析解的对比图,通过对该方程的模拟结果进行分析,发现在t =0.002之前,模拟解和解析解吻合的程度较高,但是随着时间的推移,模拟解与解析解存在一定的偏离,主要的原因是含有扰动项O ε4(),当然也可能是该类波的传播速度极快,长时间模拟就会产生偏差㊂表1给出了该方程在不同时刻的误差㊂表1 方程在不同模拟时刻的误差比较时刻/tt =0.0005t =0.0010t =0.0015t =0.0020L ɕ3.2371ˑ10-43.8721ˑ10-42.0518ˑ10-33.6037ˑ10-3L 27.4783ˑ10-55.8732ˑ10-56.4976ˑ10-47.1335ˑ10-4G R E3.8516ˑ10-44.7039ˑ10-42.5450ˑ10-34.6530ˑ10-3 从表1可以看出,在L B 模型下得到的模拟解和解析解非常逼近,无论是L ɕ误差,还是均方根误差L 2和整体相对误差G R E ,其两者之间的误差数量级都达到了,说明该数值结果是比较理想的㊂5结论现在,微分方程无处不在,各个科学领域的研究都伴随着微分方程模型㊂由于实际生活中的微分方程模型形式日趋复杂,为了与实际问题相匹配,微分方程解的形式越来越多样化㊂本文对两个特殊的K D V 方程,利用格子B o l t z m a n n 模型求解并与其精确解进行比较,得出使用格子B o l t z m a n n 方法对非线性偏微分方程求解取得了较好的效果㊂在未来的工作中,将尝试继续改进格子B o l t z m a n n 模型,并对更加复杂的偏微分901 第1期 陈梦涵,等:K d V 方程的格子B o l t z m a n n 模型求解011华北理工大学学报(自然科学版)第46卷方程或者浅水波方程进行模拟㊂希望本文可以为其他学者在求解偏微分方程方面的研究工作提供一定的参考价值㊂参考文献:[1]张海军.求解浅水波方程的熵稳定格式研究[D];西安:长安大学,2018.[2]陈文文,张文欢,汪一航,等.浅水波方程的一类改进的格子B o l t z m a n n模型[J].宁波大学学报(理工版),2020,33(01):72-79.[3] R I C K,S A L MO N.T h e l a t t i c eB o l t z m a n nm e t h o d a s ab a s i s f o r o c e a n c i r c u l a t i o nm o d e l i n g[J].J o u r n a l o fM a r i n eR e s e a r c h,1999,57(3).[4]冯士德,赵颖,茑原道久,等.旋转流场中的格子波耳兹曼模型[J].地球物理学报,2002,45(2):170-175.[5] Y U L,Z C AC,X G D,e t a l.A l a t t i c eB o l t z m a n nm o d e l f o r t h e v i s c o u s s h a l l o ww a t e r e q u a t i o n sw i t h s o u r c e t e r m s[J].J o u r n a l o fH y-d r o l o g y,2021.[6] Z H O UJG.L a t t i c eB o l t z m a n nm o d e l f o r t h e s h a l l o w w a t e r e q u a t i o n s[J].C o m p u t e rM e t h o d s i nA p p l i e d M e c h a n i c s a n dE n g i n e e r i n g,2002,191(32):3527-39.[7]张宗宁.基于格子B o l t z m a n n方法求解若干非线性偏微分方程[D];银川:北方民族大学,2022.[8]戴厚平,郑洲顺,段丹丹.一类偏微分方程的格子B o l t z m a n n模型[J].计算机工程与应用,2016,52(3):21-26.[9]王慧敏.非线性偏微分方程中孤波解的格子B o l t z m a n n模拟[D];长春:吉林大学,2014.[10]何郁波,林晓艳,董晓亮.应用格子B o l t z m a n n模型模拟一类二维偏微分方程[J].物理学报,2013,62(19):290-296.[11]乐励华,高云,刘唐伟.偏微分方程求解的一种新颖方法--格子B o l t z m a n n模型[J].大学数学,2011,27(03):75-82.[12]张春泽,程永光,李勇昌.二维浅水波方程格子B o l t z m a n n算法的G P U并行实现[J].水动力学研究与进展A辑,2011,26(2):194-200.[13]赫万恒,钱跃竑.浅水波方程的格子B o l t z m a n n模拟[C].中国力学学会北方七省市区第十三届学术大会论文集.郑州.2010:42-45.[14]赖惠林,马昌凤.非线性偏微分方程的高阶格子B G K模型[J].中国科学(G辑:物理学力学天文学),2009,39(07):913-922.[15]施卫平,胡守信,阎广武.用格子B o l t z m a n n方程模拟浅水波问题[J].力学学报,1997,(05):7-11.S o l u t i o n M e t h o d f o rK D VE q u a t i o nb y L a t t i c eB o l t z m a n n M o d e lC H E N M e n g-h a n,WA N G X i-y i n,L I J i n(C o l l e g e o f S c i e n c e,N o r t hC h i n aU n i v e r s i t y o f S c i e n c e a n dT e c h n o l o g y,T a n g s h a nH e b e i063210,C h i n a)K e y w o r d s:K D Ve q u a t i o n;D1Q5m o d e l;l a t t i c eB o l t z m a n nm e t h o dA b s t r a c t:S h a l l o w w a t e rw a v em o d e l i sw i d e l y u s e d t o s i m u l a t e t h e d y n a m i c b e h a v i o r o fw a t e rw a v e p r o p a-g a t i o n i n o c e a n a n d a t m o s p h e r e f i e l d.M a n y p r o b l e m s,s u c h a s s t r o n g n o n l i n e a r p r o b l e m s,n o n-e q u i l i b r i u m p r o b l e m s,p r o b l e m s i n p r a c t i c a l a p p l i c a t i o n s,m a k e t h e t r a d i t i o n a l t h e o r e t i c a l r e s e a r c hm e t h o d s a r e u s u a l l y p o w e r l e s s.I n t h i s p a p e r,t h e b a s i c t h e o r y o f l a t t i c eB o l t z m a n nm e t h o d(L B M)w a s f i r s t l y g i v e n.T h e n,t h ef o r m u l a o f l a t t i c eB o l t z m a n n(L B)m o d e lw i t hc o r r e c t i o n t e r mi nK o r t e w e g-d eV r i e s(K d V)e q u a t i o nw a sg i v e nb y u s i n g t h e c l a s s i c a l o n e-d i m e n s i o n a l f i v e-v e l o c i t y(D1Q5)d i s c r e t ev e l o c i t y m o d e l.F i n a l l y,t h e a c-c u r a t e s o l u t i o no f t h eK d Ve q u a t i o nw a s c o m p a r e dw i t ht h e s i m u l a t e ds o l u t i o n,a n d t h e nt h ea c c u r a c y o f t h em o d i f i e dm o d e l i s v e r i f i e d.T h e e x p e r i m e n t a l r e s u l t s s h o wt h a t t h e l a t t i c eB o l t z m a n nm e t h o d i s u s e d t o s o l v e t h eK d Ve q u a t i o n,a n d t h e s i m u l a t i o n s o l u t i o n i s i n g o o d a g r e e m e n tw i t h t h e e x a c t s o l u t i o n.。
基于格子波尔兹曼方法的回热器数值模拟

基于格子波尔兹曼方法的回热器数值模拟夏宇栋;陈曦;马诗旻;张华【摘要】Based on Lattice Boltzmann Method ( LBM ) , the flow field in the microstructure of mesh screen and etched foil regenerators were simulated. The velocity and pressure drop were obtained by LBM. The simulation results show that velocity field in etched foil regenerator is better distributed than that of mesh screen. And the etched foil regenerator has less resistance coefficient than that of mesh screen regenerator. The simulation results basically agree with the experiments.%利用格子玻尔兹曼方法,直接对蚀刻薄片和层叠丝网回热器的微观结构流场进行了模拟.得到了两种回热器填料的微观流场和两端的压差.模拟结果显示,当回热器的直径、水力直径和填充率相近情况下,不同流速下蚀刻薄片卷裹式回热器的稳态阻力系数均比层叠丝网回热器小.稳态阻力系数的模拟变化趋势与实验一致.【期刊名称】《低温工程》【年(卷),期】2012(000)005【总页数】5页(P41-45)【关键词】格子波尔兹曼方法;回热器;流阻系数;数值模拟【作者】夏宇栋;陈曦;马诗旻;张华【作者单位】上海理工大学能源与动力学院上海200093;上海理工大学能源与动力学院上海200093;上海理工大学能源与动力学院上海200093;上海理工大学能源与动力学院上海200093【正文语种】中文【中图分类】TB651回热器是低温制冷机的关键部件,也是影响制冷机性能的最重要因素。
格子Boltzmann方法三种边界格式的对比分析

格子Boltzmann方法三种边界格式的对比分析刘连国;杨帆;王宏光【摘要】采用格子Boltzmann方法(Lattice Boltzmann Method- LBM)对二维顶盖驱动方腔流动进行数值模拟.在计算中分别使用半步长反弹、非平衡反弹、以及非平衡外推三种边界处理格式,并得到了不同格式对应的流线分布,流函数最小值、涡心坐标、几何中心线速度分布等.通过将所得结果与基准解进行比较,就三种边界格式的计算效率,计算精度、以及计算稳定性等方面进行了讨论和分析,为LBM计算中边界格式的选择提供了有益的参考.【期刊名称】《机械研究与应用》【年(卷),期】2012(000)001【总页数】5页(P18-22)【关键词】格子Boltzmann方法;边界处理格式;半步长反弹格式;非平衡反弹格式;非平衡外推格式【作者】刘连国;杨帆;王宏光【作者单位】上海理工大学能源与动力工程学院,上海200093;上海理工大学能源与动力工程学院,上海200093;上海理工大学能源与动力工程学院,上海200093【正文语种】中文【中图分类】O357.11 引言格子Boltzmann方法(LBM)是近年来迅速发展的一种新型数值计算方法。
边界条件的处理是LBM实施中一项非常关键的内容。
实际计算表明:选取不同的边界条件会对数值计算的精度、稳定性以及效率产生很大影响。
作为LBM的一个基本问题,边界条件的处理一直是流体力学一个重要的研究方面。
根据边界条件的类型,可将之分为两类:压力边界和速度边界[1],其中的速度边界又可细分为:平直边界和曲面边界。
笔者从经典的流体力学问题二维顶盖方腔流模拟入手,对三种平直边界格式进行对比和分析,为LBM计算中边界格式的选择提供了有益的参考。
2 二维九点格子Boltzmann模型目前最常用的格子Boltzmann模型为LBGK模型,通过引入“单一弛豫时间”来简化Boltzmann方程中碰撞项的计算[2]。
九点格子LBGK模型的演化方程为:式中:(x,t)是在t时刻、x处的平衡态分布函数;τ为单一弛豫时间因子;eα为网格点各方向上的粒子速度。
格子Boltzmann方法三种边界格式的对比分析

中分别使 用半步长反 弹、 非平衡反 弹、 以及 非平衡 外推 三种边界处理格 式 , 并得 到 了不同格 式对应 的流 线分
布 , 函数最 小值、 流 涡心 坐标、 几何 中心线速度分布等。通过将所得结果与基 准解进行 比较 , 三种边界格 式 就
的计算效率 , 计算精 度、 以及计算稳 定性等方 面进行 了讨论和分析 , L M 计算 中边界格 式 的选择 提供 了有 为 B
Ab t a t n t i p p r wo d me so a i - r e a i o s s lae y lt c l ma n meh d T r e b u d r sr c :I h s a e ,t - i n in ll d i n c vt f w i i d v yl mu td b at e Bot n to . h e o n a y i z
益 的参考。 关键 词 : 格子 B hman方法; oz n 边界处理格 式; 半步长反弹格式 ; 非平衡 反弹格 式; 非平衡 外推格式 中图分类号 : 77 1 0 5 . 文献标识码 : A 文章编号 :07 4 1 (0 2 0 — 0 8 0 10 - 4 4 2 1 ) 1 0 1 — 5 T ei mp r t ea ay i ft r eb u d r c e a e i It e l ma n meh d h o a ai n ls o e o n a y s h meb sd Ol a f eBot : v s h i z n to
b u d r c e e e t g i BM i l t n o n ay s h me s l ci L n n s muai . o
Ke r s:lt c otma n meh d y wo d at e B l n t o ;b u d r c e ;h l — y b u c —b c c e ;n n q i b i m o c —b c i z o n a y s h me a wa o n e a k s h me o e u l ru b u e a k f i n s h me;n n q i b i m xr p lt n s h me ce o e ul ru e  ̄ o o ua o o iee tsh me ae dsus d tbH fcmp tt n frdf rn c e r ic se .whc r x etd t mdd au be rfrn e t i f ih ae ep ce o p e v a l eee c o l
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
1. 引言
随着数值计算方法成为研究科学与工程问题的重要手段,对计算效率的要求也越来越高。为了确保 在尽可能短的时间内完成一项计算任务,在单机中我们主要有两个选择,第一个选择是设计一个更好的 算法,可以用更少的步骤达到相同的结果。第二种方法是在指令级别执行优化,该方法是在需要较多时 间完成的指令的位置上使用一个或多个更快的指令。
收稿日期:2020年7月22日;录用日期:2020年8月5日;发布日期:2020年8月12日
*通讯作者。
文章引用: 李阳贵, 吕莹, 梁大成, 黄笑冲. 格子 Boltzmann 数值计算的性能优化[J]. 软件工程与应用, 2020, 9(4): 272-277. DOI: 10.12677/sea.2020.94030
第二个优化是注意到在迁移期间为每个节点存储的分布的值是在碰撞步期间读取的那些值。因此,
四个函数 stream、computeRhoU、collide 和 save 可以组合成一个函数 stream_computeRhoU_collide_save, 该函数访问内存的频率显著降低。
格子 Boltzmann 方法是模拟微观模型的数值计算方法。格子 Boltzmann 方法除了被应用于一般的流 体力学问题之外,还在湍流[1]、多孔介质流[2] [3]、粒子悬浮流[4]、磁流体力学[5]、多相流[6]等相关领 域也取得了比较成功的应用。但格子 Boltzmann 方法的计算量较大,为提高计算效率,本文将探讨泰勒– 格林涡流(Taylor-Green Vortex)格子 Boltzmann 数值计算的性能优化。
e−t td
(13)
DOI: 10.12677/sea.2020.94030
274
软件工程与应用
李阳贵 等
( ) td
= υ
1
k
2 x
+
k
2 y
(14)
( ) p ( x,t)
= p0 − ρ u402 kkxy cos (2kx x)
+
kx ky
cos
2ky y
(15)
其中 u0 为初始速度大小, kx, y = 2π x, y 为波矢量分量, td 为涡流衰减时间,平均压力 p0 可以是任意的, 计算区域为 x × y = 640 × 640 (格子单位)。计算采用标准平衡分布,平衡分布初始化采用 Mei [9]的初始 化方案。泰勒–格林涡流的速度 ux 和 uy 如图 1 所示。
李阳贵 等
摘要
格子Boltzmann方法是科学与工程计算的重要模型与数值方法,该方法能够从底层描述和计算细尺度特 征,同时计算量也较大。为了提高计算效率,本文主要从指令级别对格子Boltzmann数值计算进行性能 优化。数值实验结果表明,优化后计算性能显著提高。
关键词
格子Boltzmann方法,性能优化,泰勒–格林涡流
Copyright © 2020 by author(s) and Hans Publishers Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). /licenses/by/4.0/
ρ ( x,t ) = ∑ fi ( x,t ) , ρu ( x,t ) = ∑ ci fi ( x,t )
(11)
i
i
(IV) 执行碰撞(collide)
( ) fi∗
(
x,
t
)
= fi ( x,
t
)
−
∆t τ
fi ( x,t ) − fieq ( x,t )
(12)
(V) 保存数据(save) (每 saveN 步保存一次) (VI) 重复(II)~(V)直到满足终止条件。
(3)
同时,粒子受到碰撞算符 Ωi ( x,t ) 的影响,该算符通过在每个位置的分布中重新分配粒子来模拟粒子
碰撞。虽然有许多不同的碰撞算符可用,但最常用的一个算符是 Bhatnagar-Gross-Krook (BGK)运算符[7]:
DOI: 10.12677/sea.2020.94030
273
软件工程与应用
∂t fi + ci ⋅ ∇fi = Ωi ( f )
(9)
的时间分裂方法[8],具体步骤如下:
(I) 初始化分布函数 fi ( x, 0)
(II) 执行迁移(stream)
fi ( x + ci∆t,t + ∆t ) = fi∗ ( x,t )
(10)
(III) 计算宏观量(computeRhoU)
(7)
第二部分是迁移(或传播):
fi ( x + ci∆t,t + ∆t ) = fi∗ ( x,t )
(8)
碰撞只是一个代数局部运算,首先要计算密度,找到平衡分布后的宏观速度,碰撞后,将得到相邻 节点的分布 fi∗ ,当碰撞和迁移两个操作完成时,经过一个时间步,再重复这些操作。
格子 Boltzmann 数值计算的实现过程可分为碰撞和迁移两个子过程,通常格子 Boltzmann 方法的程 序结构有两种形式:碰撞–迁移结构和迁移–碰撞结构。本文程序设计采用迁移–碰撞结构,这种程序 结构可以看作是求解离散速度方程
3. 计算与性能优化
本文计算泰勒–格林涡流问题,速度和压力的初始状态解析地设定,泰勒–格林涡流在 x × y 区域 内是非定常的全周期流动,其速度场和压力场在二维空间上表示为
(( )) u ( x,t) = u0
−
ky kx
kx cos (kx x)sin ky y ky sin (kx x) cos ky y
2. 格子 Boltzmann 模型
格子 Boltzmann 方法的一个主要优点是基于 Boltzmann 方程而不是连续方程和动量方程,它的实现
比传统方法简单。格子 Boltzmann 方法的基本量是离散速度分布函数 fi ( x,t ) ,通常称为粒子分布。可以
通过此分布函数来求解质量密度和动量密度:
ρ ( x,t ) = ∑ fi ( x,t )
(1)
i
ρu ( x,t ) = ∑ ci fi ( x,t )
(2)
i
fi 函数表示所有的参数变量都是离散的,下标 i 代表一个离散速度集,定义 fi 在空间中为正方形格子各
个方向的分布。速度集通常用 DdQq 表示,d 是速度集覆盖的空间维数,q 是速度方向的个数。最常用的
速度集是 D1Q3、D2Q9、D3Q15、D3Q19 和 D3Q27。本文采用 D2Q9 模型,i 取值范围是 0,1,2,3,4,
5,6,7,8 九个方向。
通过在物理空间、速度空间和时间上离散 Boltzmann 方程,得到了格子 Boltzmann 方程:
fi ( x + ci∆t,t= + ∆t ) fi ( x,t ) + Ωi ( x,t )
李阳贵 等
Ωi
(
f
)
=
−
fi
− τ
f i eq
⋅ ∆t
(4)
它使分布 fi 以由松弛时间τ 决定的速率趋于平衡分布 fieq ,平衡分布为
fieq ( x,t ) =
wi ρ 1+
u ⋅ ci cs2
+
(u ⋅ ci )2
2cs4
−
u
⋅
u
2cs2
(5)
格子 BGK (LBGK)方程(即用 BGK 碰撞算符完全离散的 Boltzmann 方程)可以表示为:
Performance Optimization of Lattice Boltzmann Numerical Calculation
Yanggui Li1,2, Ying Lv3*, Dacheng Liang4, Xiaochong Huang5
1School of Mathematics and Statistics, Lingnan Normal University, Zhanjiang Guangdong 2Institute of Synthetic Biology, Shenzhen Institute of Advanced Technology, Chinese Academy of Sciences, Shenzhen Guangdong 3School of Electromechanical Engineering, Lingnan Normal University, Zhanjiang Guangdong 4Department of Mathematics, Guangdong Preschool Normal College in Maoming, Maoming Guangdong 5Aizhou Senior High School of Zhanjiang, Zhanjiang Guangdong
Software Engineering and Applications 软件工程与应用, 2020, 9(4), 272-277 Published Online August 2020 in Hans. /journal/sea https:///10.12677/sea.2020.94030
Received: Jul. 22nd, 2020; accepted: Aug. 5th, 2020; published: Aug. 12th, 2020
Abstract
Lattice Boltzmann method is an important model and numerical method in scientific and engineering calculation. It can describe and calculate the fine scale characteristics from the bottom layer, and the calculation cost is also large. In order to improve the computational efficiency, this paper mainly optimizes the performance of lattice Boltzmann numerical calculation from the instruction level. The results of numerical experiments show that the computational performance is improved significantly after optimization.