力学学报, 2019, 51(2): 392-404 DOI: 10.6052/0459-1879-18-278

流体力学

浸没边界-简化热格子Boltzmann方法研究及其应用 1)

李桥忠*, 陈木凤††,2), 李游*,3), 牛小东*,4), AdnanKhan*

* 汕头大学工学院机械电子工程系,广东汕头 515063

†† 龙岩学院物理与机电工程学院,福建龙岩 364012

IMMERSED BOUNDARY-SIMPLIFIED THERMAL LATTICE BOLTZMANN METHOD FOR FLUID-STRUCTURE INTERACTION PROBLEM WITH HEAT TRANSFER AND ITS APPLICATION 1)

Li Qiaozhong*, Chen Mufeng††,2), Li You*,3), Niu Xiaodong*,4), Adnan Khan*

* College of Engineering, Shantou University Shantou, Shantou 515063, Guangdong, China

†† College of Electromechanic Engineering, Longyan University, Longyan 364012, Fujian, China

收稿日期: 2018-08-24   接受日期: 2018-11-7   网络出版日期: 2019-03-18

基金资助: 国家自然科学基金资助项目.  11372168
国家自然科学基金资助项目.  11772179

Received: 2018-08-24   Accepted: 2018-11-7   Online: 2019-03-18

作者简介 About authors

2)陈木凤,讲师,主要研究方向:流固耦合传热,磁流体力学.E-mail:903194866@qq.com

3)李游,博士,主要研究方向:计算流体力学.E-mail:16yli10@stu.edu.cn

4)牛小东,教授,主要研究方向:LBM,计算和实验流体力学.E-mail:xdniu@stu.edu.cn

摘要

针对流固耦合传热问题,本文提出了一种基于浸没边界-简化热格子玻尔兹曼方法(immersed boundary method-simplified thermal lattice Boltzmann method,IB-STLBM)的耦合模型.不同于传统的格子玻尔兹曼方法使用分布函数演化流场和温度场,简化热格子玻尔兹曼方法(simplified thermal lattice Boltzmann method,STLBM)的演化过程不需要依赖分布函数,只涉及平衡态分布函数和非平衡态分布函数,能够直接演化宏观量,极大减小了计算过程中所占用的虚拟内存,简化了边界条件的实现方式,同时具有较高的稳定性.传统的浸没边界法对流场的计算采用欧拉网格,对固体边界采用拉格朗日网格,认为固体边界是对流场产生某种体积力.在应用浸没边界法时,汲取介观的思想,把固体的介入看作是对流场的干扰,打破了固体附近流体介观微团颗粒原始的平衡状态,这种干扰可以看作是在耦合边界上产生的一个非平衡项,可用非平衡态分布函数来表示.基于此,在模型中浸没边界法与简化热格子玻尔兹曼方法更紧密联系在一起,更大程度发挥二者的优点,整个计算过程更加简单直观,符合物理特性.通过对热圆柱绕流和内含热颗粒的封闭方腔自然对流问题的模拟以及对其结果的分析,验证了该算法在求解流固耦合传热问题的有效性和可行性.

关键词: 简化热格子Boltzmann方法 ; 浸没边界法 ; 流固耦合传热 ; 自然对流

Abstract

An Immersed boundary-simplified thermal lattice Boltzmann method(IB-STLBM) for fluid-structure interaction problem with heat transfer is developed in this work. In the IB-STLBM, an effective simplified thermal lattice Boltzmann method without the evolution of distribution is used for the intermediate flow field. Different from the stander thermal lattice Boltzmann method, STLBM directly updates the macroscopic variables instead of the distribution functions, which offers several distinct benefits:lower cost in virtual memories, simpler implementation of physical boundary condition and higher numerical stability. In addition, from the mesoscopic view, the existence of solid boundary in the field is considered as an interference of system, which breaks the original equilibrium state of fluid particle, and a non-equilibrium state occurs on the fluid-structure interaction physics boundary. On this basis, in the present IB-STLBM, fluid-structural interaction duo to Immersed boundary appearance in the fluid can be expressed by the non-equilibrium distribution function, which is calculated by the popular non-equilibrium bounce-back boundary condition of the LBM. Hence, the solution procedure of present IB-STLBM can satisfy the non-slip boundary by a simpler way. Numerical experiments for the forced convection over a stationary heated circular cylinder and natural convection in a square cavity with a circle particle are presented to verify the stability, the capability and the flexibility of IB-STLBM for fluid-structure interaction problem with heat transfer. In the case of a stationary heated circular cylinder, quantitative and qualitative comparisons are carried out with previous study. The results of the drag coefficient and the avenge Nusselt numbers on the cylinder are in accordance with the results of previous study. From the case of natural convection in a square cavity with a circle particle, some interesting phenomena can be found. First, the temperature field is clearly stirred by the suspended particle. Second, the temporal trajectories of the particle exhibited regular changes. Third, the particle enhances heat transfer and the average Nusselt numbers periodically oscillate with time.

Keywords: LBM ; IBM ; fluid-structure interaction ; natural convection

PDF (3125KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

李桥忠, 陈木凤, 李游, 牛小东, AdnanKhan. 浸没边界-简化热格子Boltzmann方法研究及其应用 1). 力学学报[J], 2019, 51(2): 392-404 DOI:10.6052/0459-1879-18-278

Li Qiaozhong, Chen Mufeng, Li You, Niu Xiaodong, Adnan Khan. IMMERSED BOUNDARY-SIMPLIFIED THERMAL LATTICE BOLTZMANN METHOD FOR FLUID-STRUCTURE INTERACTION PROBLEM WITH HEAT TRANSFER AND ITS APPLICATION 1). Chinese Journal of Theoretical and Applied Mechanics[J], 2019, 51(2): 392-404 DOI:10.6052/0459-1879-18-278

引 言

流固耦合传热问题是流体力学领域研究者关注的热点,并广泛存在于生活和工业中[1-4].在研究流固耦合传热问题时,如何高效处理流体和固体耦合界面是一个关键节点.在实际的应用中,物体的外形结构往往十分复杂,传统的计算方法大多采用的是贴体网格,需要花费大量的计算用于捕捉物体的形状.虽然合理的网格设计以及生成高质量的网格等技术可以从某种程度上简化边界网格的生成过程,但是对于移动边界仍需要花费大量的计算进行网格重构.为了克服这一障碍,一种简洁高效的流固耦合边界处理方法---浸没边界法(immersed boundary method, IBM)应运而生[5-6].

浸没边界法最初是由Peskin等[5-6]提出的,并成功应用于模拟心脏内的血液流.IBM中,流场的计算采用欧拉网格,固体边界采用拉格朗日网格,其核心思想是把流固边界之间的干扰描述成一种体积力,并采用插值的方法得到流场与固体边界相互影响的信息.IBM不但适用于任意形状的边界,而且在追踪动边界时不需要实时对复杂的边界进行网格重构,极大简化了网格处理方式.格子玻尔兹曼方法(lattice Boltzmann method, LBM)是一种基于介观思想发展而来的方法,通过简单的迁移和碰撞模型实现对流体复杂运动的模拟,具有物理意义清晰、程序易于实施与并行、边界处理简单等优点[7].IBM简洁高效的边界处理方式及其使用欧拉节点来模拟流场等特性,使得IBM与LBM具有天然的亲和性[8-11].Feng等[8-9]首先把IBM与LBM结合起来并成功应用于流固耦合模拟,先后提出惩罚力IB-LBM耦合模型[8]和直接力IB-LBM耦合模型[9].IBM和LBM的耦合模型兼具两种方法的优点,使二者优势互补,LBM发挥了处理流场时的简洁高效特点,同时IBM处理边界的优势弥补LBM在这一方面上的缺陷[10-11].

IB-LBM耦合模型不断改进与发展,并在等温流动领域取得很好成效,但在实际的问题中,流体流动伴随着温度或是能量的转移过程.对热流动问题的研究,推动了浸没边界法与热格子玻尔兹曼方法(thermal lattice Boltzmann method, TLBM)[12-13]耦合模型的发展,并成功应用在颗粒流、自然对流、强制对流和混合对流等问题[14-18]上.IBM和TLBM的耦合模型虽然具有许多优良特性,但是缺陷也同样伴随着TLBM存在.在现有的IB-TLBM耦合模型中,大都是基于耦合双分布函数TLBM模型发展而来.耦合双分布函数TLBM模型采用两个分布函数分别模拟速度场和温度场,在数值模拟中具有较高的稳定性,但是在迁移和碰撞步中都需要储存分布函数,这就导致在计算的过程中会占用大量的虚拟内存;同时,边界条件的实现也需要依靠分布函数来转化,使得边界条件的处理变得更加复杂.鉴于传统TLBM所存在的这些不足,Chen等[19]最近提出一种简化的热格子玻尔兹曼方法(STLBM).在STLBM中,通过把控制方程的解在LBM的框架内进行重构,使得整个计算过程不再依赖分布函数,直接演化宏观量,从根本上改进了LBM,极大地减小了计算过程所占用的虚拟内存,简化了边界条件的实现方式.再者,通过Neumann稳定性分析,证明STLBM比传统LBM具有更高的稳定性[19].

基于对以上的研究分析,本文提出一种基于浸没边界法和简化热格子玻尔兹曼方法(IB-STLBM)的耦合模型对流固耦合传热机制进行探究.在IB-STLBM中,STLBM应用于流场和温度场的计算.STLBM在计算中只涉及平衡态分布函数和非平衡态分布函数(通过两个时间步上的平衡分布函数计算),整个计算过程直接对宏观量进行演化,故而只需要少量虚拟内存来储存宏观量,并且宏观边界条件也直接通过宏观量实现.同时,在处理浸没边界时,借鉴LBM的介观思想,把流场中的固体视为某种外部干扰项,干扰项的存在使得固体边界附近流体微团偏离了原来的平衡状态,而这一干扰项可以用一个非平衡分布函数来描述.依据这一思想,IBM的具体过程为:首先把流场和温度场的平衡态分布函数通过狄拉克函数插值到拉格朗日点上,再根据固体边界条件求出边界上的平衡态分布函数,接着通过反弹原理求出非平衡项,最后再把非平衡项分配到流场计算的欧拉网格上,修正固体边界.这样,在修正流固边界时也只是涉及平衡态分布函数和非平衡态分布函数,延续了IBM与STLBM的优越性,同时在介观层面上描述了流固界面的相互影响,具有清晰的物理含义.与传统的IB-LBM模型相比,本文所提出的IB-STLBM耦合模型具有更高的稳定性,适用于任意形状的边界,极大减小了计算过程中对虚拟内存的占用,使计算过程更加简洁直观.

1 简化热格子玻尔兹曼方法

1.1 热格子玻尔兹曼方法

基于BGK近似的一般不可压热格子玻尔兹曼模型的演化方程如下[13,20]

$$ f_\alpha \left( {x + e_\alpha \Delta t,t + \Delta t} \right) = f_\alpha \left( {x,t} \right) + \frac{f_\alpha ^{\rm eq} \left( {x,t} \right)-f_\alpha \left( {x,t} \right)}{\tau_{\rm f}} $$

$$g_\alpha (x + e_\alpha \Delta t,t + \Delta t) = g_\alpha (x,t) + \frac{g_\alpha ^{\rm eq} (x,t)-g_\alpha (x,t)}{\tau_{\rm g}} $$

其中,$f_\alpha$和$g_\alpha$分别代表在$\alpha$方向上的密度分布函数和温度分布函数,$f_\alpha ^{\rm eq}$和$g_\alpha ^{\rm eq}$代表其平衡态分布函数,$\Delta t$是迁移时间步长,$\tau _{\rm f}$和$\tau_ {\rm g}$分别是与运动粘度和热扩散系数相关的松弛系数,$e_\alpha$为格子离散速度,在D2Q9模型中,离散格子速度被定义为[21]

$ { e}_\alpha \left\{ {\begin{array}{llll} (0,0),& \alpha = 0 \\ (\pm 1,0),(0,\pm 1) ,&\alpha = 1,2,3,4 \\ (\pm 1,\pm 1),& \alpha = 5,6,7,8 \\ \end{array}} \right.$

平衡态分布函数定义为

$$ f_\alpha ^{\rm eq} = \omega _\alpha \rho \left[ {1 + \frac{e_\alpha \cdot u}{c_{\rm s}^2 } + \frac{(e_\alpha \cdot u)^2}{2c_{\rm s}^4 }-\frac{\left| u \right|^2}{2c_{\rm s}^2 }} \right]$$

$$ g_\alpha ^{\rm eq} = \omega _\alpha T\left[ {1 + \frac{e_\alpha \cdot u}{c_{\rm s}^2 } + \frac{(e_\alpha \cdot u)^2}{2c_{\rm s}^4 }-\frac{\left| u \right|^2}{2c_{\rm s}^2 }} \right]$$

上式中,$\rho$、$u$和$T$分别代表流体密度,速度和温度;$c_{\rm s}$为LBM中的声速[21];$\omega$$_{\alpha }$是权系数,在D2Q9 模型中$\omega _0 = 4 / 9$,$\omega _{1 \sim 4} = 1 / 9 , \omega _{5 \sim 8} = 1 / 36$;在LBM中,密度、速度和温度等宏观量根据分布函数求出,具体关系如下

$ \rho \mbox{ = }\sum\limits_{\alpha = 0}^N {f_\alpha ^{\rm eq} } , \quad T\mbox{ = }\sum\limits_{\alpha = 0}^N {g_\alpha ^{\rm eq} } , \quad \rho u = \sum\limits_{\alpha = 0}^N {e_\alpha f_\alpha ^{\rm eq} }$

在LBM中,Chapman-Enskog多尺度分析构建起连接LB方程和宏观控制方程的桥梁.通过Chapman-Enskog多尺度分析,LBM可以恢复到宏观控制方程[19]

$$\frac{\partial \rho }{\partial t} + \nabla \cdot \bigg(\sum\limits_\alpha {e_\alpha } f_\alpha ^{\rm eq} \bigg) = 0$$

$$\frac{{\partial \rho u}}{{\partial t}} + \nabla \cdot \sum\limits_\alpha {{{({e_\alpha })}_\chi }{{({e_\alpha })}_\gamma }\bigg[} f_\alpha ^{\rm eq}\bigg (1-\frac{1}{{2{\tau _f}}}\bigg)f_\alpha ^{\rm neq}\bigg] = {\bf 0}$$

$$\frac{\partial T}{\partial t} + \nabla \cdot \sum\limits_\alpha {e_\alpha \bigg[} g_\alpha ^{\rm eq} + \bigg(1-\frac{1}{2\tau _{\rm g}}\bigg)g_\alpha ^{\rm neq} \bigg] = 0 $$

其中,$(e_\alpha )_\chi$和$(e_\alpha )_\gamma$分别代表在$\chi$和$\gamma$ 方向上的格子速度,$f_\alpha ^{\rm neq}$和$g_\alpha ^{\rm neq}$分别是是流场和温度场的非平衡分布函数.与此同时,通过Chapman-Enskog分析,可以得出非平衡分布函数和平衡态分布函数之间的联系[19,22]

$$f_\alpha ^{\rm neq} = f_\alpha-f_\alpha ^{\rm eq} =-\tau _{\rm f} \Delta t\bigg(\frac{\partial }{\partial t} + e_\alpha \cdot \nabla\bigg )f_\alpha ^{\rm eq} $$

$$g_\alpha ^{\rm neq} = g_\alpha-g_\alpha ^{\rm eq} =-\tau _{\rm g}\Delta t\bigg(\frac{\partial }{\partial t} + e_\alpha \cdot \nabla \bigg)g_\alpha ^{\rm eq} $$

根据动量和能量守恒原则,非平衡分布函数满足以下方程

$ \sum\limits_\alpha {f_\alpha ^{\rm neq} } = { 0},\quad \sum\limits_\alpha {e_\alpha f_\alpha ^{\rm neq} } = {\bf0}, \quad \sum\limits_\alpha {g_\alpha ^{\rm neq} } = 0$

1.2 简化热格子玻尔兹曼方法(STLBM)

传统的热格子玻尔兹曼方法集成了标准LBM特性,但其缺陷也是明显的.首先,传统的热格子玻尔兹曼方法在计算的过程中需要消耗大量的虚拟内存.无论是在迁移还是碰撞步中,每个节点上所有方向的分布函数都需要分别储存起来.例如,在D2Q9模型中,需要在迁移和碰撞步中分别储存每个格子9个方向上的分布函数值,这就造成了LBM在模拟的过程中对虚拟内存有着较高的要求,特别是对于三维问题或是多组分多相流问题.其次,如何处理宏观边界也是传统热格子玻尔兹曼方法所面临的一个重要问题.在传统的热格子玻尔兹曼方法中,内部流场中每个格子节点上的分布函数可在计算的过程中直接得出,但是四周边界点上向内的分布函数是不能由演化过程直接得出,而下一步的计算必须在边界上的分布函数确定之后才能够顺利进行,所以需要先处理边界上未知的分布函数.同时,物理边界条件的实现需要依靠分布函数边界条件的转化.因此,边界分布函数的处理在LBM方法中是一个不可避免的问题,并会严重影响计算的稳定性和精度.

为了克服这些缺点,Chen等[19]提出了一种简化的热格子玻尔兹曼方法.此方法是通过重构基于LB方程恢复的宏观控制方程的解得到的,并且通过predictor-corrector两个步骤进行求解,整个计算过程不需要依赖于分布函数,直接演化宏观量.基于格子特性以及Chapman-Enskog分析所给出的关系,简化热格子玻尔兹曼方法的方程如下[19].

预测步骤

$$\rho ^ * \mbox{ = }\sum\limits_\alpha {f_\alpha ^{\rm eq} } (x-e_\alpha \Delta t,t-\Delta t) $$

$$\rho ^ * u^ * \mbox{ = }\sum\limits_\alpha {e_\alpha f_\alpha ^{\rm eq} } (x-e_\alpha \Delta t,t-\Delta t) $$

$$T^ * \mbox{ = }\sum\limits_\alpha {g_\alpha ^{\rm eq} } (x-e_\alpha \Delta t,t-\Delta t) $$

校正步骤

$$\rho \mbox{ = }\rho ^ * $$

$$\rho u\mbox{ = }\rho ^ * u^ * + \bigg (1-\frac{1}{\tau _{\rm f} }\bigg)\cdot\\ \sum\limits_\alpha {e_\alpha f_\alpha ^{\rm neq} } (x-e_\alpha \Delta t,t) + F_E \Delta t $$

$$T\mbox{ = }T^ * + \bigg(1-\frac{1}{\tau _{\rm g}}\bigg)\sum\limits_\alpha {g_\alpha ^{\rm neq} } (x-e_\alpha \Delta t,t) $$

其中非平衡分布函数的计算表达式为[19, 22-23]

$$f_\alpha ^{\rm neq} (x,t) =-\tau _{\rm f} \Delta t\bigg(\frac{\partial }{\partial t} + e_\alpha \cdot \nabla\bigg )f_\alpha ^{\rm eq} (x,t)= \\-\tau _{\rm f} [f_\alpha ^{\rm eq} (x,t)-f_\alpha ^{\rm eq} (x-e_\alpha \Delta t,t-\Delta t)] $$

$$g_\alpha ^{\rm neq} (x,t) =-\tau _{\rm g}\Delta t\bigg(\frac{\partial }{\partial t} + e_\alpha \cdot \nabla \bigg)g_\alpha ^{\rm eq} (x,t) =\\-\tau _{\rm g}[g_\alpha ^{\rm eq} (x,t)-g_\alpha ^{\rm eq} (x-e_\alpha \Delta t,t-\Delta t)] $$

其中当前时刻的平衡态分布函数是由预测步中的宏观量计算得到.基于Boussinesq近似,外力项$F_E$近似为与温度有关的浮力[24]

$ F_E \mbox{ = }\left( {\begin{array}{c} 0 \\-\rho g\beta (T-T_0 ) \\ \end{array}} \right)$

其中,$\beta$为热膨胀系数,g代表重力加速度,$T_0$是参考温度.

显然,从以上方程中可知,在STLBM计算过程中只涉及平衡态分布函数和非平衡态分布函数.平衡态分布函数是由宏观量直接计算,同时非平衡分布函数又是通过两个平衡分布函数在不同位置和时间水平上相减得到.基于此,STLBM的演化过程直接与宏观量相关联,使得STLBM具有许多优良特性[19].其一,在计算的过程中直接对宏观量进行演化,分布函数不复存在,所以只需要储存宏观量,极大减小了对虚拟内存的占用.其二,边界条件也可以直接实现,不需要经过分布函数的转化.其三,经过纽曼稳定性分析,STLBM具有无条件稳定性.基于STLBM的优良特性,在本文中采用STLBM来计算流场和热场中的宏观量.

1.3 浸没边界法

浸没边界的基本思想是利用某种体积力来模拟流场中边界对周围流体产生的影响,把边界作为某种力源项来处理,使复杂边界的处理变得简单.对于流场中的传热问题,基于IBM的控制方程如下

$$\frac{\partial \rho }{\partial t} + \nabla (\rho u) = 0 $$

$$\frac{\partial \rho u}{\partial t} + \nabla \cdot (\rho uu) =-\nabla p + \nabla \cdot [\rho \upsilon (\nabla u + (\nabla u)^{\rm T})] + f $$

$$\frac{\partial T}{\partial t} + \nabla \cdot (Tu) = \nabla \cdot (\kappa \nabla T) + q $$

其中,f表示浸没边界作用在流场上的力密度,$q$表示浸没边界上的热源作用在流场上的热密度,可通过以下公式求出

$$f(x,t) = \int\limits_\varGamma {F(X_l ,t)} \delta (x-X_l ){\rm d}s_l $$

$$q(x,t) = \int\limits_\varGamma {Q(X_l ,t)} \delta (x-X_l ){\rm d}s_l $$

式中,$x,X_l$分别代表欧拉点和拉格朗日点,$\delta (x-X_l )$ 是联系流场和浸没边界的狄克拉函数,$F(X_l ,t)$ 为浸没边界的力密度,$Q(X_l ,t)$ 是浸没边界的热通量.通过对方程(13)进行离散,可得到其离散格式

$$f(x,t) = \sum\limits_{i,j} {F(X_l ,t)D(x-X_l )\Delta s_l } $$

$$q(x,t) = \sum\limits_{i,j} {Q(X_l ,t)D(x-X_l )\Delta s_l } $$

其中,$D(x-X_l )$ 是一个逼近狄拉克函数的脉冲函数[25]

$$D(x-X_l )\mbox{ = }\frac{1}{h^2}\delta _h \left(\frac{x-X_l }{h}\right)\delta _h \left(\frac{y-X_l }{h}\right)$$

$$\delta (b) = \left\{ {\begin{array}{lll} 0.25[1 + \cos (0.5π b)],&\left| b \right| < 2 \\ 0,&\mbox{others}\\ \end{array}} \right.$$

式中$h = {\rm d}x$.

以上就是浸没边界法的主要计算流程,从中可知IBM的核心就是计算力密度和热通量.在求解力密度的研究中,Niu等[26]提出了一种简单易行的基于动量交换IB-LBM,把介观思想引入到IBM中.在此方法中,首先把流场中的分布函数通过狄拉克函数插值获得浸没边界上的分布函数,再利用"反弹"的方式获得新的满足无滑移边界条件的密度分布函数,然后基于动量交换原理求出浸没边界上力密度,最后利用式(14a)求出欧拉网格上的力密度.基于Niu等[26]的工作,Chen等[27-29]提出了一种基于分布函数修正的浸没边界格子玻尔兹曼方法.在Chen等[28]的方法中,用分布函数的非平衡项求解流固耦合作用,使整个演化的过程中与分布函数相结合,很好地保留了LBM的介观特性,同时又满足无滑移边界条件.基于以上工作,本文提出一种基于平衡态分布函数修正的IB-STLBM,具体细则将会在下一章中说明.

2 简化热格子玻尔兹曼方法和浸没边界法的耦合模型

2.1 基于平衡态分布函数修正的IB-STLBM的框架

在1.3部分,式(12a)~式(12c)给出基于浸没边界力密度的宏观控制方程.采用分步法(fractional step technique)[30-34],式(12a)~式(12c)可在predictor-corrector-IB corrector方案内进行重构.

预测步骤

$$\frac{\partial \rho }{\partial t} + \nabla \cdot \bigg(\sum\limits_\alpha {e_\alpha } f_\alpha ^{\rm eq} \bigg) = 0 $$

$$\frac{\partial \rho u}{\partial t} + \nabla \cdot \bigg[\sum\limits_\alpha {(e_\alpha )_\beta } (e_\alpha )_\beta f_\alpha ^{\rm eq} +\\ \frac{1}{2\tau _{\rm f}}\sum\limits_\alpha {(e_\alpha )_\beta } (e_\alpha )_\beta f_\alpha ^{\rm neq} ] = {\bf 0} $$

$$\frac{\partial T}{\partial t} + \nabla \cdot \bigg(\sum\limits_\alpha {e_\alpha } g_\alpha ^{\rm eq} + \frac{1}{2\tau _{\rm g}}\sum\limits_\alpha {e_\alpha } g_\alpha ^{\rm neq} \bigg) = 0 $$

校正步骤

$$\frac{\partial \rho }{\partial t} = 0 $$

$$\frac{\partial \rho u}{\partial t}\mbox{ = }f $$

$$\frac{\partial T}{\partial t} = q $$

对以上控制方程进行求解,并在LBM的框架内重构,便可得到基于平衡态分布函数修正的IB-STLBM的公式.

预测步骤

$$\rho ^ * \mbox{ = }\sum\limits_\alpha {f_\alpha ^{\rm eq} } (x-e_\alpha \Delta t,t-\Delta t) $$

$$\rho ^ * u^ * \mbox{ = }\sum\limits_\alpha {e_\alpha f_\alpha ^{\rm eq} } (x-e_\alpha \Delta t,t-\Delta t) $$

$$T^ * \mbox{ = }\sum\limits_\alpha {g_\alpha ^{\rm eq} } (x-e_\alpha \Delta t,t-\Delta t) $$

校正步骤

$$\rho ^{ * * }\mbox{ = }\rho ^ * $$

$$\rho ^{ * * }u^{ * * }\mbox{ = }\rho ^ * u^ *+ \bigg (1-\frac{1}{\tau _f }\bigg)\cdot\\ \sum\limits_\alpha {e_\alpha f_\alpha ^{\rm neq} } (x-e_\alpha \Delta t,t) + F_E \Delta t $$

$$T^{ * * }\mbox{ = }T^ * + \bigg(1-\frac{1}{\tau _{\rm g}}\bigg)\sum\limits_\alpha {g_\alpha ^{\rm neq} } (x-e_\alpha \Delta t,t) $$

IB校正步骤

$$\rho \mbox{ = }\rho ^{ * * } $$

$$\rho u\mbox{ = }\rho ^{ * *} u^ {* *} + \Delta t\sum\limits_\alpha {\delta f_{\alpha ,b} e_\alpha } $$

$$T\mbox{ = }T^ {* *} + \Delta t\sum\limits_\alpha {\delta g_{\alpha ,b} } $$

其中,非平衡函数$f_\alpha ^{\rm neq}$和$g_\alpha ^{\rm neq}$的计算与式(10a)和式(10b)相同.整个计算流程由3个步骤构成,先在预测步骤和校正步骤求出整个流场的宏观量,再通过浸没边界法求出固体边界对流场的影响($\delta f_{\alpha ,b}$和$\delta g_{\alpha ,b} )$,最后再在IB校正步骤中耦合到STLBM的计算当中,对固体边界进行修正.

2.2 基于平衡函数修正的浸没边界法

在介观理论体系中,流体被离散成一系列的微团(介观微粒),这些微团介于宏观和微观尺度之间,相对于宏观体系无穷小,但又比微观颗粒大.LBM正是通过演化计算这一系列介观微粒,模拟出与物理现象相符合的流体运动规律.根据此介观思想,当固体边界出现在流场中时,可以看作是对流场中固体边界附近的介观微粒产生了干扰,使其原来的平衡状态发生了变化.基于此,本文提出一种基于平衡函数修正的浸没边界法(equilibrium distribution function-correction IBM, EDC-IBM),把这种变化作为浸没边界对流场平衡态分布函数的干扰,从而产生一个非平衡项,再用非平衡项修正流场.

图1是浸没边界法所用到的网格,其中$\Omega _{1}$区域表示的是流场的欧拉网格,$\Omega _{2}$表示固体边界的拉格朗日点.在EDC-IBM中,首先通过狄拉克函数把$\Omega _{1}$边界周围的欧拉网格点的平衡态分布函数插值到拉格朗日点上(图1中黑点)

$$f_\alpha ^{\rm eq} (X_l ,t) = \sum\limits_{i,j} {f_\alpha ^{\rm eq} (x,t)D(x-X_l )h^2} $$

$$g_\alpha ^{\rm eq} (X_l ,t) = \sum\limits_{i,j} {g_\alpha ^{\rm eq} (x,t)D(x-X_l )h^2} $$

图1

图1   基于浸没边界的网格模型

Fig.1   Immersed boundary illustration


依据无滑移边界条件,把拉格朗日边界点上的速度、密度、温度等宏观量代入到平衡态分布方程3(a)和3(b)中,计算出其平衡态分布函数$F_\alpha ^{\rm eq} (X_l ,t)$,$G_\alpha ^{\rm eq} (X_l ,t)$;之后借鉴反弹原理,计算出所有拉格朗日点上的非平衡项

$$f_\alpha ^{\rm non-eq} ( X_l ,t) = f_{\overline \alpha }^{\rm eq} (X_l ,t)-F_{\overline \alpha }^{\rm eq} (X_l ,t) $$

$$g_\alpha ^{\rm non-eq} (X_l ,t) = g_{\overline \alpha }^{\rm eq} (X_l ,t)-G_{\overline \alpha }^{\rm eq} (X_l ,t) $$

其中,$\overline \alpha$是$\alpha$的反方向.对于Neumann边界条件(边界法线是给定的热通量$Q(X_l ,t)\mbox{ = }\frac{\partial T(X_l ,t)}{\partial n}$的边界条件,$\kappa$为热扩散系数).通过图1可以清晰看出整个计算域分成了浸没边界内部和外部两个区域,但是在浸没边界法中,并没有内外边界的区分,所以需要考虑边界内外对热源项所产生的影响,因而每个边界点上都存在两个法线方向,一个指向流场,一个指向浸没物体内部.边界上的两个热通量(即$Q(X_l ,t)$和$-\kappa \frac{\partial T^ * (X_l ,t)}{\partial n}$之间的差异所形成的热流,影响了周边欧拉网格温度场.因此,这类问题边界对温度场的非平衡项可以表示为[35]

$ g_\alpha ^{\rm non-eq} (X_l ,t)\mbox{ = }2\left[Q(X_l ,t) + \kappa \frac{\partial T^ * (X_l ,t)}{\partial n}\right]$

最后,通过狄拉克函数,边界对流场和温度场的非平衡项便可在欧拉节点上体现

$$f_{\alpha ,b} (x,t) = \sum\limits_\alpha {f_\alpha ^{\rm non-eq} (X_l ,t)D(x-X_l )\Delta s_l } $$

$$g_{\alpha ,b} (x,t) = \sum\limits_\alpha {g_\alpha ^{\rm non-eq} (X_l ,t)D(x-X_l )\Delta s_l } $$

把式(26a)和式(26b)代入IB校正步骤中,便可对边界点上的速度和温度等宏观量进行修正.

3.3 IB-STLBM算法流程

IB-STLBM的算法流程归纳如下:

(1) 选择迁移时间步长 ,给定各个初始参数值.

(2)根据式(20a)~式(20c)预测中间时刻的速度、密度和温度.

(3)根据式(10a)~式(10b)计算流场和温度场的非平衡态分布函数.

(4)根据式(21a)~式(21c)修正中间时刻的宏观量.

(5)根据式(23)~式(25)求出流固边界的非平衡干扰项.

(6)根据式(22a)~式(22c)和式(26a)和式(26b)修正流固耦合边界和宏观量.

(7)重复计算步骤(2)~(6)至计算结果达到收敛条件或者到最大计算时间.

4 数值方法的验证与分析

强制对流和自然对流问题是研究传热过程的两个重要方向,为了验证本文所提出的IB-STLBM的有效性与正确性,本章将针对Taylor-Green涡流、热圆柱绕流问题和含热源颗粒的封闭方腔内自然对流问题算例进行探讨.

4.1 Taylor-Green涡流

在本文中,采用STLBM计算流场的温度场,在Chen等[19]的研究中已经证实其具有二阶精度.同时本文的IBM方法也是汲取LBM的思想进行改进,而LBM本身具有二阶精度,在理论上来说本文发展的理论模型也具有二阶精度.但是由于计算的过程中存在许多不可预测的因素,理论上的推论并不可靠.为了进一步验证本模型的全局精度,在本节中将模拟Taylor-Green 涡流.计算域采用$[-L,L]\times [-L,L]$的正方形,其中$L = 1.0$,在计算域的中心有一个半径$R = 0.5$的圆柱.计算域的四周均采用周期性边界条件. 精确解的计算如下[36]

$$u =-\cos (xπ )\sin (yπ ){\rm e}^{-(2tπ ^2 / Re)} $$

$$v = \sin (xπ )\cos (yπ ){\rm e}^{-(2tπ ^2 / Re)} $$

$$\rho = \rho _0-\rho _0 \frac{U^2}{4c_{\rm s}^2 }[\cos (2xπ ) + \cos (2yπ )]{\rm e}^{-\frac{4tπ ^2}{Re}} $$

在本文的模拟中,单松弛时间$\tau \mbox{ = }0.65$,雷诺数$Re = 10$,采用四种不同的网格$(20\times 20, 40\times 40, 80\times 80, 160\times 160)$作为对比.由于在计算的过程中结果会随着时间演化而变化,所以在模拟中采用$t = 0$时刻的结果为初始值,$t = 1.0$为计算的结束时间. 误差通过下式求得

$ {\rm error} = \sqrt {\frac{\sum\limits_{N^2} {(u-u^{\rm exact})^2} }{N^2}}$

其中,$u$为模拟速度,$u^{\rm exact}$为精确解,$N$为网格长度.图2中展示的是数值误差在对数坐标下的分布情况,由图可以看到数值分布在斜率为1.86的直线附近,所以可以认为本模型具有二阶精度.

图2

图2   对数坐标下速度误差随网格变化分布图

Fig.2   L2 norms of relative error of u versus mesh steps in Taylor-Green vortex


4.2 热圆柱绕流问题

作为一个经典的强制对流算例,热圆柱绕流问题深受研究者的关注.对于热强制对流模型,速度场对温度场的影响是单向的,所以在研究中只考虑速度对温度影响.在研究热强制对流问题时,需要考虑两个重要的无量纲参数---雷诺数($Re)$和普朗特数($Pr)$.其中雷诺数用来表征流体流动情况,普朗特数体现了流体物理性质对对流传热过程的影响,其定义如下

$ Re = \frac{\rho dU}{\mu },Pr = \frac{\mu C_{\rm p} }{\kappa }$

其中,$\rho$表示流体密度,$U$为流体速度,d是圆柱的直径,$\mu$表示动力黏性系数,$C_{\rm p}$为等压比热容,$\kappa$表示物体的热导率.为了减小边界对流体流动的影响,本文选择$30d\times 20d$的计算区域,其中$d = 1$.圆柱的中心位置为($10d,10d)$.整个流场的初始速度$U$=0.01,初始温度$T_{0}$=0.整个计算区域的边界条件如表1所示. 表中$ u$为x方向的速度,v 为y方向的速度. 在此算例中热通量$Q_{\rm b} = 1$.

表1   等热通量圆柱绕流问题的边界条件

Table 1  The boundary condition of forced convection over a stationary heated circular cylinder

BoundaryBoundary condition
inflow boundaryu = U, v = 0, T = 0
outflow boundarydu/dx = 0,dv/dx = 0,dT/dx = 0
up boundarydu/dy = 0, v = 0,dT/dy = 0
down boundarydu/dy = 0, v = 0,dT/dy = 0
cylindricalT = T0,-KdT/dn = Qb

新窗口打开| 下载CSV


对于圆柱绕流问题,在雷诺数小于47的情况下,流场的流动以及所形成的涡街都是相对稳定的.本文的研究中,采用20和40的雷诺数作为对比,同时两种情况下普朗特数都取0.7.针对每一种情况,在模拟中都计算了阻力系数$C_{\rm d}$、圆柱尾流的回旋长度$L_{\rm w}$和圆柱表面的平均努塞尔数$\overline {Nu}$用于对比 ,$C_{\rm d}$和$\overline {Nu}$的计算公式为

$ C_{\rm d} = \frac{2F_{\rm d} }{\rho U^2d},\quad \overline {Nu} \mbox{ = }\frac{1}{2π d}\int {\frac{Q_{\rm b} }{\kappa (T-T_0 )}} {\rm d}s$

式中,$F_{\rm d}$表示圆柱的阻力,$s$为圆柱表面的弧长.表2中列出在不同网格下的$C_{\rm d}$和$\overline {Nu}$,通过表中数据可以反映出本模型具有较好的网格独立性,在不同的网格下相对误差较小.同时,在$900\times 600$与$1200\times 800$两种网格中相对误差在1%之内,能够得到较为精确的结果,为了减小计算量,在以下的算例中都采用$900\times 600$的网格数. 表3中分别列出阻力系数和回旋长度在不同$Re$时的模拟值.通过与文献[28,37-38]对比,显然本文的模拟结果与现有的研究结果十分吻合.在雷诺数相对较大的情况下,阻力系数和升力系数会随着涡街的脱落呈现周期性变化,图3中为$Re=200$时,阻力系数和升力系数随时间的变化图,可以看出二者呈现出很清晰的周期性.表4中展示本模型以及参考文献在雷诺数为20和40 下的圆柱表面$\overline {Nu}$的模拟值,本模型的结果与参考文献所提供的数值相吻合.

表2   Re= 20时阻力系数Cd和回流长度Nu在不同网格下的对比

Table 2  Comparisons of Cd and Nu at different mesh sizes

Mesh600 x 400900 x 6001200 x 800
Cd2.2432.1862.175
Nu2.6572.7142.727

新窗口打开| 下载CSV


表3   Re=40时阻力系数Cd和回流长度iw的对比

Table 3  Comparisons of Cd and Lw at Re=40

Reynold numberRefs.CdLw
present2.1861.021
Re=20Chen et al.[28]2.2250.992
Hu et al.[37]2.2131.016
Dennis et al. [38]2.0490.940
present1.6582.488
Re=40Chen et al.[28]1.6632.468
Hu et al. [37]1.6602.410
Dennis et al. [38]1.5222.345

新窗口打开| 下载CSV


表4   在不同雷诺数下Nu的对比

Table 4  Comparisons of Nu at different Reynolds number

Refs.Re=20Re=40
present2.7143.618
Chen et al.[28]2.7293.667
Ahmad et al.[39]2.6223.472
Ren et al.[35]2.7413.741

新窗口打开| 下载CSV


图3

图3   Re=200时升力系数和阻力系数随时间变化曲线

Fig.3   Evolution of drag and lift coefficients for flow at Re=200


图4中展示了在两种雷诺数下流场的流线.可以观察到在圆柱后会形成两个对称的涡,同时,在$Re =20$时,涡的回流长度大致与圆柱直径相等;在$Re =40$时,涡的回流长度变大,大致为直径的2.5倍.图5图6为不同雷诺数下热圆柱附近以及来流后方的等温线分布图.从中可以观察到,在圆柱前方等温线十分密集,而在后方等温线相对稀疏且间距逐渐增大.这主要是受到入口边界强制对流速度的影响,使得圆柱前方温度变化急剧,等温线密集,后方由于圆柱的阻挡,使得温度变化较缓,等温线稀疏.同时,雷诺数的变化也对温度传递造成影响,在雷诺数为20时,温度的传递明显小于雷诺数为40时,这一变化也与已有的研究一致.圆柱边界的$Nu$分布在图7中展示,通过与参考文献的对比,本文的结果与其十分接近,证实本模型的可行性与准确性.

图4

图4   圆柱绕流模型在$Re\mbox{ = }20$和$Re\mbox{ = }40$时的流线图

Fig.4   Streamline for the flow over a circular at$Re = 20$ and$Re = 40$


图5

图5   $Re\mbox{ = }20$时的等温线图

Fig.5   Temperature distribution at$Re = 20$


图6

图6   $Re\mbox{ = }40$时的等温线图

Fig.6   Temperature distribution at$Re = 40$


图7

图7   时$Nu$在上半圆柱的分布(顺时针方向)

Fig.7   Comparison of local$Nu$ distributes around the cylinder surface at$Re = 20$


4.3 含热源颗粒的封闭方腔内自然对流

自然对流现象是在没有外力驱动的情况下,由于物体体系内存在温差,冷热不均使得物体内部产生密度差,由浮升力引起物体内部流动,使得体系内产生传热现象.自然对流模型是流动和传热问题的耦合,瑞利数$Ra$对系统的流动传热有重要影响,当瑞利数较小时,能量传递主要以热传导为主,当瑞利数超过一定临界值时,主要受对流影响.瑞利数$Ra$的表达式如下

$ Ra = \frac{\beta g\Delta TL_0^3 Pr }{\nu }$

其中,$\Delta T$表示温度差,$L_{0}$是参考长度,$\nu$表示动力黏性系数.在自然对流模型中,其四周为固壁边界,方腔左边为热壁面($T_{\rm H}$=1.0),右边为冷壁面($T_{\rm L}$=0),流场的初始温度为左右壁面的平均温度,且在整个模拟过程中左右壁面温度保持不变.

为了探究流固干扰下方腔内传热性能,在上述模型内放置一个直径为0.05,密度为1.005的热颗粒.热颗粒的初始坐标为(0.5, 0.9),热通量边界条件为$-\kappa \partial T / \partial n\mbox{ = }5.0$. 在计算域为$1\times 1$的封闭方腔内,采用$200\times 200$的网格进行模拟,并对比不同瑞利数($Ra$=10$^{3}$,10$^{4}$,10$^{5})$的现象.

图8图9分别是无颗粒自然对流和加入热颗粒后($t = 180$ s)的自然对流模型在3个$Ra$下的等温线图. 对比图中的变化,可以归纳出$Ra$对整个体系传热的影响.在$Ra$=10$^{3}$时,能量传播主要以热传导为主,等温线几乎处于竖直状态,对流作用很弱;在$Ra$=10$^{4}$,10$^{5}$时,对流传热显著,方腔中部的等温线几乎是趋于水平,只有在四周壁面附近的等温线趋于垂直状态.在加入颗粒后,相同Ra情况下,流场等温线的变化趋势与图8基本相同.但是颗粒的加入影响了颗粒附近温度场的演化,同时,由于颗粒在方腔内运动,起到了"搅拌"的作用,加速了场内流体的流动,传热效果也显著提高.在$Ra$=10$^{3}$时,原来近乎垂直的等温线在颗粒的影响下有了大幅倾斜,随着$Ra$的增大($Ra$=10$^{4}$,10$^{5})$,颗粒对周边温度场的影响逐渐变小,这主要是因为在$Ra$较大的情况下,对流起着主要驱动作用.此外,在未加入颗粒前,流场经过一段时间后逐渐达到稳定状态,但在颗粒加入后,这种稳态被打破.

图8

图8   不同$Ra$的条件下自然对流等温线

Fig.8   Isotherms of natural convection in a square cavity at different Rayleigh number


图9

图9   不同$Ra$的条件下内含颗粒的自然对流等温线图

Fig.9   Isotherms of natural convection in a square cavity with a buoyant particle at different Rayleigh number


为了更好说明温度场的演化情况以及颗粒对系统传热的影响,图10中对比了自然对流和加入热颗粒后的自然对流的冷壁面平均$Nu$数在不同的$Ra$数下随时间的变化趋势.从图中可以看出,未加颗粒的自然对流,平均$Nu$值的变化趋势在一定时间后成一条水平直线;加颗粒后平均$Nu$数随着时间演化逐渐成周期性变化,这也反映出流场在后期并非是一个稳态,而是受颗粒的干扰呈现周期性演化.基于在加入颗粒后$Nu$呈周期性变化,且在50 s后变化周期趋于稳定,选取100~200 s的平均值做一个量化对比.表5中展示了与未加入颗粒自然对流模型的$Nu$值的对比.可以得出对于不同的瑞利数,$Nu$数都明显提高,这也表明在加入颗粒后由于颗粒搅拌作用,打破整个系统的稳态,导致传热效果有明显的提高.图10中还列举与Chen等[28]的结果对比,$Nu$数变化趋势与本模型基本一致,进一步验证了本模型是有效可行的.

表5   两种自然对流在不同此数下心值的对比

Table 5  Comparisons of the average Nusselt number of natural convection in a square cavity with and without buoyant particle at different Rayleigh number

RaNu(no particle)Nu (with particle)Increase
1031.1161.57741.29%
1042.2412.64818.15%
1054.5024.9269.42%

新窗口打开| 下载CSV


图10

图10   单纯自然对流与加入热颗粒后的自然对流的冷壁面平均$Nu$在不同的$Ra$下随时间的变化(图中虚线为单纯自然对流,实线为加颗粒后的自然对流,圆为Chen等[28]的结果)

Fig.10   Temporal variations of the average Nusselt number of natural convection in a square cavity with and without a buoyant particle at different Rayleigh number[28]


颗粒在3种$Ra$下的运动轨迹如图11所示.在图中可以清晰看出颗粒在方腔内做周期性运动,且与图9中$Nu$的变化周期相同.当$Ra = 10^3$时,颗粒的运动轨迹近似于一个类带圆角的正方形,随着时间推演,其周长有细微缩小;当$Ra = 10^4$时,颗粒似乎绕着某一点做近似圆周运动,其运动曲线呈现出一个旋涡转,每一个运动周期周长逐渐变小,这也和图10中$Nu$的变化趋势类似.当$Ra =10^5$时,前期的运动规律较为不规则,在经过一段时间后,几乎是沿着固定轨迹运动做周期运动,仅仅是有一些及其细微的变化.

图11

图11   在不同的$Ra$下颗粒运动轨迹

Fig.11   Temporal trajectory of the particle at different Rayleigh number


5 结 论

本文在简化热格子玻尔兹曼的基础上,针对流固耦合传热问题提出了一种IB-STLBM模型.基于分步技术的拓展应用,把控制方程的解在预测-校正-IB校正三个阶段内进行重构.在预测-校正步中STLBM被用于计算温度场和流场.STLBM能够对宏观量直接演化计算,消除了对分布函数的依赖,节省了大量虚拟内存,简化了边界条件的应用方式,同时相对于LBM有较高的稳定性.在IB 校正步中,采用浸没边界法对流固界面宏观量进行修正.在浸没边界法的应用中,借鉴介观的思想,把固体的介入看作是对附近流体微团流场平衡态的干扰,并把这种干扰用一个非平衡项来表示.通过对Taylor-Green涡流、热圆柱绕流问题和内含热源颗粒的自然对流模型的模拟,得到与现有文献研究相吻合的结果,很好的证明了这种方法的稳定性.同时,在内含热颗粒的封闭方腔自然对流的研究中,观察到颗粒的介入,打破了原始自然对流的稳态,颗粒在流场中的运动,使得自然对流最后形成一种周期性运动的状态,有效的提高其传热性能.当然,需要说明的是本模型虽然在各个方面表现出优异的特性,但是由于方法本身的限制,在涉及高雷诺数问题时仍然存在一些缺陷.对于这一问题,传统的做法有采用多松弛格式或者加入特殊高阶展开等方法,但是考虑到方法的简便与有效性,本模型并未使用这些方法,同时我们也在探寻一种简单的能够在高雷诺数下保持稳定的浸没边界与STLBM耦合模型,这是我们下一步研究的重心.但是,本模型仍不失为一种有效的模型可供研究者选用.

The authors have declared that no competing interests exist.
作者已声明无竞争性利益关系。

参考文献

邹勇, 朱桂平, 李来 .

液桥内热质耦合对流不稳定性及旋转磁场法控制

力学学报, 2017,49(6):1280-1289

DOI      URL     [本文引用: 1]

浮区法因具有无坩埚接触污染的生长优点而成为生长高完整性和高均匀性单晶材料的重要技术.但熔体中存在的毛细对流会给浮区法晶体生长带来极大挑战,这是由于对流的不稳定会导致晶体微观瑕疵的产生和宏观条纹等缺陷的形成.为了提高浮区法生长单晶材料的品质,研究浮区法晶体生长中毛细对流特性及如何控制其不稳定性显得尤为重要.本文采用数值模拟的方法对半浮区液桥内SixGe1-x体系中存在的热质毛细对流展开研究并施加旋转磁场对其进行控制.结果表明:纯溶质毛细对流表现为二维轴对称模式,温度场主要由热扩散作用决定,而浓度场则由对流和溶质扩散共同支配;纯热毛细对流呈现三维稳态非轴对称流动,浓度分布与熔体内热毛细对流的流向密切相关,等温线在对流较大的区域发生弯曲;耦合溶质与热毛细对流则为三维周期性旋转振荡流.施加旋转磁场后,熔体周向速度沿径向向外增大,熔体内浓度场和流场均呈现二维轴对称分布.

( Zou Yong, Zhu Guiping, Li Lai , et al.

Instability of coupled thermo-solute capillary convection in liquid bridge and control by rotating magnetic field

Chinese Journal of Theoretical and Applied Mechanics, 2017,49(6):1280-1289 (in Chinese))

DOI      URL     [本文引用: 1]

浮区法因具有无坩埚接触污染的生长优点而成为生长高完整性和高均匀性单晶材料的重要技术.但熔体中存在的毛细对流会给浮区法晶体生长带来极大挑战,这是由于对流的不稳定会导致晶体微观瑕疵的产生和宏观条纹等缺陷的形成.为了提高浮区法生长单晶材料的品质,研究浮区法晶体生长中毛细对流特性及如何控制其不稳定性显得尤为重要.本文采用数值模拟的方法对半浮区液桥内SixGe1-x体系中存在的热质毛细对流展开研究并施加旋转磁场对其进行控制.结果表明:纯溶质毛细对流表现为二维轴对称模式,温度场主要由热扩散作用决定,而浓度场则由对流和溶质扩散共同支配;纯热毛细对流呈现三维稳态非轴对称流动,浓度分布与熔体内热毛细对流的流向密切相关,等温线在对流较大的区域发生弯曲;耦合溶质与热毛细对流则为三维周期性旋转振荡流.施加旋转磁场后,熔体周向速度沿径向向外增大,熔体内浓度场和流场均呈现二维轴对称分布.

刘成, 叶正寅, 叶坤 .

转捩位置对全动舵面热气动弹性的影响

力学学报, 2017,49(4):802-810

DOI      URL    

高超声速附面层的转捩预测一直是流体力学研究中的难点,转捩前后物面的摩擦系数和传热系数会发生改变,转捩位置的不同会影响到飞行器表面热环境,进而使得飞行器的气动弹性特性发生显著变化.鉴于高超声速附面层转捩预测的不确定性,本文分析了转捩位置对高超声速全动舵面热气动弹性的影响.首先分别用层流模型和湍流模型求解N-S方程,得到气动热环境,并对气动热进行参数化;然后在不同转捩位置情况下构造出不同转捩位置的热分布模型,基于此种温度分布,结合热应力和材料属性的影响分析结构的热模态,将结构模态插值到气动网格上,采用基于CFD的当地流活塞理论进行气动弹性分析.以M=6,H=15km的某舵面为对象进行计算,结果表明:(1)随着转捩位置向后缘移动,结构频率上升,结构颤振速度呈增大趋势,转捩位置的变化能够带来颤振临界速度最大6%的变化量;(2)当转捩位置位于舵轴附近时,结构的颤振特性变化剧烈.通过刚度特性的分解和分析发现,导致颤振特性变化的主要因素在于舵轴的刚度特性变化,舵轴的影响量占整个结构刚度特性变化量的80%以上.

( Liu Cheng, Ye Zhengyin, Ye Kun .

The e ect of transiton location on aerothermoelasticity of a hypersonic all-movable centrol surface

Chinese Journal of Theoretical and Applied Mechanics, 2017,49(4):802-810 (in Chinese))

DOI      URL    

高超声速附面层的转捩预测一直是流体力学研究中的难点,转捩前后物面的摩擦系数和传热系数会发生改变,转捩位置的不同会影响到飞行器表面热环境,进而使得飞行器的气动弹性特性发生显著变化.鉴于高超声速附面层转捩预测的不确定性,本文分析了转捩位置对高超声速全动舵面热气动弹性的影响.首先分别用层流模型和湍流模型求解N-S方程,得到气动热环境,并对气动热进行参数化;然后在不同转捩位置情况下构造出不同转捩位置的热分布模型,基于此种温度分布,结合热应力和材料属性的影响分析结构的热模态,将结构模态插值到气动网格上,采用基于CFD的当地流活塞理论进行气动弹性分析.以M=6,H=15km的某舵面为对象进行计算,结果表明:(1)随着转捩位置向后缘移动,结构频率上升,结构颤振速度呈增大趋势,转捩位置的变化能够带来颤振临界速度最大6%的变化量;(2)当转捩位置位于舵轴附近时,结构的颤振特性变化剧烈.通过刚度特性的分解和分析发现,导致颤振特性变化的主要因素在于舵轴的刚度特性变化,舵轴的影响量占整个结构刚度特性变化量的80%以上.

徐飞彬, 周全, 卢志明 .

二维方腔热对流系统中纳米颗粒混合及凝并特性的数值模拟

力学学报, 2015,47(5):740-750

DOI      URL     Magsci    

<p>采用泰勒展开矩方法对二维瑞利-贝纳德热对流系统(1&times;10<sup>6</sup> &le;Ra &le;1 &times;10<sup>8</sup>) 中纳米颗粒群的混合和凝并特性进行了数值模拟. 结果显示颗粒群随时间演化经历了扩散阶段、混合阶段、充分混合阶段3 个阶段, 随着颗粒群混合和凝并的进行, 颗粒数目浓度减少, 颗粒群的平均体积增大; 得到了颗粒分布函数各特征量与温度相关系数以及各特征量的空间分布标准偏差在3 个阶段的不同特征; 得到了颗粒分布函数各阶矩以及平均体积长时间演化的渐近行为, 结果与零维渐近解析解一致. 最后, 本文进一步研究了无量纲数(包括瑞利数<em>Ra</em>, 斯密特数<em>Sc</em><sub>M</sub>, 达姆科勒数<em>Da</em>) 对颗粒群达到自保持分布时间的影响, 发现该时间随着<em>Ra</em>和<em>Sc</em><sub>M</sub>的增大呈对数率减小, 随着<em>Da</em>的增大呈线性增大</p>

( Xu Feibin, Zhou Quan, Lu Zhiming .

Numerical simulation of Brownian coagulation and mixing of Nanoparticles in 2-D Rayleigh-B'enard convection

Chinese Journal of Theoretical and Applied Mechanics, 2015,47(5):740-750(in Chinese))

DOI      URL     Magsci    

<p>采用泰勒展开矩方法对二维瑞利-贝纳德热对流系统(1&times;10<sup>6</sup> &le;Ra &le;1 &times;10<sup>8</sup>) 中纳米颗粒群的混合和凝并特性进行了数值模拟. 结果显示颗粒群随时间演化经历了扩散阶段、混合阶段、充分混合阶段3 个阶段, 随着颗粒群混合和凝并的进行, 颗粒数目浓度减少, 颗粒群的平均体积增大; 得到了颗粒分布函数各特征量与温度相关系数以及各特征量的空间分布标准偏差在3 个阶段的不同特征; 得到了颗粒分布函数各阶矩以及平均体积长时间演化的渐近行为, 结果与零维渐近解析解一致. 最后, 本文进一步研究了无量纲数(包括瑞利数<em>Ra</em>, 斯密特数<em>Sc</em><sub>M</sub>, 达姆科勒数<em>Da</em>) 对颗粒群达到自保持分布时间的影响, 发现该时间随着<em>Ra</em>和<em>Sc</em><sub>M</sub>的增大呈对数率减小, 随着<em>Da</em>的增大呈线性增大</p>

Mills ZG, Aziz B, Alexeev A .

Beating synthetic cilia enhance heat transport in microfluidic channels

Soft Matter, 2012,8(45):11508-11513

DOI      URL     [本文引用: 1]

\0

Peskin CS .

Flow pat terns around heart valves: a numerical method

Journal of Computational Physics, 1972,10(2):225-271

DOI      URL     [本文引用: 2]

The subject of this paper is the flow of a viscous incompressible fluid in a region containing immersed boundaries which move with the fluid and exert forces on the fluid. An example of such a boundary is the flexible leaflet of a human heart valve. It is the main achievement of the present paper that a method for solving the Navier-Stokes equations on a rectangular domain can now be applied to a problem involving this type of immersed boundary. This is accomplished by replacing the boundary by a field of force which is defined on the mesh points of the rectangular domain and which is calculated from the configuration of the boundary. In order to link the representations of the boundary and fluid, since boundary points and mesh points need not coincide, a semi-discrete analog of the 未 function is introduced. Because the boundary forces are of order h 1, and because they are sensitive to small changes in boundary configuration, they tend to produce numerical instability. This difficulty is overcome by an implicit method for calculating the boundary forces, a method which takes into account the displacements that will be produced by the boundary forces themselves. The numerical scheme is applied to the two-dimensional simulation of flow around the natural mitral valve.

Peskin CS, Printz ABF .

Improved volume conservation in the computation of flows with immersed elastic boundaries

Journal of Computational Physics, 1993,105(1):33-46

DOI      URL     [本文引用: 2]

This paper introduces a recipe for the construction of a finite-difference divergence operator the coefficients of which are completely determined once an interpolation scheme has been chosen. Substitution of this divergence operator for the previously used divergence based on central differences makes a dramatic improvement in the overall volume conservation that is observed in an immersed-boundary computation. This improvement is particularly important for computations in which an elastic boundary separates chambers containing fluid at substantially different pressures, a situation which is prominent in cardiac fluid dynamics during the contraction of the ventricles.

Chen S, Doolen GD .

Lattice boltzmann method for fluid flows

Annu. Rev. Fluid. Mech, 1998,30:329-364

DOI      URL     [本文引用: 1]

Abstract We present an overview of the lattice Boltzmann method (LBM), a parallel and efficient algorithm for simulating single-phase and multiphase fluid flows and for incorporating additional physical complexities. The LBM is especially useful for modeling complicated boundary conditions and multiphase interfaces. Recent extensions of this method are described, including simulations of fluid turbulence, suspension flows, and reaction diffusion systems.

Feng Z, Michaelides E .

The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems

Journal of Computational Physics, 2004,195(2):602-628

DOI      URL     [本文引用: 3]

A new computational method, the immersed boundary-lattice Boltzmann method, is presented. This method is a combination and utilizes the most desirable features of the lattice Boltzmann and the immersed boundary methods. The method uses a regular Eulerian grid for the flow domain and a Lagrangian grid to follow particles that are contained in the flow field. The rigid body conditions for the fluid and the particles are enforced by a penalty method, which assumes that the particle boundary is deformable with a high stiffness constant. The velocity field of the fluid and particles is solved by adding a force density term into the lattice Boltzmann equation. This novel method preserves the advantages of LBM in tracking a group of particles and, at the same time, provides an alternative and better approach to treating the solid-fluid boundary conditions. The method also solves the problems of fluctuation of the forces and velocities on the particles when the "bounce-back" boundary conditions are applied. This method enables one to simulate problems with particle deformation and fluid-structure deformation. Its results are validated by comparison with results from other methods

Feng Z, Michaelides E .

Proteus: A direct forcing method in the simulations of particulate flows

Journal of Computational Physics, 2005,202(1):20-51

DOI      URL     [本文引用: 2]

A new and efficient direct numerical method for the simulation of particulate flows is introduced. The method combines desired elements of the immersed boundary method, the direct forcing method and the lattice Boltzmann method. Adding a forcing term in the momentum equation enforces the no-slip condition on the boundary of a moving particle. By applying the direct forcing scheme, Proteus 1 In the Greek mythology Proteus is a hero, the son of Poseidon. In addition to his ability to change shapes and take different forms at will, Zeus granted him the power to make correct predictions for the future. One cannot expect better attributes from a numerical code. 1 eliminates the need for the determination of free parameters, such as the stiffness coefficient in the penalty scheme or the two relaxation parameters in the adaptive-forcing scheme. The method presents a significant improvement over the previously introduced immersed-boundary-lattice-Boltzmann method (IB-LBM) where the forcing term was computed using a penalty method and a user-defined parameter. The method allows the enforcement of the rigid body motion of a particle in a more efficient way. Compared to the ounce-back scheme used in the conventional LBM, the direct-forcing method provides a smoother computational boundary for particles and is capable of achieving results at higher Reynolds number flows. By using a set of Lagrangian points to track the boundary of a particle, Proteus eliminates any need for the determination of the boundary nodes that are prescribed by the ounce-back scheme at every time step. It also makes computations for particles of irregular shapes simpler and more efficient. Proteus has been developed in two- as well as three-dimensions. This new method has been validated by comparing its results with those from experimental measurements for a single sphere settling in an enclosure under gravity. As a demonstration of the efficiency and capabilities of the present method, the settling of a large number (1232) of spherical particles is simulated in a narrow box under two different boundary conditions. It is found that when the no-slip boundary condition is imposed at the front and rear sides of the box the particles motion is significantly hindered. Under the periodic boundary conditions, the particles move faster. The simulations show that the sedimentation characteristics in a box with periodic boundary conditions at the two sides are very close to those found in the sedimentation of two-dimensional circular particles.

Cheng YG, Zhu LD, Zhang CZ .

Numerical study of stability and accuracy of the immersed boundary method coupled to the Lattice Boltzmann BGK model

Communications in Computational Physics, 2014,16(1):136-168

DOI      URL     [本文引用: 1]

This paper aims to study the numerical features of a coupling scheme between the immersed boundary (IB) method and the lattice Boltzmann BGK (LBGK) model by four typical test problems: the relaxation of a circular membrane, the shearing flow induced by a moving fiber in the middle of a channel, the shearing flow near a non-slip rigid wall, and the circular Couette flow between two inversely rotating cylinders. The accuracy and robustness of the IB-LBGK coupling scheme, the performances of different discrete Dirac delta functions, the effect of iteration on the coupling scheme, the importance of the external forcing term treatment, the sensitivity of the coupling scheme to flow and boundary parameters, the velocity slip near non-slip rigid wall, and the origination of numerical instabilities are investigated in detail via the four test cases. It is found that the iteration in the coupling cycle can effectively improve stability, the introduction of a second-order forcing term in LBGK model is crucial, the discrete fiber segment length and the orientation of the fiber boundary obviously affect accuracy and stability, and the emergence of both temporal and spatial fluctuations of boundary parameters seems to be the indication of numerical instability. These elaborate results shed light on the nature of the coupling scheme and may benefit those who wish to use or improve the method.

Kang SK, Hassan YA .

A comparative study of direct-forcing immersed boundary-lattice Boltzmann methods for stationary complex boundaries

International Journal Numerical Methods Fluids, 2011,66(9):1132-1158

DOI      URL     [本文引用: 2]

In this study, we assess several interface schemes for stationary complex boundary flows under the direct-forcing immersed boundary-lattice Boltzmann methods (IB-LBM) based on a split-forcing lattice Boltzmann equation (LBE). Our strategy is to couple various interface schemes, which were adopted in the previous direct-forcing immersed boundary methods (IBM), with the split-forcing LBE, which enables us to directly use the direct-forcing concept in the lattice Boltzmann calculation algorithm with a second-order accuracy without involving the Navier tokes equation. In this study, we investigate not only common diffuse interface schemes but also a sharp interface scheme. For the diffuse interface scheme, we consider explicit and implicit interface schemes. In the calculation of velocity interpolation and force distribution, we use the 2- and 4-point discrete delta functions, which give the second-order approximation. For the sharp interface scheme, we deal with the exterior sharp interface scheme, where we impose the force density on exterior (solid) nodes nearest to the boundary. All tested schemes show a second-order overall accuracy when the simulation results of the Taylor reen decaying vortex are compared with the analytical solutions. It is also confirmed that for stationary complex boundary flows, the sharper the interface scheme, the more accurate the results are. In the simulation of flows past a circular cylinder, the results from each interface scheme are comparable to those from other corresponding numerical schemes. Copyright 2010 John Wiley & Sons, Ltd.

Shan XW .

Simulation of Rayleigh-Bénard convection using a lattice Boltzmann method

Physical Review E, 1997,55(3):2780

[本文引用: 1]

He XY, Chen SY, Doolen GD .

A novel thermal model for the lattice Boltzmann method in incompressible limit

Journal of Computational Physics, 1998,146(1):282-300

DOI      URL     [本文引用: 2]

A novel lattice Boltzmann thermal model is proposed for studying thermohydrodynamics in incompressible limit. The new model introduces an internal energy density distribution function to simulate the temperature field. The macroscopic density and velocity fields are still simulated using the density distribution function. Compared with the multispeed thermal lattice Boltzmann models, the current scheme is numerically more stable. In addition, the new model can incorporate viscous heat dissipation and compression work done by the pressure, in contrast to the passive-scalar-based thermal lattice Boltzmann models. Numerical simulations of Couette flow with a temperature gradient and Rayleigh B nard convection agree well with analytical solutions and benchmark data.

Jeong HK, Yoon HS, Ha MY , et al.

An immersed boundary-thermal lattice Boltzmann method using an equilibrium internal energy density approach for the simulation of flows with heat transfer

Journal of Computational Physics, 2010,229(7):2526-2543

DOI      URL     [本文引用: 1]

In present paper, a novel immersed boundary-thermal lattice Boltzmann method by the name of ldquoan equilibrium internal energy density approachrdquo is proposed to simulate the flows around bluff bodies with the heat transfer. The main idea is to combine the immersed boundary method (IBM) with the thermal lattice Boltzmann method (TLBM) based on the double population approach. The equilibrium internal energy density approach based on the equilibrium velocity approach [X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815] is used to combine IBM with TLBM. The idea of the equilibrium internal energy density approach is that the satisfaction of the energy balance between heat source on the immersed boundary point and the amount of change of the internal energy density according to time ensures the temperature boundary condition on the immersed boundary. The advantages of this approach are the simple concept, easy implementation and the utilization of original governing equation without modification. The simulation of natural convection in a square cavity with various body shapes for different Rayleigh numbers has been conducted to validate the capability and the accuracy of present method on solving heat transfer problems. Consequently, the present results are found to be in good agreement with those of previous studies. [All rights reserved Elsevier].

Kang SK, Hassan YA .

A direct-forcing immersed boundary method for the thermal lattice Boltzmann method

Computational Fluids, 2011,49(1):36-45

DOI     

In this study, a direct-forcing immersed boundary method (IBM) for thermal lattice Boltzmann method (TLBM) is proposed to simulate the non-isothermal flows. The direct-forcing IBM formulas for thermal equations are derived based on two TLBM models: a double-population model with a simplified thermal lattice Boltzmann equation (Model 1) and a hybrid model with an advection iffusion equation of temperature (Model 2). As an interface scheme, which is required due to a mismatch between boundary and computational grids in the IBM, the sharp interface scheme based on second-order bilinear and linear interpolations (instead of the diffuse interface scheme, which uses discrete delta functions) is adopted to obtain the more accurate results. The proposed methods are validated through convective heat transfer problems with not only stationary but also moving boundaries the natural convection in a square cavity with an eccentrically located cylinder and a cold particle sedimentation in an infinite channel. In terms of accuracy, the results from the IBM based on both models are comparable and show a good agreement with those from other numerical methods. In contrast, the IBM based on Model 2 is more numerically efficient than the IBM based on Model 1.

Seta T .

Implicit temperature-correction-based immersed-boundary thermal lattice Boltzmann method for the simulation of natural convection

Physical Review E, 2013,87(6):063304

DOI      URL     PMID     

In the present paper, we apply the implicit-correction method to the immersed-boundary thermal lattice Boltzmann method (IB-TLBM) for the natural convection between two concentric horizontal cylinders and in a square enclosure containing a circular cylinder. The Chapman-Enskog multiscale expansion proves the existence of an extra term in the temperature equation from the source term of the kinetic equation. In order to eliminate the extra term, we redefine the temperature and the source term in the lattice Boltzmann equation. When the relaxation time is less than unity, the new definition of the temperature and source term enhances the accuracy of the thermal lattice Boltzmann method. The implicit-correction method is required in order to calculate the thermal interaction between a fluid and a rigid solid using the redefined temperature. Simulation of the heat conduction between two concentric cylinders indicates that the error at each boundary point of the proposed IB-TLBM is reduced by the increment of the number of Lagrangian points constituting the boundaries. We derive the theoretical relation between a temperature slip at the boundary and the relaxation time and demonstrate that the IB-TLBM requires a small relaxation time in order to avoid temperature distortion around the immersed boundary. The streamline, isotherms, and average Nusselt number calculated by the proposed method agree well with those of previous numerical studies involving natural convection. The proposed IB-TLBM improves the accuracy of the boundary conditions for the temperature and velocity using an adequate discrete area for each of the Lagrangian nodes and reduces the penetration of the streamline on the surface of the body.

Zhang H, Yuan HZ, Trias FX , et al.

Particulate immersed boundary method for complex fluid-particle interaction problems with heat transfer

Computers & Mathematics with Applications, 2016,71:391-407

DOI      URL    

In our recent work (Zhang et02al., 2015), a Particulate Immersed Boundary Method (PIBM) for simulating fluid-particle multiphase flow was proposed and assessed in both two- and three-dimensional applications. In this study, the PIBM was extended to solve thermal interaction problems between spherical particles and fluid. The Lattice Boltzmann Method (LBM) was adopted to solve the fluid flow and temperature fields, the PIBM was responsible for the no-slip velocity and temperature boundary conditions at the particle surface, and the kinematics and trajectory of the solid particles were evaluated by the Discrete Element Method (DEM). Four case studies were implemented to demonstrate the capability of the current coupling scheme. Firstly, numerical simulation of natural convection in a two-dimensional square cavity with an isothermal concentric annulus was carried out for verification purpose. The current results were found to have good agreement with previous references. Then, sedimentation of two-and three-dimensional isothermal particles in fluid was numerically studied, respectively. The instantaneous temperature distribution in the cavity was captured. The effect of the thermal buoyancy on particle behaviors was discussed. Finally, sedimentation of three-dimensional thermosensitive particles in fluid was numerically investigated. Our results revealed that the LBM–PIBM–DEM is a promising scheme for the solution of complex fluid–particle interaction problems with heat transfer.

Hu Y, Li DC, Shu S , et al.

An efficient immersed boundary-lattice boltzmann method for the simulation of thermal flow problems

Communications in Computational Physics, 2016,20:1210-1257

DOI      URL     [本文引用: 1]

In this paper, a diffuse-interface immersed boundary method (IBM) is proposed to treat three different thermal boundary conditions (Dirichlet, Neumann, Robin) in thermal flow problems. The novel IBM is implemented combining with the lattice Boltzmann method (LBM). The present algorithm enforces the three types of thermal boundary conditions at the boundary points. Concretely speaking, the IBM for the Dirichlet boundary condition is implemented using an iterative method, and its main feature is to accurately satisfy the given temperature on the boundary. The Neumann and Robin boundary conditions are implemented in IBM by distributing the jump of the heat flux on the boundary to surrounding Eulerian points, and the jump is obtained by applying the jump interface conditions in the normal and tangential directions. A simple analysis of the computational accuracy of IBM is developed. The analysis indicates that the Taylor-Green vortices problem which was used in many previous studies is not an appropriate accuracy test example. The capacity of the present thermal immersed boundary method is validated using four numerical experiments: (1) Natural convection in a cavity with a circular cylinder in the center; (2) Flows over a heated cylinder; (3) Natural convection in a concentric horizontal cylindrical annulus; (4) Sedimentation of a single isothermal cold particle in a vertical channel. The numerical results show good agreements with the data in the previous literatures.

Chen Z, Shu C, Tan D .

A simplified thermal lattice Boltzmann method without evolution of distribution functions

International Journal of Heat Mass Transfer, 2017,105:741

DOI      URL     [本文引用: 9]

In this paper, a simplified thermal lattice Boltzmann method (STLBM) without evolution of the distribution functions is developed for simulating incompressible thermal flows. With the assistance of the fractional step technique, the macroscopic governing equations recovered from Chapman–Enskog (C–E) expansion analysis are resolved through a predictor–corrector scheme. Then in both the predictor and corrector steps, using the isentropic properties of lattice tensors and relationships of C–E analysis, the macroscopic flow variables are explicitly calculated from the equilibrium and non-equilibrium distribution functions. In STLBM, the equilibrium distribution functions are calculated from the macroscopic variables, while the non-equilibrium distribution functions are evaluated from the differences between two equilibrium distribution functions at different locations and time levels. Therefore, STLBM directly updates the macroscopic variables during the computational process, which lowers the virtual memory cost and facilitates the implementation of physical boundary conditions. Through von Neumann stability analysis, the present method is proven to be unconditionally stable, which is further validated by numerical tests. Three representative examples are presented to demonstrate the robustness of STLBM in practical simulations and its flexibility on different types of meshes and boundaries.

Peng Y, Shu C, Chew Y .

Simplified thermal lattice Boltzmann model for incompressible thermal flows

Physical Review E, 2003,68:026701

DOI      URL     PMID      [本文引用: 1]

Abstract Considering the fact that the compression work done by the pressure and the viscous heat dissipation can be neglected for the incompressible flow, and its relationship with the gradient term in the evolution equation for the temperature in the thermal energy distribution model, a simplified thermal energy distribution model is proposed. This thermal model does not have any gradient term and is much easier to be implemented. This model is validated by the numerical simulation of the natural convection in a square cavity at a wide range of Rayleigh numbers. Numerical experiments showed that the simplified thermal model can keep the same order of accuracy as the thermal energy distribution model, but it requires much less computational effort.

Qian Y, Humieres D. Lallemand P .

Lattice BGK models for Navier-Stokes equation

Europhysicas Letters, 1992,17(6):479-484

[本文引用: 2]

Guo ZL, Shu C .

Lattice Boltzmann Method and Its Applications in Engineering

World Scientific, 2013

DOI      URL     [本文引用: 2]

Lattice Boltzmann method and its applications in engineering Zhaoli Guo, Chang Shu (Advances in computational fluid dynamics / editors-in -chief, Chi-Wang Shu and Chang Shu, v. 3) World Scientific, c2013

Wang Y, Shu C, Teo C .

Thermal lattice Boltzmann flux solver and its application for simulation of incompressible thermal flows

Computer & Fluids, 2014,94:98-111

DOI      URL     [本文引用: 1]

A thermal lattice Boltzmann flux solver (TLBFS) is developed in this work for simulation of incompressible thermal flows. In TLBFS, the thermal lattice Boltzmann method (TLBM) is only applied to reconstruct the local solution of TLBM for evaluation of fluxes at the cell interface. Meanwhile, the macroscopic flow variables at cell centers are obtained by using the finite volume method to solve conservative differential equations recovered from Chapman–Enskog analysis of thermal lattice Boltzmann equation. The physical boundary conditions in TLBFS can be directly implemented using the same way as in conventional Navier–Stokes (N–S) solvers. The present solver eliminates the constraints associated with conventional lattice Boltzmann method such as limitation to uniform Cartesian mesh, tie-up between the time step and the mesh spacing, as well as implementation of boundary conditions for distribution functions. TLBFS is validated through numerical examples of natural convection in enclosures, including a square cavity and a cylindrical annulus, and mixed convection from a heated cylinder. Through numerical validation, it is shown that TLBFS can be effectively and flexibly applied to solve thermal flow problems with curved boundaries.

Gray DD, Giorgini A .

The validity of the Boussinesq approximation for liquids and gases

International Journal of Heat Mass Transfer, 1976,19:545

DOI      URL     [本文引用: 1]

A new method for obtaining approximate equations for natural convection flows is presented. The systematic application of this method leads to explicit conditions for the neglect of various terms. It is shown that this method allows the specification of the conditions under which the traditional Boussinesq approximation applies to a given Newtonian liquid or gas. The method is applied to room temperature water and air.

Peskin CS .

The immersed boundary method

Acta Numer, 2002,11:479-517

[本文引用: 1]

Niu X, Shu C, Chew Y , et al.

A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows

Physics Letters A, 2006,354(3):173-182

DOI      URL     [本文引用: 2]

A momentum exchange-based immersed boundary-lattice Boltzmann method is presented in this Letter for simulating incompressible viscous flows. This method combines the good features of the lattice Boltzmann method (LBM) and the immersed boundary method (IBM) by using two unrelated computational meshes, an Eulerian mesh for the flow domain and a Lagrangian mesh for the solid boundaries in the flow. In this method, the non-slip boundary condition is enforced by introducing a forcing term into the lattice Boltzmann equation (LBE). Unlike the conventional IBM using the penalty method with a user-defined parameter or the direct forcing scheme based on the Navier tokes (NS) equations, the forcing term is simply calculated by the momentum exchange of the boundary particle density distribution functions, which are interpolated by the Lagrangian polynomials from the underlying Eulerian mesh. Numerical examples show that the present method can provide very accurate numerical results.

Chen MF, Niu XD .

An improved momentum-exchanged immersed boundary-based lattice Boltzmann method for incompressible viscous thermal flows

Int. J. Modern Phys, 2016,42:1660161

DOI      URL     [本文引用: 1]

An improved momentum-exchanged immersed boundary-based lattice Boltzmann method (MEIB-LBM) for incompressible viscous thermal flows is presented here. MEIB-LBM was first proposed by Niu et al, which has been shown later that the non-slip boundary condition is not satisfied. Wang. et al. and Hu. et al overcome this drawback by iterative method. But it needs to give an appropriate relaxation parameter. In this work, we come back to the intrinsic feature of LBM, which uses the density distribution function as a dependent variable to evolve the flow field, and uses the density distribution function correction at the neighboring Euler mesh points to satisfy the non-slip boundary condition on the immersed boundary. The same idea can also be applied to the thermal flows with fluid-structure interference. The merits of present improvements for the original MEIB-LBM are that the intrinsic feature of LBM is kept and the flow penetration across the immersed boundaries is avoided. To validate the present method, examples, including forced convection over a stationary heated circular cylinder for heat flux condition, and natural convection with a suspended circle particle in viscous fluid, are simulated. The streamlines, isothermal contours, the drag coefficients and Nusselt numbers are calculated and compared to the benchmark results to demonstrate the effective of the present method.

Chen MF, Niu XD, Yamaguchi H , et al.

A lattice Boltzmann modeling fluid-structure interaction problem and its applications in natural convections in a square cavity with particles suspended inside

Adv App Math Mech, 2018,10(2):303-328

[本文引用: 8]

陈木凤, 李翔, 牛小东 .

两个非磁性颗粒在磁流体中的沉降现象研究

物理学报, 2017,66(16):164703

URL     Magsci     [本文引用: 1]

在磁场作用下,在磁流体里添加非磁性颗粒(non-magnetic particles,NPs),可以使得NPs形成不同的结构,操控NPs的运动从而影响磁流体的特性,这种应用逐渐受到了研究者的关注.为了更好地操控磁流体里NPs的运动,本文采用一种多物理模型研究在外加磁场作用下,磁流体中两个NPs沉降的运动过程.其中,用格子玻尔兹曼方法模拟磁流体的运动,外加磁场对磁流体的影响用一种自修正方法求解泊松方程,这个自修正方法可以使欧姆定律满足守恒定律.NPs之间的偶极干扰力采用偶极力模型,同时采用一种相对过渡平滑的共轭边界条件处理NPs与磁流体交界面的流固干扰以避免磁场密度过渡的突变.本文主要探究两个NPs在磁流体中的沉降,揭示磁场作用下NPs的相互干扰原理;同时,对控制NPs运动时的参数进行调节,得到NPs不同的运动轨迹,达到操控颗粒运动的目的.本研究可对NPs在磁流体中的应用提供定量的分析结果,对NPs在工业上的应用提供有力的理论支撑.

( Chen Mufeng, Li Xiang, Niu Xiaodong , et al.

Sedimentation of two non-magnetic particles in magnetic fluid

Acta Phy. Sin, 2018,66(16):164703 (in Chinese))

URL     Magsci     [本文引用: 1]

在磁场作用下,在磁流体里添加非磁性颗粒(non-magnetic particles,NPs),可以使得NPs形成不同的结构,操控NPs的运动从而影响磁流体的特性,这种应用逐渐受到了研究者的关注.为了更好地操控磁流体里NPs的运动,本文采用一种多物理模型研究在外加磁场作用下,磁流体中两个NPs沉降的运动过程.其中,用格子玻尔兹曼方法模拟磁流体的运动,外加磁场对磁流体的影响用一种自修正方法求解泊松方程,这个自修正方法可以使欧姆定律满足守恒定律.NPs之间的偶极干扰力采用偶极力模型,同时采用一种相对过渡平滑的共轭边界条件处理NPs与磁流体交界面的流固干扰以避免磁场密度过渡的突变.本文主要探究两个NPs在磁流体中的沉降,揭示磁场作用下NPs的相互干扰原理;同时,对控制NPs运动时的参数进行调节,得到NPs不同的运动轨迹,达到操控颗粒运动的目的.本研究可对NPs在磁流体中的应用提供定量的分析结果,对NPs在工业上的应用提供有力的理论支撑.

Chorin AJ .

Numerical solution of the Navier-Stokes equations

Math. Compo, 1968,22:745-762

[本文引用: 1]

Chorin AJ .

Numerical solution of incompressible flow problems

Stud. Numer. Anal, 1968,2:64-71

URL    

Temam R .

Sur l'approximation de la solution des equations de Navier-Stokes par la methode des pas fractionnaires (I)

Arch. Ration. Mech. Anal, 1969,32:135-153

Temam R .

Sur. l'approximation de la solution des equations de Navier-Stokes par la methode des pas fractionnaires (II)

Arch. Ration. Mech. Anal, 1969,33:377-385

Kim J, Moin P .

Application of a fractional-step method to incompressible Navier-Stokes equations

Journal of Computational Physics, 1985,59:308-323

DOI      URL     [本文引用: 1]

A numerical method for computing three-dimensional, time-dependent incompressible flows is presented. The method is based on a fractional-step, or time-splitting, scheme in conjunction with the approximate-factorization technique. It is shown that the use of velocity boundary conditions for the intermediate velocity field can lead to inconsistent numerical solutions. Appropriate boundary conditions for the intermediate velocity field are derived and tested. Numerical solutions for flows inside a driven cavity and over a backward-facing step are presented and compared with experimental data and other numerical results.

Ren W, Shu C, Yang W .

An efficient immersed boundary method for thermal flow problems with heat flux boundary conditions

Heat & Mass Transfer, 2013,64:694-705

DOI      URL     [本文引用: 2]

A boundary condition implemented-immersed boundary method (IBM) involving velocity correction and heat flux correction is presented in this paper. In the framework of IBM, the velocity correction is made with Dirichlet condition (non-slip), and the temperature correction is made with Neumann (heat flux) condition. The main feature of present approach is to accurately satisfy the governing equations and boundary conditions through velocity and heat flux correction, which is performed by introducing a forcing term in the momentum equation and a heat source/sink term in the energy equation to consider the effect of the immersed boundary. The forcing term and heat source/sink are determined in such a way that the physical boundary conditions for velocity and temperature can be accurately satisfied. Numerical experiments for both forced convection and natural convection problems are conducted to validate the capability and efficiency of present method. Good agreements with available data in the literature have been achieved.

Chen DJ, Lin KH, Lin CA .

Immersed boundary method based lattice Boltzmann to simulate 2D and 3D complex geometry flows

International Journal of Modern Physics C, 2007,18:585-594

DOI      URL     [本文引用: 1]

In this paper, the lattice Boltzmann method is combined with the immersed boundary technique to simulate complex geometry flows. The complex geometry is represented by Lagrangian markers and forces are exerted at the Lagrangian markers in order to satisfy the prescribed velocity of the boundary. This force at the Lagrangian markers is then distributed to the Eulerian grid by a well-chosen discretized delta function. With the known force field in the Eulerian grid to mimic the boundary, the lattice Boltzmann method is used to compute the flow field where the complex geometry is immersed inside the Cartesian computational domain. Numerical experiments show that the second-order accuracy of the adopted numerical scheme is degraded to 1.8 order. The proposed method is examined by computing decaying vortex, lid driven cavity flow and 2D and 3D flows over asymmetrically placed cylinder. All the numerical results are compatible with the benchmark solutions.

Hu Y, Yuan H, Shu S , et al.

An improved momentum exchanged-based immersed boundary-lattice Boltzmann method by using an iterative technique

Computers & Mathematics with Applications, 2014,68(3):140-155

DOI      URL     [本文引用: 3]

A novel immersed boundary–lattice Boltzmann method (IB–LBM) is proposed for incompressible viscous flows in complex geometries. Based on the momentum exchanged-based IB–LBM, an iterative technique is introduced to enforce the non-slip boundary condition at the boundary points. Moreover, the proposed IB–LBM overcomes the drawback that the numerical results of the previous work (Wu and Shu, 2009) which is affected by the distribution of Lagrangian points. A simple theoretical analysis is developed to obtain the optimal iteration parameters. Numerical results show that the present scheme has second-order accuracy and is not affected by the distribution of Lagrangian points. It also shows that the non-slip boundary condition is satisfied on the boundary. This verifies that our IB–LBM is capable of simulating flow problems with complex boundaries.

Dennis SCR, Chang GZ .

Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100

Journal of Fluid Mechanics, 1970,42(3):471-489

DOI      URL     [本文引用: 3]

Finite-difference solutions of the equations of motion for steady incompressible flow around a circular cylinder have been obtained for a range of Reynolds numbers from R = 5 to R = 100. The object is to extend the Reynolds number range for reliable data on the steady flow, particularly with regard to the growth of the wake. The wake length is found to increase approximately linearly with R over the whole range from the value, just below R = 7, at which it first appears. Calculated values of the drag coefficient, the angle of separation, and the pressure and vorticity distributions over the cylinder surface are presented. The development of these properties with Reynolds number is consistent, but it does not seem possible to predict with any certainty their tendency as R [rightward arrow] [infty infinity]. The first attempt to obtain the present results was made by integrating the time-dependent equations, but the approach to steady flow was so slow at higher Reynolds numbers that the method was abandoned.

Ahmad RA, Qureshi ZH .

Laminar mixed convection from a uniform heat flux horizontal cylinder in a crossflow

Journal of Thermophysics & Heat Transfer, 1992,6(2):277-287

DOI      URL     [本文引用: 1]

An analysis of combined forced and natural convection (mixed convection) heat transfer from a horizontal cylinder dissipating a uniform heat flux is conducted by solving the full two-dimensional steady-state Navier-Stokes and energy equations. For Pr = 0.7 flow fields and transport results are determined in the ranges of Reynolds numbers 1-60 and the modified Grashof numbers 0-16,000. The results are presented as local and average values of the Nusselt number, the vorticity, the pressure distribution and the coefficient of drag around the cylinder.

Bharti R, Chhabra RP, Eswaran V .

A numerical study of the steady forced convection heat transfer from an unconfined circular cylinder

Heat & Mass Transfer, 2007,43:639-648

DOI      URL    

Forced convection heat transfer from an unconfined circular cylinder in the steady cross-flow regime has been studied using a finite volume method (FVM) implemented on a Cartesian grid system in the range as 1002≤02 Re 02≤0245 and 0.702≤02 Pr 02≤02400. The numerical results are used to develop simple correlations for Nusselt number as a function of the pertinent dimensionless variables. In addition to average Nusselt number, the effects of Re, Pr and thermal boundary conditions on the temperature field near the cylinder and on the local Nusselt number distributions have also been presented to provide further physical insights into the nature of the flow. The rate of heat transfer increases with an increase in the Reynolds and/or Prandtl numbers. The uniform heat flux condition always shows higher value of heat transfer coefficient than the constant wall temperature at the surface of the cylinder for the same Reynolds and Prandtl numbers. The maximum difference between the two values is around 15–20%.

/