中国科学院流固耦合系统力学重点实验室

基于 CFD 降阶的非线性气动弹性稳定性分析周强;陈刚;李跃明【摘要】为了快速寻找基于 CFD／CSD 的非线性气弹系统的颤振边界，根据Lyapunov 稳定性理论对非线性流固耦合系统进行稳定性分析。

首先通过微扰理论对非线性流固耦合系统处理建立近似线性化状态空间方程，再通过 POD （Proper Orthogonal Decomposition）方法将高维状态空间方程降阶为便于分析的降阶系统，根据系统所有特征值即可判定原始非线性系统稳定性。

Lyapunov 稳定性理论主要针对非线性系统，在实现过程中采用了 POD 降阶的方法，与大多数对降阶系统稳定性判定不同，其数学理论上是反映原始非线性流固耦合系统稳定性。

POD 降阶方法从系统内部流场出发，可以较好反映系统内部特性。

二维三维算例仿真结果验证了该方法预测颤振边界的正确性，分析发现，在亚音速阶段，稳定性主要由于结构模态主导；而在跨音速和超音速流动阶段，气弹稳定性主要由受流体特性影响。

%In order to quickly find the flutter boundary of a nonlinear aeroelstic system based on CFD /CSD,this paper used the Lyapunov stability theory to analyze the stability of a nonlinear fluid-structure coupling system.Firstly,a linearized state space model was obtained based on a nonlinear fluid-structure coupling system through the small perturbation theory;then a reduced order model was achieved by reducing the high dimensional linearized model through a POD(proper orthogonal decomposition)method.According to the eigenvalues of the ROM system,the stability of the original nonlinear system could be determined.Lyapunov stability theory was mainly aimed at the nonlinear system,and the POD method was adopted to realize the process.Differentfrom other stability methods,the mathematical theory reflected the stability of the original nonlinear fluid-structure coupling system.As the POD /ROM used in this paper was based on the CFD fluid field,so it could better reflect the internal characteristics of the original nonlinear coupling system.Numerical cases including two-dimensional aerofoil and three-dimensional wing's models were used to validate the accuracy of the stability method.The results show that in the subsonic domain,the stability is determined by the structural model;but in transonic and supersonic domain,the aeroelastic stability is mainly affected by fluid characteristics.【期刊名称】《振动与冲击》【年(卷),期】2016(035)016【总页数】7页(P17-23)【关键词】流固耦合;李雅普洛夫;正则正交分解;状态空间;颤振【作者】周强;陈刚;李跃明【作者单位】西安交通大学机械结构强度与振动国家重点实验室，西安 710049;西安交通大学机械结构强度与振动国家重点实验室，西安 710049;西安交通大学机械结构强度与振动国家重点实验室，西安 710049【正文语种】中文【中图分类】V211;V215.3飞行器气动弹性是研究飞行器结构与气动力相互耦合关系，是典型的流固耦合振动问题。

基于相互作用积分方法的裂纹扩展分析陈旻炜;李敏;陈伟民【摘要】In the process of analyzing crack growth with traditional finite element method,crack tip singularity is an urgent problem that must be faced with,and then the direction and length of the propaga-tion are closely related to the mesh of model.In this paper,a model that is not thoroughly dependent on the mesh and can be used for analyzing problems of crack growth is set up.Because interaction integral method can relieve the effects of singularity to some extent,severe requirements for meshing caused by singularity can be avoided.Stress intensity factor at the crack tip is solved by interaction integral via finite element method,and direction of the propagation is determined by the maximum circumferential tensile stress theory.When cracks propagate,only the meshes of concerned areas are regenerated in order to ob-tain the final extension path in a relatively precise way.At last,several typical examples,including growth of mode-I crack,cracks with a hole plate and cracks subj ected to compression,are compared with refer-ences.Reliability of the method in this paper is verified and applicability in many damaged patterns is shown as well.%使用传统有限元方法对裂纹扩展问题进行分析时,裂尖奇异性是首先必须要面对的问题,其次裂纹扩展的方向、长度与模型的网格密切相关.本文通过建立一个不完全依赖于网格的模型来分析裂纹扩展问题,利用相互作用积分方法来求解裂尖参数,积分方法可以一定程度上消除奇异性影响,从而避免了裂纹尖端的奇异性对网格划分的苛刻要求.采用相互作用积分法则求解裂纹尖端的应力强度因子,由最大周向拉应力理论来判断裂纹的扩展方向.模拟裂纹扩展的过程中,只对涉及扩展区域的网格进行重新剖分来获得最终的扩展路径以保证计算的精度.最后通过I型裂纹扩展、带孔板裂纹扩展及受压裂纹扩展等几个典型算例与参考文献进行对比,证实了本文方法具有良好的可靠性,并展示了在多种破坏模式下方法的适用性.【期刊名称】《湖南大学学报（自然科学版）》【年(卷),期】2018(045)004【总页数】8页(P74-81)【关键词】裂尖奇异性;相互作用积分;应力强度因子;有限元方法;裂纹扩展【作者】陈旻炜;李敏;陈伟民【作者单位】北京航空航天大学航空科学与工程学院,北京 100191;北京航空航天大学航空科学与工程学院,北京 100191;中国科学院力学研究所流固耦合力学重点实验室,北京 100190【正文语种】中文【中图分类】O346.1关于裂纹扩展问题的研究，对于预防工程结构的断裂失效、评价结构的完整性与安全性具有非常重要的意义.在非常规油气资源的开采中，往往通过水力压裂技术增大岩层中裂纹的表面积来实现提高产量的目的，这也促进了众学者对裂纹扩展问题的研究.裂纹扩展问题属于不定边界问题，具有高度的非线性[1].它是一个连续面与非连续面不断转换的过程，而该问题的关键是如何实现对不连续面的准确刻画.通常用来分析裂纹扩展的方法主要有边界元方法和有限单元方法.边界元法(BEM)[2-3]在模拟复杂的模型时，严重依赖于问题的基本解，大大降低了处理非线性、非均质等问题时的优势.有限元方法中的单元消除(杀死)技术也可以简单、有效地实现有限单元法对于断裂破坏过程的模拟[4].为了解决有限元模型中重划分的问题[5]，学者们在传统有限元基础上引入了结点释放技术[6]以及扩展有限元法[7-8].结点释放技术需要预先知道裂纹潜在的扩展方向，因此在模拟复杂情况下的裂纹扩展问题时会比较麻烦.扩展有限元法将裂缝独立于计算网格之外，利用扩充位移项来描述不连续的位移场从而表征裂缝特征.该方法模拟精度高，裂缝处理方便，但针对不同的裂纹问题需要选取不同的富集函数.Belytschko和Black[9]针对二维裂纹问题引入Westergaard函数作为富集函数(Enriched functions)；Moes等[10]则引入阶跃函数作为富集函数；对于一些尚无解析解的问题，Waisman和Belytschko[11]还提出了一种带未知参数的富集函数方法来提高计算精度.另外，裂缝条数过多时会引入附加自由度，计算速度慢，收敛性较差.本文采用有限元方法来进行裂纹扩展的数值模拟，无论采用何种模式，对于裂纹扩展数值模拟问题，裂尖应力场的分析均是关键问题.由于裂纹尖端是奇异点，所以在利用有限元法求解裂尖应力场时，其结果的准确性与单元类型和网格尺寸都密切相关，最好在特定位置采用奇异单元与足够小的单元尺寸.在模拟裂纹扩展的过程中完全满足此类要求是非常困难的，因为计算时需要跟踪裂纹几何位置的变化不断重构网格并变更单元类型.本文不需要直接面对有限元算不准的裂尖奇异应力场，而是通过裂纹尖端的应力强度因子衡量裂纹的开裂特性，并且它还可以判断出裂纹的开裂方向，避免了需要预先定义裂纹扩展路径的问题.应力强度因子是采用相互作用积分方法计算得到，而积分方法可以减轻计算结果对裂尖处单元类型与单元尺寸的依赖.该方法只需要通过远离裂尖处某一积分路径上的应力场解与位移场解，便可计算出裂尖的应力强度因子.远场处的有限元结果往往具有很高的精度，所以该方法能够获得较为满意的结果，这样就给裂纹扩展的分析带来了极大的方便.本文给出了进行裂纹扩展分析所需使用的策略，力求在计算含奇异点的裂纹扩展问题时寻求一种较为可靠的方法.为了提高可用与易用性，在通用有限元软件上进行简单的改进，实现了ANSYS软件的二次开发，计算了几个比较有代表性的断裂力学问题，比较单元杀死方法和扩展有限元的结果以及网格划分的影响，验证了本文方法的正确性与有效性.1 裂纹扩展的理论基础1.1 应力强度因子的计算应力强度因子常见的计算方法有相互作用积分方法、结点位移外推法[12]以及复变函数方法[13].本文的应力强度因子是通过相互作用积分方法得到的.然而，相互作用积分方法是基于准静态条件下的计算方法，所以无法求解裂纹扩展速度等动态参数.具体的基本理论如下[14].如图1所示，相互作用积分是包含裂纹尖端回路G的能量积分.它的定义如下[7]: I=(1)式中考虑了两个独立的平衡状态，状态为真实应力-变形场的变量，状态为附加应力-变形场的变量;为相互作用应变能密度;nj为积分回路上的外法线单位矢量.附加场可以是任何满足变形方程和本构方程的应力-变形场，一般情况下取为裂尖渐近场的解析解[15].图1 相互作用积分示意图Fig.1 Schematic diagram of interaction integration 当回路G 接近裂纹尖端时，相互作用积分与真实变形场和附加变形场的应力强度因子之间有如下关系:(2)式中E*为平面应变或平面应力下的弹性模量.在实际应用中，为方便有限元的计算，通常把相互作用积分转换成一个面积分，在回路G 的外面再设置另外一个回路C，两个回路之间所围成的区域为S，通过散度定理将式(1)的线积分转换为面积分，另外再考虑裂纹面上的面力作用，则相互作用积分就可以写为[16]I=(3)式中:q是平滑的权函数，在区域S内取值为1，其余地方为0;my是裂纹面上的外法线.通过以上分析可以知道，选取特定的辅助场就可以得到真实状态下的应力强度因子.为了求解参数可以假设辅助场为纯I型的状态，并且满足如下条件(4)由于纯I型状态下的应力场和位移场都有理论解，因此很容易通过式(3)得到相互作用积分的值，最后结合式(2)和式(4)就可以得到参数的表达式(5)同理假设辅助场为纯II型状态就可以得到参数的结果.1.2 裂纹扩展准则典型的断裂准则包括能量释放率、J积分和应力强度因子等相关准则.本文采用应力强度因子准则，在相互作用积分方法所求得裂纹尖端i应力强度因子和的基础上，定义等效应力强度因子[17]:(6)当裂尖的等效应力强度因子超过临界应力强度因子时，裂尖就开始扩展.对于复合型的裂纹扩展，应用最大周向拉应力强度因子理论来确定出裂纹扩展的方向角[18]: θi=(7)1.3 模拟裂纹扩展的过程在裂纹开裂过程的计算中，主要是通过计算确定出裂纹的起裂条件与裂纹的开裂角.另外一个问题就是裂纹扩展时模型网格的重新划分，重新划分的区域只需要选择所关心的扩展部分即可，因此只是增加了部分的结点自由度，详细的划分尺寸可以由下文的算例得到.由于利用相互作用积分方法求解裂纹尖端的应力强度因子，所以在网格划分时不需要采用特殊单元，这样就降低了网格重新划分的难度.具体的计算流程如下：1) 通过相互作用积分方法计算出裂纹尖端i处的应力强度因子和然后利用式(5)得到裂纹开裂方向角θi，最后根据复合裂纹断裂准则确定开裂点.2) 设置裂纹尖端的开裂步长为Δl，沿裂纹扩展方向角劈裂模型.关于开裂步长的设置，参考文献[19]中给出的建议增长步长为其中h为单元尺寸.3) 重新划分模型网格，需要特别注意裂纹尖端附近的网格尺寸，并防止重新划分的网格出现畸变.4) 重复过程1) ～3)，进行下一个开裂点的操作.如果开裂过程中不满足开裂准则那么裂纹扩展过程中止.2 典型问题的分析与比较在ANSYS有限元程序中，模型均采用八结点的PLANE183平面应变单元.分析接触问题时，接触单元采用三结点CONTA172单元，目标单元选择三结点TARGE169单元.在选择的这些典型问题中，将单元杀死方法、相互作用积分方法和扩展有限元方法进行了比较.首先比较了单元杀死方法和文章中的相互作用积分方法，目的是说明本文方法要比单元杀死方法计算更正确.其次比较了相互作用积分方法与扩展有限元方法，进一步证明了本文方法是有效可靠的.下面详细讨论这些算例.2.1 应力强度因子的计算应力强度因子的计算是本文方法中判断裂纹开裂的主要参数，其对模型网格的敏感性是首先必须要说明的问题.只有部分裂纹的应力强度因子有封闭解，所以本文选择有限宽中心裂纹板的解析解为例来说明应力强度因子的敏感性计算问题.该计算模型如图2所示，单向受均匀拉应力σ的有限尺寸板，宽为2W，高为2H，中心存在一条长为2a的裂纹.假设板为各向均匀的线弹性材料，裂纹尖端的应力强度因子表达式为[20]:(8)图2 含中心裂纹的有限大板Fig.2 The finite plate with a central crack其中F为有限宽度的修正系数，代表性系数包括以下3种.Isida近似系数：对a/W<0.7时误差为0.3%.修正的Koiter公式：误差为0.3 %.修正的Feddersen公式：误差为0.1 %.本算例中a/W的取值为0.2.对计算出来的应力强度因子进行无量纲化，记为Kn，Kn=KI/(σ表1中列出了不同的应力强度因子修正结果与本文计算得到的结果.从表1中可以看出，计算得到的结果与手册中的解析结果非常接近.表明了本文求解应力强度因子的方法是可靠的.表1 裂纹尖端的无量纲应力强度因子Tab.1 Normalized stress intensity factors at the crack tip求解方法Kn相对误差/%Isida方法1.02520.18Koiter方法1.0224-0.10Feddersen方法1.02440.10本文方法1.0234—通常在采用奇异单元求解应力强度因子的方案时，其计算结果对网格的尺寸有较高的要求[21]，这样在某些问题中就会因追求精度而降低了计算效率，并且不利于对裂纹扩展的计算.因此表2针对不同网格尺寸与类型所得到的结果进行了对比.结果表明网格的改变对最终结果的影响不大，并且对计算值的影响也很小，建议网格尺寸的取值为0.1a.由此可见网格划分对该方法的影响相比加密网格的传统有限元方法来说是不敏感的，并且不依赖于裂纹尖端近场应力的数据，所以本文方法具有一定的优越性.2.2 I型单裂纹的扩展模拟例证2为最简单的I型裂纹扩展问题.图3为该问题的模型简图，板的边缘有一条长度为I的初始裂纹，板的左右两侧承受σ的拉应力作用.最终的计算结果还与利用单元杀死方法的结果进行了对比并绘制出图4，试图通过这两种方法的对比说明本文方法预测裂纹扩展路径的可靠性.其中的实线为使用本文方法所计算出来的结果，而带三角符号的线(映射划分网格)、带四边形符号的线(自由划分网格1)和带十字星符号的线(自由划分网格2)均为单元杀死方法的结果.自由划分网格1与网格2的区别是网格尺寸的不同，网格1的尺寸是网格2尺寸的二倍.其余结果的网格尺寸与自由划分网格1的尺寸相同.表2 不同网格尺寸下裂纹尖端的无量纲应力强度因子Tab.2 Normalized stressintensity factors at the crack tip with different mesh size网格类型求解方法应力强度因子参考值0.02a0.05a0.10a0.20a计算值相对误差/%计算值相对误差/%计算值相对误差/%计算值相对误差/%三角形网格Isida方法1.02521.02450.071.02310.201.02340.181.0220.31Koiter方法1.02241.0245-0.211.0231-0.071.0234-0.101.0220.04Feddersen方法1.02441.0245-0.011.02310.131.02340.101.0220.23四边形网格Isida方法1.02521.02440.081.02240.271.02310.201.02150.36Koiter方法1.02241.0244-0.201.022401.0231-0.071.02150.09Feddersen方法1.02441.024401.02240.201.02310.131.02150.28注：a为裂纹半长.图3 I型裂纹拉伸的模型示意图Fig.3 Model of mode-I crack subjected to tension图4 I型裂纹拉伸的扩展结果图Fig.4 Propagation results of mode-I subjected to tension从图4中可以看出虽然各条曲线的趋势都是向前扩展，但是具体细节上还是存在差异.裂纹扩展过程中带四边形的线与带十字星的线都出现了偏折，这是因为单元杀死方法只能按照预先设定好的单元分布情况进行扩展，在自由划分的网格中往往不能保证网格线平直.然而带四边形的线与带十字星的线的扩展路径并不重合，这说明单元杀死方法在网格尺寸不同的情况下计算结果出现了差异.图中带三角形的线与实线基本重合，在映射网格下单元杀死方法可以捕捉到简单的裂纹开裂规律.虽然单元杀死方法可以预测出裂纹扩展的总体趋势，但是它严重依赖于网格划分的方式，并且由于网格划分的问题会使裂纹扩展路径不同，不利于后续问题的研究，因此不是一种稳定可靠的计算方法.利用本文方法计算出来的虚线不受网格划分的影响，才最接近真实的裂纹扩展结果.2.3 带孔板边缘多裂纹扩展模拟在验证完简单的裂纹扩展模型后，通过较复杂的例证对本文方法进行深度的校验.用来验证的模型示意图如图 5所示，旨在研究受圆孔缺陷的影响下裂纹的扩展规律.具体描述为含两个中心对称圆孔的矩形平板，底部固定，顶端承受预先给定的位移加载，使它们只能在竖直方向上有移动.两条长度均为1 mm的预制裂纹在孔洞离边界的较远一侧，并且呈中心对称分布.模型的材料参数如下：弹性模量E=200 GPa，泊松比ν=0.3，材料的临界应力强度因子KIC =1 500 N/mm3/2. 图5 带孔板边缘多裂纹扩展的模型示意图(单位：mm)Fig.5 Model of multiple edge cracks propagation with hole plate为了验证本文方法的可靠性，本文将通过改变单元类型的相互作用积分方法计算得到的裂纹扩展路径与文献中Khoei等人[22]的路径进行对比.裂纹的扩展均需要预先假定一个阈值，但是由于裂纹尖端具有奇异性，因此一般方法在处理裂纹尖端问题时对裂尖网格的划分要求比较高，即使单元类型相同，而尺寸的差异也会导致裂纹扩展路径的变化.依据之前介绍的计算步骤，本文方法所计算出来的裂纹扩展径图如图 6中带四边形实线与带三角形实线所示，而Khoei等人的计算结果则用实线标识.从图 6中可以看出本文的计算结果与Khoei的结果[22]都显示出相似的裂纹扩展路径.裂纹均因圆孔的影响而在开裂初期朝圆孔的方向扩展，在后期由于两条裂纹之间的相互作用而互相靠拢.由此可知本文的计算结果与文献的计算结果吻合较好，并且不同的网格划分方式对于路径的影响小.裂纹扩展不同阶段的网格划分图如图 7所示，裂纹每次向前扩展一定距离，在裂纹附近区域的网格都会进行重构，并且裂纹附近的网格尺寸比其他区域的尺寸更小，原因除了保证一定的精度外还因为能较好地保证网格的划分质量.图6 裂纹扩展路径对比图Fig.6 Graph of crack growth with different methods综上所述，本文利用的相互作用积分方法，对网格的划分方式并没有特别严格的要求，既适用于三角形网格又适用于四边形网格，而且扩展路径也与文献中的结果相差不大，因此是一种更可靠的方法.2.4 受压裂纹扩展的模拟在地层中有很多岩层都处于受压状态，然而其中的裂纹依然会产生错动而发生开裂.例如在水力压裂工程中，由于其具有围压的特殊性，天然裂纹会受水压裂纹以及地应力的影响而扩展，若在进行数值模拟时不考虑天然裂纹的压裂纹情况则会造成偏差.本文从这一角度出发，探讨利用同样的方法是否可以解决裂纹在受压条件下的扩展问题.图7 裂纹扩展不同阶段的单元划分图Fig.7 Element partition plot of crack growth path at different stage裂纹面在受压时会相互接触并且产生摩擦力作用，因此必须考虑接触力作用，而且为了防止裂纹面之间发生相互渗透，模拟时需要在裂纹面之间设置接触分析，一般来说就是用接触单元来模拟接触面的作用.其中较为常见的分析方法有罚函数法和增广型拉格朗日乘子法等.但是罚函数法存在收敛性差的问题[23]，因此本文采用的是增广型拉格朗日乘子法，通过判断双面约束来进行计算.计算过程中实时监测了接触面上的接触压力，保证接触面上不会产生破坏.模型采用与岩层类似的材料参数，具体数据为：弹性模量E =0.33 GPa，泊松比ν= 0.19，Keq=1.4N/mm3/2.模型的几何尺寸如图 8所示，模型长60 mm，高120 mm，裂纹长为a=12 mm，裂纹之间的竖向距离为b=20 mm，裂纹的倾斜角均为45°，在模型的底部和顶端均有压应力的作用.为了验证本文方法的正确性，将计算结果与石路杨等人[17]利用扩展有限元方法所得到的结果进行对比.图9(a)中的黑折线为本文方法计算出来的裂纹扩展路径，而图9(b)中的白折线为文献中的计算结果，可以看出两种方法所得到的结果基本一致.从图9中的计算路径可以得出，模型在受压状态下，裂纹的扩展过程始终是沿着模型的最大主应力进行的，这与初始裂纹考虑了摩擦力有关[18].除了初始裂纹面处于受压力作用外，扩展生成的新裂纹面是处于受拉作用的，因此这种破坏形式属于劈裂破坏.图 8 受压裂纹扩展的模型示意图Fig. 8 Model of crack subjected to compression(a)本文计算结果 (b)文献计算结果图9 裂纹扩展路径图Fig.9 Crack growth path如果受压的裂纹面之间没有摩擦，则裂纹的开裂路径则与带摩擦的情形有显著区别，角度大小与摩擦系数有关.裂纹的起裂方向不会沿着模型的最大主应力方向，在本例中与竖直方向有一个大约35°的夹角.在扩展一段时间后，裂纹才会逐渐转向并沿着最大主应力方向继续扩展.裂纹扩展路径图可详见图10，其中虚线代表无摩擦条件下的扩展路径，而实线则表示有摩擦的情况.该方法可以嵌入通用的有限元软件中，例如ANSYS，ABAQUS等.由于ANSYS擅长进行多物理场的分析，而水力压裂工程的模拟就是流固耦合的计算，因此基于本文所述的方法结合通用有限元平台进行流固耦合分析就较容易实现水力压裂中的裂纹扩展.这也是笔者今后要进一步深入研究的方向，因为依托大型通用有限元软件，所以计算方法也很容易推广.图10 无摩擦与有摩擦下的裂纹扩展路径图Fig.10 Crack growth paths under friction and frictionless condition3 结论本文提供了一种模拟裂纹扩展的计算方法.按照本文提出的计算流程与策略，便可以较正确地得到裂纹的扩展路径.利用相互作用积分方法来求解裂尖参数可以一定程度上消除奇异性影响，从而避免了裂纹尖端的奇异性对网格划分的苛刻要求.在传统有限元方法的基础上，通过相互作用积分法则求解裂纹尖端的应力强度因子，由最大周向拉应力理论来判断裂纹的扩展方向.模拟裂纹扩展的过程中，对涉及扩展区域的网格进行重新剖分并获得了较正确的扩展路径与较高的计算精度.应用该方法在模拟裂纹扩展的过程中，兼顾计算精度与效率的情况下建议网格尺寸取0.1倍的裂纹半长.该方法既能模拟裂纹拉伸时的情况也可以模拟受压时的扩展.相比裂纹受拉伸情况下的扩展，裂纹受压时的扩展需要设置接触，因此在计算时间上相对较长，但在网格尺寸方面并没有差别.在分析岩石破坏的过程中，由于材料存在较大的分散性，因此一般只关心岩层中裂纹的开裂方向，而对精度往往没有过高的要求.今后希望可以应用这种可靠性的方法来做进一步的研究.参考文献[1] 陶伟明, 班勇婷, 王慰军,等. 一种裂纹扩展的有限元仿真分析方法及其实现[J].力学季刊,2009, 30(4): 612-617.TAO Weiming,BAN Yongting,WANG Weijun,et al. A finite element simulation procedure and its realization for crack growth[J].Chinese Quarterly of Mechanics,2009,30(4):612-617. (In Chinese)[2] PORTELA A, ALIABADI M H, ROOKE D P. The dual boundary element method: effective implementation for crack problems[J]. International Journal for Numerical Methods in Engineering, 1992,33(6):1269-1287. [3] 雷钧, 黄祎丰, 杨庆生. 裂纹动态扩展的边界元模拟[J] . 北京工业大学学报, 2013, 39(6): 806-810.LEI Jun,HUANG Yifeng,YANG Qingsheng. Dynamic crack propagation by boundary element method[J].Journal of Beijing University of Technology,2013,39(6):806-810. (In Chinese)[4] 王强, 宋雪迪, 郝中华,等. 基于ABAQUS纤维梁单元的钢筋混凝土柱受力破坏全过程数值模拟[J]. 土木工程学报, 2014, 47(12):16-26.WANG Qiang, SONG Xuedi, HAO Zhonghua,et al.Numerical simulation of failure process of RC columns based on fiber beam element in ABAQUS[J]. China Civil Engineering Journal, 2014,47(12):16-26. (In Chinese)[5] BAO J Q, FATHI E, AMERI S. A coupled finite element method for the numerical simulation of hydraulic fracturing with a condensation technique[J]. Engineering Fracture Mechanics, 2014,131:269-281.[6] FU P, JOHNSON S M, CARRIGAN C R. An explicitly coupled hydro-geomechanical model for simulating hydraulic fracturing in arbitrary discrete fracture networks[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2012，37:278-300.[7] DOLBOW J, BELYTSCHKO T. A finite element method for crack growth without remeshing[J].International Journal for Numerical Methods in Engineering,1999,46:131-150.[8] 曾青冬, 姚军. 水平井多裂缝同步扩展数值模拟[J]. 石油学报，2015,36(12):1571-1579.ZENG Qingdong,YAO Jun.Numerical simulation of multiple fractures simultaneous propagation in horizontal well[J].Acta PetroleiSinica,2015,36(12):1571-1579. (In Chinese)[9] BELYTSCHKO T, BLACK T. Elastic crack growth in finite elements with minimal remeshing[J]. International Journal for Numerical Methods in Engineering, 1999, 45: 601-620.[10] MOES N, DOLBOW J, BELYTSCHKO T. A finite element method for crack growth without remeshing[J]. International Journal for NumericalMethods in Engineering, 1999, 46: 131-150.[11] WAISMAN H, BELYTSCHKO T. Parametric enrichment adaptivity by the extended finite element method[J]. International Journal for Numerical Methods in Engineering, 2008, 73: 1671-1692.[12] 肖汉宁, 彭苏华, 高朋召. 支柱瓷绝缘子表面裂纹应力强度因子的有限元分析[J]. 湖南大学学报(自然科学版), 2014, 41(12):53-58.XIAO Hanning, PENG Suhua, GAO Pengzhao. Stress intensity factor analysis of surface crack on rod porcelain insulator through finite element method[J]. Journal of Hunan University(Natural Sciences), 2014, 41(12):53-58. (In Chinese)[13] 肖万伸, 欧阳成, 张春雨. 一维六方准晶中界面裂纹与位错的交互作用[J]. 湖南大学学报(自然科学版), 2016, 43(12):71-75.XIAO Wanshen, OUYANG Cheng, ZHANG Chunyu. Interaction between circular interfacial crack and screw dislocation in 1d hexagonal quasicrystals[J]. Journal of Hunan University(Natural Sciences), 2016,43(12):71-75. (In Chinese)[14] 庄茁，柳占立，成斌斌,等. 扩展有限单元法[M]. 北京：清华大学出版社, 2012:20-21.ZHUANG Zhuo,LIU Zhanli,CHENG Binbin,et al.Extended finite element method[M].Beijing:Tsinghua University Press,2012:20-21.(In Chinese) [15] YAU J F, WANG S S, CORTEN H T. A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity[J]. Journal of Applied Mechanics, 1980, 47(2):335.[16] DOLBOW J, N, BELYTSCHKO T. An extended finite element method。

气-液-固三相流混合建模与求解方法*范兴华 谭大鹏† 李霖 殷梓超 王彤(浙江工业大学机械工程学院, 特种装备制造与先进加工技术教育部/浙江省重点实验室, 杭州　310014)(2020 年12 月14日收到; 2021 年1 月6日收到修改稿)气-液-固三相流混合过程是一个复杂的多重流固耦合动力学问题, 颗粒参数与流道物理空间尺度之间的关系直接影响计算收敛性, 强剪切区域的流固双向耦合作用数值建模与网格处理具有较高难度. 针对上述问题,提出了一种气-液-固三相流混合的建模与求解方法. 基于流体体积-离散单元耦合模型, 建立考虑颗粒运动的三相动力学模型, 通过求解动量方程, 实现两相流体与颗粒的双向耦合. 自主开发用户自定义函数(UDF)通信接口, 得到流体与颗粒间的相互作用力, 提出了一种多孔相间耦合解法来描述颗粒运动轨迹. 以带强剪切的三相流混合过程为例, 使用该方法研究了不同充气条件对流道物理空间内自由表面、速度分布和颗粒悬浮特性的影响规律. 结果表明, 强剪切和壁面作用可以将流体的切向速度转化为轴向和径向的速度; 选择合适的充气速度可以消除自由液面的不稳定性; 增加流体的流动速度, 对于部分区域颗粒的悬浮提升作用有限. 研究结果可为复杂多相流相间作用机理研究提供有益借鉴, 也可为气-液-固三相颗粒混合生产调控提供技术支持.关键词：三相流, 固体颗粒, 流体体积-离散单元耦合, 强剪切PACS：45.70.Mg, 47.57.–s, 47.57.E–, 64.75.Ef　DOI: 10.7498/aps.70.202021261 引　言气-液-固三相流混合是高端化工、锂电生产的关键工艺环节, 混合执行构件和流道物理空间需要提供较高的传质和高湍流能力, 且伴随强剪切过程. 上述要求使得气-液-固三相流混合过程非常复杂, 且难于观察整体流场及关键区域的颗粒分布[1−3]. 混合物理空间几何尺度相对于颗粒要高多个数量级, 其内部三维循环流动和湍流多相流的复杂性, 给混合执行机构优化设计、混合过程边界条件调控提出了重要技术挑战[4,5].当前计算流体力学(computational fluid dyna-mics, CFD)方法广泛应用于液-固混合过程的模拟计算, 一般基于欧拉-欧拉模型, 将颗粒固相视为连续相, 来描述相间相互渗透过程, 该模型占用计算资源比较少, 但模拟精度较低, 且无法获得颗粒运动状态[6−10]. 基于欧拉-拉格朗日模型的离散单元法(discrete element method, DEM)可以获得颗粒的运动和相互作用, 可以与不同的流体动力学计算方法相结合, 来模拟流体-颗粒流[11,12], 如格子-玻尔兹曼方法(lattice Boltzmann method, LBM).该方法在离散的晶格网格上使用代表流体相的虚拟颗粒, 并通过求解离散的Boltzmann方程模拟流体的流动[13]. 相关学者已经对三维LBM-DEM 耦合进行了尝试[14], 但是由于要求流体尺寸要比固体颗粒尺寸小得多, 对计算能力要求非常高, 多数针对三维问题的LBM-DEM耦合解法仍在开发中. 将CFD与离散单元法(CFD-DEM)耦合使用,可以预测颗粒尺度的变化, 颗粒-颗粒和颗粒-壁之间的相互作用都通过牛顿运动方程求解, 而颗粒-流体之间的相互作用则通过源相交换来实现. Bastien等[15]使用CFD-DEM模型, 以非常好的可靠性再现了固-液混合系统中发生的各种现象,* 国家自然科学基金(批准号: 51775501)和浙江省重点科学基金(批准号: LZ21E050003)资助的课题.† 通信作者. E-mail: tandapeng@© 2021 中国物理学会 Chinese Physical Society 研究证明, 在非惯性参考系下进行CFD-DEM模拟是可行的, 这为CFD-DEM的应用开辟了广阔的前景. Shao等[16]利用CFD-DEM耦合的模型研究了三维混合过程中的固体悬浮行为, 与实验测量和基于欧拉-欧拉方案的CFD仿真相比, CFD-DEM仿真可以提供更多流场信息. Blais等[17−19]开发了一种半解析CFD-DEM模型, 并利用该模型进行固-液混合操作的设计以及优化, 提高了颗粒悬浮比例, 改善了流型和颗粒分布, 但未考虑自由液体表面的分布及其稳定性.对于自由液体表面流, 常用流体体积(volume of fluid, VOF)模型进行模拟, Xu等[20]研究了表面涡流对混合物理空间中颗粒分散的影响, 研究发现表面涡流的产生降低了颗粒分布的均匀性, 混合挡流物理构件可使得表面涡流得到有效控制, 提高了颗粒分布的均匀性. Sun和Sakai[21]开发了一种基于欧拉-拉格朗日方法的模型, 可以模拟复杂的三相流行为, 包括自由表面的变形和颗粒引起的液体位移. Wu等[22]利用虚拟双重网格孔隙度模型改进了DEM-VOF模型, 新模型克服了计算过程中的不稳定性, 可用于固-液混合过程的流场模拟计算. Kang等[23]利用DEM-VOF模型结合雷诺压力模型(RSM), 对具有自由表面的无挡流构件的混合物理空间中的颗粒悬浮动力学进行了模拟,得到了混合流道的几何形状、叶轮转速、颗粒密度和直径等对自由表面涡流、流型和颗粒悬浮动力学的影响规律.当前关于多相流的研究主要以两相为主, 考虑三相耦合过程的多集中于自由液面, 却鲜有考虑内部充气对液-固两相的影响, 尤其是视固相为离散颗粒相的情况. 因此, 建立流体和颗粒的动力学模型, 通过求解动量方程, 实现流体与颗粒的双向耦合, 进行气-液-固三相流混合的研究, 揭示三相耦合在复杂混合过程中的作用规律非常必要.针对上述目标, 本文首先建立了气-液-固三相流动力学模型, 分别包括VOF模型、DEM模型以及两者的耦合模型; 然后进行了建模与网格划分以及边界条件和参数设置, 并进行了网格无关性验证, 进而选取最终使用的网格数量并开始不同案例的数值计算; 自主开发用户自定义函数(user defined function, UDF)通信接口, 得到流体与颗粒间的相互作用力, 提出了一种多孔相间耦合解法来描述颗粒运动轨迹; 最后通过计算结果讨论了充气对自由液面、流体速度以及颗粒悬浮的影响, 揭示了不同充气条件下混合物理空间内多相流的演变规律并得出了相应结论.2 三相流动力学模型对于VOF-DEM模型耦合的多相流, 通过对体积平均的Navier-Stokes方程进行求解, 进而对连续相的流体进行描述, 针对VOF模型还存在自由液面的问题, 气-液两相存在明显交界面, 可通过2.1节模型求解气-液两相问题. 同时使用离散单元法对固体相颗粒进行建模, 可通过2.2节模型求解. 这两个模型以一定的时间步间隔进行双向耦合交换数据, 一般CFD的时间步长明显大于DEM时间步长, 以便正确获得接触作用, 颗粒通常不会在单个DEM时间步长中移动很远的距离.因此, 不需要两者1∶1的时间步长比, 对于DEM, CFD的时间步长比一般从1∶10到1∶100不等, 两者耦合时选择合适的时间步长比, 不仅可以节省计算时间, 还可以避免计算的发散.2.1 VOF模型混合过程中存在复杂的气-液两相耦合现象,所以应该用多相流模型来描述. VOF模型是基于欧拉网格下的表面追踪模型, 通过求解单相或多相的体积分数来追踪和捕捉互不相融流体间的交界面. 通过计算各相的控制方程, 能够准确地模拟混合过程中多相组分的动态演化和瞬态特征的捕捉,其中流体的连续性方程和动量方程表示如下[21]:ρfεfµF pfF st式中, 是流体密度, 是空隙率, u是流体速度, p是压力, 是流体动力黏度, 是颗粒流体间的相互作用项的反作用力, 是自由液面附近的表面张力. 为了提高模拟精度, 使计算更加接近实际情况, 本文采用连续表面张力(continuum surface force, CSF)模型处理表面张力, 其表达式如下:σκ式中是流体的表面张力系数; 而是气-液相界面的曲率, 其表达式为n =∇·α2其中是次相体积分数的法向量.对于VOF 模型气-液两相的交界面可通过界面传输方程求解:αα=10<α<1α=0其中 是液体的体积分数, 若 , 代表该计算单元全是液相, 不含气相; 若 , 则说明该计算单元内同时包含气、液两相; 若 , 代表该相全部是气相. 因此, 基于VOF 模型可以求解气-液交界面的相互作用过程, 尤其是在强剪切作用下的气相剥离过程.k -ε此外, 多相流环境中尤其是强剪切区域处于湍流状态下, 为了能够得到精确的模拟结果, 流体控制方程选择标准的湍动能-耗散( )模型, 该模型在剧烈变化的流场中有较好的计算性能, 其控制方程如下[24]:εG k G b Y M C 1εC 2εC 3εi,j x,y,z σk σε式中, k 是湍动能; 是湍动能耗散率; 是由于平均速度梯度引起的湍动能产生项; 是由于浮力影响引起的湍动能产生项; 为可压缩湍流脉动膨胀对总的耗散率的影响; , , 为经验常数, 张量下标 表示 三个方向分量; ,为湍动能和湍动耗散率对应的普朗特数.2.2 DEM模型在考虑流体相时, 已经引入了颗粒与流体间的相互作用, 因此在该部分将分析颗粒模型. DEM 是一种可以用于计算非连续颗粒的运动规律, 并且可以分析离散颗粒接触力以及运动的分析模型. 此模型中, 将颗粒相视为离散相, 相比于其他方法更加贴近实际, 更能还原颗粒的真实运动情况, 颗粒的运动是基于牛顿第二定律的计算得出的, 通过计算可以得到颗粒的平移、旋转的速度和位置随时间的变化关系, 颗粒的平移运动取决于作用在其上的力的总和, 而旋转运动则由接触的转矩控制, 其控制方程可表示为[16,25]m i x i F c ,ij F pf ,i I i θp ,i T t ,ij T r ,ij 式中, 是颗粒i 的质量; 是颗粒的位移; 是颗粒i 和j 之间的接触力; 是颗粒i 与流体间的相互作用力; g 是重力加速度; 是颗粒i 的惯性力矩; 是颗粒i 的角位移; 和 是作用在颗粒上的切向和滚动摩擦力矩, (8)式和(9)式中包含的一些力和扭矩的详细计算见如下公式.接触力:法向接触力:.切向接触力:切向摩擦力矩:滚动摩擦力矩:颗粒-流体作用力:Y ∗R ∗δcn ,ij δct ,ij S cn ,ij S ct ,ij m ∗v cn v ct L ij n ij µr ωij ∆V n p 这里 表示接触颗粒间等效杨氏模量; 表示接触颗粒间等效半径; , 分别表示接触颗粒间法向、切向重叠大小; , 分别表示接触颗粒间法向、切向接触刚度; 表示接触颗粒间等效质量; , 分别表示接触颗粒间的法向、切向相对速度; 表示接触颗粒间中心距离; 表示接触颗粒间单位矢量; 表示颗粒的滚动摩擦系数; 表示颗粒接触平面上的角速度矢量; 是颗粒i 所在的计算网格的体积; 表示颗粒数量;f pf ,i f d ,i f ∇p,i f ∇τ,i f st ,i f vm ,i f B ,i f Saff ,i f Mag ,i F pf 表示颗粒i 受到的所有流体、固体对颗粒作用之和, 包括曳力 [26,27]、压力梯度力 、黏性力、表面张力 、虚拟质量力 、Basset 力 、Saffman 升力 [28]和Magnus 升力 [23].由于动量方程表达式已经直接包含了压力、黏性力和表面张力, 所以在 中相应地把这些力减去了.2.3 耦合模型为了获得流体与颗粒间准确的相互作用力, 提出了一种相间耦合解法—多孔模型来描述精确的颗粒运动轨迹, 其计算表达式如下:εps ,i 式中, 代表流体单元内第i 个多孔球单位体积内的颗粒体积. 该模型克服了传统方法中当颗粒粒径接近网格尺寸时引起的计算不稳定性问题, 改善了颗粒-流体的相互作用, 并且在计算过程中, 把颗粒的体积考虑在内, 因此, 用该方法求解出的流场也更加精确.将2.1节和2.2节两个模型以及多孔模型的控制方程写入接口程序, 进行编译, 最终通过动量方程中的动量源交换项实现双向耦合, 这样VOF-DEM 耦合通过编好的用户自定义函数(UDF)进行数据通信, 实现欧拉双流体相和拉格朗日颗粒相的双向耦合.整体模型的计算流程如上图1所示. 首先对流场和颗粒场进行初始化, 该过程通过CFD 计算接口文件实现; 然后, 开始计算, 通过2.1节(1)式和(2)式迭代求解得到流场的速度和力等信息, 求解(5)式得到关于自由液面的变化信息, 通过(6)式和(7)式计算流体的湍动能-耗散; DEM 通过利用2.2节的(8)式和(9)式迭代计算获得颗粒的速度和位置等信息, 并进行更新; 接着通过判断收敛与否, 进行选择, 如果不收敛通过求解2.3节(17)式得到流体单元的空隙率, 继续前述流场的计算, 如此循环实现双向耦合, 互相交换数据, 直到收敛停止计算, 完成模拟.3 数值计算3.1 物理模型与网格划分研究所选取的物理混合空间为带挡流板的半圆底容器, 混合执行构件为叶轮, 结构如图2所示,具体物理参数如下: 直径T = 200 mm, 高度H =3T /2, 液位高度h l = T , 叶轮直径D = T /2, 桨叶长度L = 45 mm, 宽度W = T /10, 厚度t = 2 mm,倾斜角为45°, 安装高度C = 93 mm, 搅拌轴直径d 1 = 14 mm, 底部进气口直径d 2 = 14 mm, 高度h a = 4 mm, 挡流板高度h b = 11 T /10, 宽度为T /10.图 2 混合空间结构示意图Fig. 2. Diagram of mixed space structure.首先建立混合物理空间三维模型, 对流体域进行网格划分, 最终生成的网格如图3所示. 为了方便计算, 流体域被划分为包含混合叶轮的转子区域和除此以外的静子区域两部分. 然后对流体域进行离散化处理, 由于转子区域的变化梯度较大, 尤其是混合叶轮这种小尺寸, 强剪切区域变化梯度最大, 划分时要特别注意, 所以要进行局部网格的加密处理, 划分后的静子区域网格尺寸为7 mm, 转子区域网格尺寸为5 mm, 叶面网格尺寸进一步细化为3 mm, 进气口也进行适当的加密, 网格尺寸为3 mm, 划分网格后, 网格的正交质量保持在0.5以上, 最终用于计算的网格数目为31万.图 1 VOF-DEM 耦合计算流程图Fig. 1. VOF-DEM coupling calculation flowchart.3.2 边界条件及参数选择k -ε混合过程模拟采用了瞬态的VOF 模型计算,选择显式的时间离散格式, 湍流模型使用标准的湍动能-耗散( )模型, 该模型具有较高的物理可靠性, 可为混合湍流过程提供精确解, 近壁区域选择标准壁面函数. 混合容器顶部设置为压力出口边界条件, 容器壁面为无滑移壁面边界条件, 底部为充气管道, 并采用速度入口边界条件. 数值求解方法使用coupled 方案, 该解法耦合了动量、能量及组分方程, 能比较快地得到收敛解, 动量离散格式、湍动能和湍动能耗散率离散格式均采用二阶迎风以获得精确的解[29−31], 体积分数离散格式使用分段线性界面重构(piecewise linear interface construction, PLIC)算法, 这种方法是精度最高的一种[24,29], 监视器收敛残差均为10–6.叶轮旋转模型常用到滑移网格(sliding-grid,SG)方法和多重参考系(multiple reference frames,MRF)方法. 其中, SG 方法常用于瞬态模拟, 而MRF 方法通常用于稳态模拟, 文献[32]采用两种方法对比得到的最终结果非常相似. MRF 方法也能够用于瞬态仿真, 此种情况是以伪稳态方式进行计算, 与更准确、但更耗时的SG 方法相比, 可节省大量计算机资源, 精度能满足多数场景的需要[33−35]. 因此本次模拟采用瞬态MRF 方法进行.对于MRF 方法, 需要将流体区域划分为内部动区域和外部静止区域两部分, 两部分通过交界面(interface)进行数据传递. 对于本研究所涉及的其他物理参数设置, 见表1.3.3 网格无关性分析网格的大小会直接影响数值模拟的结果, 一般情况下, 网格越小, 计算结果也越精确, 但是, 随之带来的是网格数量也越来越多, 需要更多的计算时间才能收敛, 给计算带来了很大的困难. 因此, 有必要进行网格独立性研究, 以确保计算误差在可接受的范围内, 得到准确的模拟结果. 本次以倾角45°的桨式叶轮进行网格无关性验证, 模拟转速为400 r/min 时的混合流场, 使用了四种网格尺寸,总网格数分别为219000, 311000, 425000, 661000,考察网格数量对模拟结果的影响, 对5 s 时槽内Z = 150 mm, Y = 0 mm, X 从–100到100 mm 的轴向速度场进行对比, 结果如图4所示.可以发现流场在不同位置的轴向速度具有相似的趋势, 但网格数量为219000时整体具有较大的误差, 而其他三种网格计算误差在5%以内, 满足网格独立性要求. 基于上述结果, 后面计算模型所采用的网格数量均为311000, 这样既减少了计算量, 又可以得到比较准确的结果.表 1 流体和颗粒特性设置Table 1. Characteristics settings of fluid and particle.参数值ρa kg ·m −3空气密度 /( )1µa Pa ·s 空气黏度 /( )1 × 10–5ρw kg ·m −3水密度 /( )1000µw Pa ·s 水黏度 /( )0.001ρp kg ·m −3颗粒密度 /( )1100d p 颗粒直径 /mm 1n p颗粒数目 10000Y P MPa颗粒杨氏模量 / 1νP 颗粒泊松比 0.25Y w MPa壁面杨氏模量 / 70000νw 壁面泊松比 0.3µr 滚动摩擦系数 0.01µs 静摩擦系数 0.5e r恢复系数 0.5ωr ·min −1搅拌桨速度 /( )400∆t CFD s CFD 时间步 / 2×10−4 ∆t DEM s DEM 时间步 / 2×10−5 ∆t coupling s耦合时间步 / 2×10−4(a)(c)(b)图 3 网格划分　(a) 静子区域网格; (b) 转子区域网格;(c) 叶轮网格Fig. 3. Grids division: (a) Grids of stator region; (b) grids of rotor region; (c) grids of impeller.4 结果与讨论如前所述, 气-液-固三相混合物理空间内部是一个复杂的湍流环境, 挡流构件的存在增加了搅拌速度下湍流场的无序性和非线性, 为了获得其中的流体-固体多相耦合和相间传质特性, 将深入研究对比不同充气条件下对混合空间内自由表面、流体流动和固体颗粒悬浮的影响. 在数值算例中, 颗粒直径均为1 mm, 颗粒数目均为10000, 且颗粒随机分布在混合容器的底部区域内, 在初始条件下仅受重力作用.4.1 对自由液面的影响基于数值模拟的结果, 首先研究充气状态对自由液面的影响, 图5为t = 5 s 时四种充气条件下混合空间内的实际自由液面图, 蓝色为自由液面的演化形态, 分别对应充气速度为0, 0.01, 0.05和1 m/s. 为了能够清晰地看出自由液面的变化, 隐藏了底部液体和颗粒, 只保留了自由液面.通过图5(a)—(d)可以看出, 在图5(a)未充气的混合空间中液面存在小幅的波动, 因为底部搅拌对流场的扰动能量向上传输至自由液面时, 能量的衰减不足以对液面造成较大的扰动, 而这种轻微的波动主要是颗粒和液体的相互作用造成的; 对比之下, 充气速度v = 0.01和v = 0.05 m/s 时, 如图5(b)和图5(c)时两者液面波动都很小, 显然传输的湍动能依然达不到自由液面的扰动阈值; 而当充气速度v = 1 m/s 时, 其自由液面如图5(d)所示, 可以看出液面波动非常大, 尤其是挡流板之间的区域, 液面上升明显, 有漩涡的产生, 这主要是因为除了颗粒对液面的影响之外, 底部充气速度较大, 增加了湍流场的流体上冲动能, 对三相流系统造成了较大的扰动, 且气体上浮溢出造成了液面的不稳定, 进而有较大的振荡.针对上述充气扰动液面的现象, 考虑与搅拌下流场的流动模式密切相关, 图6给出了不同充气条件下的切向速度矢量图. 从图6可以看出, 四种条件下挡流构件附近的切向速度都是该截面内最小的, 是因为挡流构件将流体的切向速度转换为了轴向速度. 而当高速旋转的流体与挡流构件接触时,快速接触过程必然引起局部湍流涡团的耗散, 增加了流动模式的随机性, 故挡流构件附近的液面振荡相对明显. 同时, 通过对比可以发现, v = 1 m/s 时的流体切向速度分布极不均匀并且最大, 较强的充气强度对底部分布的颗粒相起到扰动作用, 在搅拌速度和底吹作用下, 整个流场处于非线性湍流状态, 这也是引起自由液面漩涡和波动的主要原因; 而其他三种状态下的切向速度则相对较小, 其中 v = 0.05 m/s 时的切向速度最均匀, 相对比较A x i a l v e l o c i t y /(m S s -1)Position/m图 4 四种网格尺寸在t = 5 s 时的轴向速度分布Fig. 4. Axial velocity distribution of four grid sizes at t = 5 s.(a)(b)(c)(d)图 5 t = 5 s 时四种充气条件下混合空间内的自由液面　(a) v = 0 m/s; (b) v = 0.01 m/s; (c) v = 0.05 m/s;(d) v = 1 m/sFig. 5. Free surface under four aeration conditions at t =5 s: (a) v = 0 m/s; (b) v = 0.01 m/s; (c) v = 0.05 m/s;(d) v = 1 m/s.稳定, 其对应的自由液面也是最稳定的, 这也与文献[23]的结论相吻合.通过与文献[22]关于自由液面的结果的对比可以发现, 本研究的混合容器虽然加装了挡流构件, 在不充气或者低充气速度时, 可以有效地消除混合过程中可能形成的漩涡的影响, 稳定混合空间内环境, 但是当充气速度过高时, 内部流场受到强烈扰动, 自由液面也会有大幅波动.4.2 对流体速度的影响为了研究不同的充气状态对混合空间内流体速度的影响, 选取了t = 5 s时, 混合容器内径向位置r = 60 mm, 轴向高度从0到200 mm的速度分布情况, 如图7所示. 从图7可以看出, 最大速度出现在靠近叶轮的区域, 底部充气会使得出现最大速度的高度上移, 但并不是充气速度越大最大速度出现的位置越靠上, 四种流体速度分布出现最大速度的高度由低到高依次是充气速度为v = 0, v = 0.05, v = 0.01和v = 1 m/s时的混合空间.上述现象主要原因为当轴向的充气速度介入流场时, 轴向动能和流场的切向动能发生了对冲能量耗散, 虽然流场的湍流随机性增加, 但流动的无序性对整个流场的总速度有所影响.2.0(a)(b)(c)(d)Tangentialvelocity/1-11.51.00.5图 6 t = 5 s时, 四种充气条件下混合空间内z = 0.15 m高度截面的切向速度矢量图　(a) v = 0 m/s; (b) v = 0.01 m/s;(c) v = 0.05 m/s; (d) v = 1 m/sFig. 6. Tangential velocity vector in height z = 0.15 m un-der four aeration conditions at t = 5 s: (a) v = 0 m/s; (b) v =0.01 m/s; (c) v = 0.05 m/s; (d) v = 1 m/s.Velocity/(mSs-1)Axialvelocity/(mSs-1)Axial height/m Axial height/mRadialvelocity/(mSs-1)Axial height/mTangentialvelocity/(mSs-1)Axial height/m图 7 径向位置r = 60 mm处的轴向高度速度分布　(a) 总速度; (b) 轴向速度; (c) 径向速度; (d) 切向速度Fig. 7. Axial velocity distribution at radial position r = 60 mm: (a) Total velocity; (b) axial velocity; (c) radial velocity; (d) tangen-tial velocity.从总速度、轴向速度和径向速度可以看出四种充气状态下的速度分布曲线趋势走向是相同的, 而对于切向速度, 在充气速度为v = 1 m/s 的桨叶下方的流体切向速度方向和其他三种情况方向完全相反. 显然, 当充气速度足够大时, 气相对整个流场的流动模式造成较大的扰动, 剪切流动变得复杂, 增加了流场的湍流混沌特性. 此外还可以看出,由于挡流构件的存在以及壁面的影响, 可以将切向速度转化为轴向速度和径向速度, 所以从图7可以看出切向速度的幅度要小于轴向速度和径向速度.4.3 对颗粒悬浮的影响图8、图9、图10和图11是四种充气状态下混合空间内三相流的模拟计算结果, 截取了部分时刻的运动状态, 这些时刻基本包含了混合过程中的所有情形, 具有一定的代表性, 其中颗粒的颜色表示颗粒的速度大小. 可以很清楚地看到在t = 1 s 时, 底部沉积的颗粒在桨叶旋转带动流体的作用下被卷吸起来. 初始状态时, 搅拌扰动流场, 增加了流场的切向流动和叶轮底部的轴向上升流运动, 颗粒以较小的速度向上升起.在t = 1.3 s 时, 颗粒群到达叶轮, 受到高速旋转的作用, 被桨叶打散高速向周围扩散, 同时可看出, 颗粒群在v = 0.01, 0.05和1 m/s 要先于v = 0情况分散开. 这是因为底部吹气在初始状态就增大了流场的轴向剪切流动能量, 诱导流体的轴向速度=1.0 s =1.3 s =1.5 s =1.7 s =2.0 s=3.0 s1.00P a r t i c l e v e l o c i t y /(m S s -1)0.750.500.250图 8 v = 0 m/s 时不同时刻的三相混合系统模拟结果Fig. 8. Simulation results of three-phase stirred system at different time when v = 0 m/s.=1.0 s =1.3 s =1.5 s =1.7 s =2.0 s=3.0 s1.00P a r t i c l e v e l o c i t y /(m S s -1)0.750.500.25图 9 v = 0.01 m/s 时不同时刻的三相混合系统模拟结果Fig. 9. Simulation results of three-phase stirred system at different time when v = 0.01 m/s.=1.0 s =1.3 s =1.5 s =1.7 s =2.0 s=3.0 s1.00P a r t i c l e v e l o c i t y /(m S s -1)0.750.500.250图 10 v = 0.05 m/s 时不同时刻的三相混合系统模拟结果Fig. 10. Simulation results of three-phase stirred system at different time when v = 0.05 m/s.增加, v = 1 m/s 时最为明显, 如图7(b)所示. t =1.5 s 时, 颗粒受到流场离心力的作用下扩散到容器壁面和挡流构件, 在两者的作用下, 颗粒流动转变为上下两个方向的运动状态, 表明当旋转流场与挡流构件接触过程时, 流场以切向为主的流动模式遇阻, 切向速度流动转变为上下剪切流动模式. 剩余流场的湍动能推动部分颗粒作上升流运动, 另一部分颗粒在重力作用向下流动, 壁面附近的湍动能难以对这部分颗粒提供足够的驱动动能.随着时间的推移, 在t = 1.7 s 时, 颗粒到达自由液面附近, 受到液面振荡的作用, 分散至容器的各个位置; 随后的t = 2和3 s 时, 除v = 1 m/s 的混合空间外, 其他变化已经不明显, 趋于稳态, 可以看出, 速度最大的颗粒分布在叶轮附近, 具有较高速度的颗粒分布在容器的下部, 而越靠近液面处, 颗粒速度越低, 到达液面时速度几乎为0, 此区域的颗粒群对气-液交界面的冲击非常微小, 颗粒停止上升, 只有挡流构件附近的颗粒群受到局部流体漩涡的裹挟作用, 会继续上升, 形成局部的凸液面. 此外, 上下两股流动实现了循环, 对颗粒分散效果有良好的促进作用, 且除v = 1 m/s 外, 其他都相对稳定, 颗粒分布都是对称的.从整个时间段上三相流的模拟计算结果对比,还可以看出底部充气使得下部伞状颗粒群伞柄处要比未充气状态的细, 而v = 1 m/s 工况下, 随着混合的进行颗粒受到的作用不均匀, 底部伞状颗粒群会被破坏, 运动状态也更加复杂, 液面波动也会更加剧烈, 难以实现稳态.为了更准确、更直观地描述颗粒分布的均匀性, 引入相对标准偏差(relative standard devia-tion, RSD)表征颗粒在液相中的分散效果[36], 在这项工作中, 将液体覆盖的区域分为12个部分(3 ×1 × 4), 即12个等体积的采样空间. 通过多个样本空间内颗粒数目的相对标准偏差随时间的变化作为评价指标. 其评价公式如下:avgX i X avg X i 其中 是采样空间i 中的颗粒数量, 是 的平均值, n 为划分的采样空间的数量. 从上面的相对标准偏差RSD 评价公式可以看出, RSD 越小代表颗粒在计算区域内的颗粒悬浮效果越好. 图12绘制了底部不同充气工作条件下的RSD 随时间的变化曲线.从图12可以看出, 在混合容器刚开始工作的一段时间内由于颗粒沉积在底部, 随时间移动较缓慢, RSD 的数值较大, 颗粒分布极不均匀,所以RSD曲线有一段上升趋势,这段时间里底部充气的混合空间由于气体的作用使得一部分颗粒上升较快, 相对未充气的颗粒率先上浮,所以RSD 值要低一些. 随着混合过程的持续进行, 颗粒到达叶轮高度, 在叶轮的作用下移动速度加快, 逐渐分散=1.0 s =1.3 s =1.5 s=1.7 s=2.0 s=3.0 s1.00P a r t i c l e v e l o c i t y /(m S s -1)0.750.500.250图 11 v = 1 m/s 时不同时刻的三相混合系统模拟结果Fig. 11. Simulation results of three-phase stirred system at different time when v = 1 m/s.R S DTime/s图 12 不同充气工作条件下的RSD 随时间的变化Fig. 12. RSD changes with time under different aeration conditions.。