数值模拟计算（精选十篇）
数值模拟计算 篇1
现阶段我国水电建设主要集中在西部地区,大多属于山区河流,两岸地形陡峻,山岩坚实,采用隧洞导流较为普遍[1],导流隧洞的水力参数、水力特性以及导流隧洞的优化设计方案一般都是通过物理模型试验得出,但模型试验难以全面量测洞内流速、压力分布,存在缩尺效应等问题,同时在体型优化时耗时、耗力。随着计算机性能的提高和数值计算方法的改进,计算流体力学在水利工程中的应用越来越广泛,大量研究表明[2,3],Fluent软件在计算外形比较规则的有压管道、隧洞以及明渠时误差较小,计算结果能够满足实际工程的要求。本文用Fluent软件对某具体工程进行了数值模拟计算研究。该工程由两条导流洞联合泄流,但由于考虑到同时计算两条导流洞时划分的网格数量很庞大,计算困难,并且两条导流洞具有一定的对称性,所以只对一条导流洞进行模拟计算。
1 数学模型
1.1 控制方程
k-ε双方程紊流模型是模拟强紊动水流的有效模型。对于有强旋流或带有弯曲壁面的流动,标准k-ε模型会出现一定的失真,而RNGk-ε模型有更好的适用性[4]。本文采用RNGk-ε模型,自由面追踪采用广泛应用的VOF方法。不可压缩非定常流的张量形式的控制方程如下:
连续方程:
VOF法[5]中在控制体内对第q相流体的容积分数规定为:αq=0表示控制体内无q相流体;αq=1 表示控制体内充满q相流体;0<αq<1表示控制体内部分充满q相流体,对所有流体的容积分数总和为1,即∑αq=1。每个控制体内混合流体的密度可表示为:ρ=∑αqρq;每个控制体内混合流体的黏性可表示为:μ=∑αqμq。
1.2 网格划分
由于主要关注导流洞内的水力特性,导流洞的进口段、出口段简化成一段明渠,计算区域包括上、下游明渠段、闸室段和洞身段。上下游明渠段和洞身段形状规则,采用结构化六面体网格,闸室段形状复杂,采用结构化与非结构化网格相结合的网格划分方式。研究表明[6],当网格间距小于0.6 m时,计算结果与网格无关,而网格间距大于0.6 m时,计算结果依赖于网格间距。综合考虑计算时间和计算精度,算例的横向断面网格间距采用0.6 m,纵向网格间距采用0.7 m。
计算区域以及网格划分见图1。
1.3 求解设置
边界条件水入口采用速度入口,空气入口采用压强入口,出口边界采用压强出口,其余都按壁面边界条件处理。湍流模型选择RNGk-ε模型,操作压力为101 325 Pa,压力速度耦合选用PISO[7], PISO算法比SIMPLE和SIMPLEC算法多了一个修正步,包含一个预测步和两个修正步,目的是使速度和压力值更好地同时满足动量方程和连续方程,以加快单个迭代步中的收敛速度。在固壁上给定法向速度为零和无滑移条件,近壁黏性底层采用标准壁面函数,用非恒定流进行计算。
2 物理模型试验
模型按重力相似准则设计,并考虑阻力相似,模型几何比尺λL=50。为了量测导流洞在有压流工况下的洞身压力,两条导流洞的洞身都布置了侧压管。其中,右岸导流洞测压管共布置了9个断面,每个断面布置四个测点,如图3、图4。
3 数值计算结果与模型试验成功的对比分析
3.1 洞身压力分布
图5为几个典型断面的压力云图。
图6(a)为导流洞各断面顶部数值计算压力和模型试验实测压力沿程分布图,通过分析对比可知,断面2至断面6的计算值与试验测试值基本一致,断面7和断面8两值相差较大,其原因为,导流洞出口为自由出流,存在自然掺气现象,而数值模拟计算无法给出自然掺气的边界条件。
图6(b)至图6(d)分别为导流洞各断面右侧中部、底部和右侧中部计算压力和试验实测压力沿程分布图,可见计算值和试验值吻合较好,从断面6到断面8压力值可以看出,由于受离心力影响,导流洞底部和右侧在弯道中部的压力比弯道起点和端点的压力高,而左侧刚好相反,弯道的内侧压力比外侧压力低。
3.2 流速、流态分布
从图7的速度云图可以看出,由于在闸室段分成了两个空腔进水,所以在导流洞的入口处形成了两股水流,最大流速集中在导流洞边壁上,最大流速达16.7 m/s,经过一段距离两股水流开始往中部靠拢并逐步融合,当到达弯道的某个断面时最终汇合成一股水流,且速度主要集中在弯道内侧附近。
4 结 语
通过对实际工程的模拟计算,把计算结果与模型试验结果分析对比,得出以下结论:
(1)采用数值模拟的手段,得出导流洞沿程压力分布、流速等水力特性的结果与模型试验得出的结果基本上吻合。
(2)从数值计算和模型试验得出的压力分布情况可看出,负压一般都只发生在导流洞的顶部,而底部和侧部一般都不发生。
(3)导流洞弯段压力分布特点是,弯道底部和外侧中部在弯道中部的压力比弯道起点和端点的压力高,而弯道内侧刚好相反,弯道的内侧压力比外侧压力低。弯道流速的分布特点为内侧的速度比外侧高。
摘要：物理模型试验是研究导流洞水力特性最普遍的方法,但存在缩尺效应等缺点。为了探索采用数值模拟的方法来研究导流洞水力特性的可行性,用Fluent商用计算软件,对某工程的单条导流洞进行了数值计算,得出了压强沿程分布、流速和流态等水力特性,选取典型断面的计算结果与模型试验结果进行对比,两者符合较好,验证了数值计算方法的可行性。数值计算方法可作为模型试验的辅助手段,并为Fluent软件在导流洞体型优化设计中的应用提供了参考依据。
关键词：导流洞,数值模拟,水力特性,RNGk-ε模型
参考文献
[1]袁光裕,胡志根.水利工程施工[M].第四版.北京:中国水利水电出版社,2005:11-12.
[2]邓军,许唯临,雷军,等.高水头岸边泄洪洞水力特性的数值模拟[J].水力学报,2005,36(10):1209-1212.
[3]周玉国,杨建东.尾水洞明满混合流的数值模拟计算[J].中国农村水利水电,2009,(12):123-126.
[4]王福军.计算流体动力学分析-CFD软件原理与应用[M].北京:清华大学出版社,2004.
[5]陈大宏,李炜.自由表面流动数值模拟方法的探讨[J].水动力学研究与进展,2001,16(2):216-224.
[6]杨晓池,张晓莉,杜淑英,等.FLUENT软件在泄洪洞体型优化中的应用[J].西北水电,2010,(2):82-84.
数值模拟计算 篇2
1大学物理教学中融入数值计算与模拟技术的必要性
随着计算机技术的发展，在大学物理教学中引入数值计算与模拟技术是计算机辅助大学物理教学的一种新形式，不失为一种比较好的教学手段。这种教学手段的优点在于:一方面它可以对物理问题进行数值计算求解，使得许多没有解析解的物理问题通过计算机求数值解而得到解决[1];另一方面，它还可以对物理问题进行模拟仿真，输出的仿真图像直观、清晰、形象、生动和真实，既可以帮助教师节约板书绘图时间，而且学生看过之后对物理知识的理解更加深刻，这种仿真还可以随意更改仿真参数，输出不同条件下的仿真图像，帮助学生全面理解所学物理知识。当前，计算机数值计算和模拟技术已日益广泛地应用于设计规划、生产制造和科学研究等各个方面，这就使得高校在校学生尤其是理工科专业学生需要具备一定的数值计算与模拟能力，以便更好地适应社会的需求。国内就有专家建议将数值计算与模拟能力写进理工科各专业培养计划中[2]，以便引起高校对理工科本科生数值计算与模拟能力培养的重视。在大学物理教学中，教师可以指导学生对物理问题进行数值计算与模拟，培养学生利用计算机解决实际物理问题的能力，使学生在大学低年级就受到这方面的教育，为学生在高年级专业课学习以及大四进行毕业设计打下良好数值计算与模拟基础。教育部曾多次发文强调加强实践教学，切实提高大学生的实践与创新能力，学生的数值计算与模拟能力无疑是其中重要一环。
2大学物理教学中融入数值计算与模拟技术的可行性
首先，教学硬件上有保障。高校一般都建设了一大批多媒体教室，教室电脑里一般都安装了各种数值计算与模拟软件，教师在大学物理教学中应用数值计算与模拟技术来辅助教学是没有问题的。同时高校一般都建设了全校性的计算机中心，理工科各专业还都建设了各自的电脑机房，学生可以在全校计算机中心或专业机房里应用计算机对研究的问题进行数值计算与模拟。此外，很多学生在大学低年级就配备了电脑，他们在宿舍里就可以运用计算机对相关问题进行数值计算与模拟。其次，教学软件上也有保障。许多大学物理教师越来越意识到计算机辅助大学物理教学的重要性，因此不断提高自己的计算机操作技能，不断学习和掌握各种办公软件，一般都熟练掌握了一、两门数值计算软件。学生在上大学物理课前，一般都学习了《计算机基础》这样的全校性基础课程，具备了一定的计算机操作能力，学习了一些常用的办公软件和数值计算软件。
3大学物理教学中融入数值计算与模拟技术的实施途径探索
3.1选择比较复杂的物理问题进行数值求解
大学物理中有些问题是比较复杂的，用简单的工具(如计算器)求解不仅耗时间，而且不一定能解出来，比如描绘麦克斯韦气体速率分布曲线，其公式f(v)=4π(m2πkT)32v2e-mv22kT非常复杂，用传统方法来描绘曲线，一般先用计算器计算出各点的值(v，f(v))，然后在作图纸上描出各点，最后将各点连成曲线，由于计算量大，花费的时间长，并且因为所取的样点不多，连成的曲线不光滑，失真度高。但是借助于计算机的数值计算，这个问题就迎刃而解了。笔者运用数值计算软件matalb[3]，先编写求解该问题的程序，程序中对自变量的步长取得很小，这样获得的样点数就很多，然后在计算机上运行程序，最后计算机输出麦克斯韦气体速率分布曲线，因为所取的样点数很多，所以计算机生成的曲线非常光滑，而且由于现在的计算机运算速度都比较快，生成曲线的时间都很短。此外，公式中的参量、也会影响曲线的分布，如果是用传统的方法，又需要重新计算，而采用计算机数值求解，只需在程序中改变参量的值即可，如图1所示。
3.2对需要可视化的.物理问题进行模拟仿真
大学物理中有些教学内容是需要可视化的，比如波动光学问题。一般教师在讲授光学内容时，有的老师直接在黑板上画图，费时又费力，画出的图像很难表现出光的明暗变化，有的老师借助于课件来展示，但是课件里的光学图像一般都是普通的画图软件制作的，跟真实的光学图像有很大距离。但是借助于数值计算与模拟软件，教师很容易对光学问题进行模拟仿真，仿真出的图像清晰、逼真。例如，采用数值计算与模拟软件matlab，笔者模拟了大学物理中的杨氏双缝干涉[4]，如图2所示，可以看出，模拟的双缝干涉图像清晰、逼真，同时还给出了光强分布曲线进行对照，学生很容易从仿真图像理解双缝干涉的特点，加深学生对双缝干涉知识的理解，这是传统的教师板书和一般的计算机辅助教学所不可比拟的。再例如，光学中的透射光干涉，实验上一般很难观察到[5]，原因是参与干涉的两束透射光光强相差较大，导致透射光干涉可见度小，这样就导致学生理解透射光干涉特征困难，为解决这个问题，笔者模拟了牛顿环透射光干涉[4]，如图3所示，仿真图像清晰、逼真，很容易看出透射光干涉条纹和反射光干涉条纹明暗互补，有助于学生全面了解牛顿环干涉。笔者模拟的牛顿环透射光干涉图片被教材《大学物理学》所选用[6]，该教材由华中科技大学出版社出版，面向全国发行。
3.3选择一些简单物理问题供学生进行数值计算与模拟作为课后作业
学习大学物理，给学生布置必要的作业是必不可少的，传统的作业基本就是要求学生做书本上的习题，这固然有助于学生理解和掌握物理学的基本概念、基本规律和基本方法，但对培养学生运用计算机解决问题的能力没什么帮助。为此，笔者在教学实践中，每学期都会布置几个简单的物理问题，要求学生运用计算机来进行数值计算和模拟。例如在讲完“真空中的静电场”这一章后，笔者要求学生模拟真空中一对点电荷的电场分布及等势线分布，学生通过数值计算与模拟软件，不仅圆满完成了任务，而且得到了锻炼，图4就是学生模拟的一对点电荷的电场分布及等势线分布，效果不错。
3.4指导学生参加数值计算与模拟方面的科技创新
为了能更进一步锻炼培养学生的数值计算与模拟能力，笔者一般会在所任教的班级中选择若干个学生组成若干个课外科技小组，参加全校每年举行的大学生科技创新活动，笔者要求这些科技小组以模拟物理问题作为课题研究方向，这样的课题研究往往是综合性的，需要学生全面理解和掌握数值计算与模拟技术。笔者一般在课后对学生进行指导，科技小组经过努力，完成了这方面的科技创新，有的科技小组带着科技作品参加科技创新答辩，得到了评委的肯定并获奖，图5、图6分别是笔者指导的届科技创新作品《基于matlab的光学仿真平台》主界面、届科技创新作品《电磁学仿真平台的设计与实现》主界面，这两项作品分别获得全校科技创新大赛三等奖和二等奖。
4结束语
数值模拟计算 篇3
关键词：爆炸与冲击 数值模拟 高精度格式 并行计算
High Presion and Large Scale Numerical Simulation of Explosion and Shock Problems
Ning Jianguo Ma Tianbao Wang Cheng
（Beijing Institute of Technology）
Abstract：Explosion and shock is a strong coupling problem， and it’s geometric， material and boundary conditions are both nonlinear. The problems often involve high speed， high pressure， high temperature， large deformation and interaction between a variety of material such as gas， liquid and solid. Numerical simulation for this problem not only need to accurately simulate the transient large deformation of various material and interface treatment between multi-material with high ratio of density under strong impact loading， also involves the modeling theory and high precision numerical method of materials with complex nonlinear mechanical behavior like dynamic constitutive and crushing， melting and vaporization under explosion and impact loading. In this paper， the conservation of the three dimensional finite difference WENO computational formats which always ensure its value positive was structured for the suspension of computing problems due to negative density or pressure in the high precision calculation of detonation wave propagation. This algorithm overcome the disadvantages of the traditional method which the negative value 0 or take the absolute value. In addition， the critical problem of the Euler parallel algorithm was studied for the three dimensional explosion and impact problems of large-scale parallel computing. Then the high performance simulation software of three-dimensional explosion and impact problem （EXPLOSION-3D） was designed based on MPI （Message Passing Interface）， and a practical plan of software testing was presented. The speedup， efficiency and scalability of the software were tested based on the numerical example of shaped charge jet， and the effect of the bottlenecks of software was discussed. At last， numerical results of typical three-dimensional explosion and impact problem show the capability of software in dealing with large scale complex engineering problems. The main innovation points of this work are to always ensure the value of high precision format positive， which is the Euler method computing problem. And then we propose parallelization strategies of Euler method， so as to realize the large-scale computing of explosion and impact problem and the computing scale can reach more than 100 million grids. The software can effectively solve the computing of large-scale complex explosion field which the commercial software cannot do.
Key Words：Explosion and shock； Numerical simulation； High precision scheme； Parallel computation
计算机数值模拟稳态磁场重联 篇4
计算机数值模拟方法在科研工作中有着广泛而重要的应用。磁场重联是等离子体能量转化和耦合过程的一种重要的机制,是空间等离子体中非常重要的物理现象和研究内容之一。关于磁场重联的数值模拟方法一般有三种:磁流体动力学(MHD)模拟方法、全粒子模拟方法、混合模拟方法。MHD研究模型和结果在重现磁层大体形状和动力学问题方面取得了成功,但是这种模拟中没有考虑到粒子运动的影响;既包括离子动力学也包括电子动力学的全粒子模拟强调了X线处的物理特性,但是对于研究出流区的大尺度(离子尺度)结构,这种模拟方法是具有数值局限性的。[1]于是一维和二维的混合模拟方法就发展起来了。混合模拟方法将离子视为粒子,运动满足牛顿方程,而作为无质量流体的电子满足动量方程。由于混合模拟主要研究与离子运动有关的低频电磁波,所以满足准中性条件:离子数密度等于电子数密度。
Petschek[2]模型是稳态磁重联的早期理论模型之一。在这个模型中,重联结构包括三部分:一个入流区,两个出流区及一个很小的中央扩散区。这个模型描绘的是电流片两侧具有相同的等离子体密度及等大反向磁场强度的对称情况。该模型认为X线两边分别出现一对慢激波,通过这些慢激波磁能快速的被转化为等离子体的热能和动能。这个模型可以应用到地球远磁尾的磁场重联研究中。
本文采用二维三分量(2.5维)[1,3]的混合模拟方法研究Petschek稳态磁重联模型下磁场重联的结构,研究发现磁场重联后形成了突出的隆起状磁场位形,在X中性点附近的一个很小的区域内,磁场分量By出现了四极分布特征。
二、数值模型
1. 基本方程
离子的运动方程为
电子的动量方程为
电子流的宏观速度由安培定律计算
式中,为电荷耦合常数。
最后,我们用法拉第电磁感应定律推进磁场
2. 初始条件
等离子体片的中心区是电流片,初始时刻位于z=0处,其法线方向沿着z轴方向。电流片两侧分布着反向平行的磁场,反平行的磁场分量沿x方向、导向场沿y方向。
磁场沿x方向分量的初始位形为
是初始电流片的半宽度。初始时电流片两侧具有相同的导向场,假定模拟区域的导向场初值处处为零,By0(28)0。为了快速的磁场重联,我们将z方向磁场分量B z设为零值。
离子数密度初始位形为
在模拟区域的中心(7)x,z(8)(28)(7)0,0(8)处放入非均匀常数电阻率用以触发和控制重联,其空间位形为
3. 边界条件
入流边界初始时无入流扰动;出流边界取自由边界条件;磁场以无旋、无散条件外推;粒子的边界条件通过在入流边界放置的缓冲区控制。
三、模拟结果
3.1重联结构
磁场重联过程中不同时刻的磁场位形如图1所示,图中的是重联层以外区域中离子的回旋频率,本模拟中时间用无量纲化;是重联层外侧区域中离子的惯性长度,本模拟中的所有长度都用其无量纲化。从图1中我们可以看出磁场位形随时间演化的情况。时,中心电流片两侧具有反向平行的磁场分量,在计算开始的初期,磁场位形随时间的演化比较缓慢。磁场重联开始发生后,在模拟区域中心部分的磁场拓扑结构发生了明显地改变,少量磁力线变得高度弯曲。随着模拟时间的增长,磁场位形改变地十分迅速,发生重联的磁力线的范围进一步沿电流片扩张,电流片附近磁场高度变形,重联区的大小也随时间的增长而变大。在计算的后期,磁场位形变化就呈现出无明显变化的现象,说明磁场重联的发展最后会进入稳态饱和阶段,磁场位形结构相对稳定。
如图2所示给出的是重联过程中不同时刻离子数相对密度的位形图。从图2中可以看出离子数密度分布随时间演化的情况。重联开始前,随着离开电流片距离的增加,磁场强度逐渐增强,等离子体密度逐渐下降。磁场重联过程中,X中性点附近的离子数相对密度与初始时的分布相比逐渐降低,但是磁场重联使得出流区中的等离子体密度变大了。这是因为等离子体在X中性点被磁张力加速,射向出流区,所以在X中性点出现了密度分布的极小值,出流区出现了密度分布的极大值。随着重联过程的进行,由磁重联而形成的新的磁力线逐渐往出流区外移动,这个过程中等离子体会随磁力线一起被带出重联区,因此在重联的演化过程中,高密度等离子体区的范围会逐渐向出流区边界推移。
3.2 By的四极分布
如图3所示的是X中性点附近,磁场的y分量By在xz平面内的等值线图,图中的实线表示的是正值,虚线表示的是负值。我们选取的区域是X中性点附近的一个很小的区域,该区域的尺度与离子的惯性长度尺度可比,但仍远大于电子的惯性长度尺度。这就导致在该尺度范围内电子仍可视为流体,但离子的粒子性特征就不能再被忽视,离子的运动会退耦脱离磁场线,在电磁场中离子和电子的运动情况不同,因而会产生Hall电流效应。考虑Hall效应后,即使初态导向场By0=0,t>0后仍会产生四极形的By结构,这是无碰撞磁场重联的一个基本特征[4]。因此从图3中我们看到了X线周围的磁场分量By出现了明显的四极形分布结构。这一模拟结果与全粒子模拟[5]结果基本一致。
四、结论
本文采用2.5维混合模拟的计算机数值模拟方法研究了Petschek稳态磁重联模型下X中性点周围关于电流片对称的矩形区域内磁场重联的过程,主要结论如下:
1. 磁场重联使得X中性点周围的磁场位形发生了明显地改变,磁力线高度弯曲,重联过程中形成了一种突出的隆起状磁场位形。等离子体在X中性点被磁张力加速,射向出流区,X中性点附近离子数相对密度减小,而出流区离子数相对密度增加。随着磁场重联的不断进行,高密度区逐渐向计算边界推移。
2. 在X中性点附近和离子惯性长度可比的一个很小的区域内,由于Hall电流效应的作用,即使初态导向场By0=0,重联发生后X线附近仍会产生四极形的By结构。
混合模拟这种计算机数值模拟方法虽然无法了解电子的动力学行为,但它可以部分地获得离子动力学的信息及其对重联发生处的等离子体的影响,并且与完全粒子模拟相比具有计算量要小得多的优势,因此在磁重联相关问题特别是离子动力学起主导作用的问题研究中具有广阔的前景。
参考文献
[1]王旭丹.磁尾磁场重联的二维混合模拟研究.大连理工大学硕士论文[D].2006.
[2]H.E.Petschek.Magnetic field annihilation[A].National Aeronautics and Space Administration,NASA Spec.Publ.SP-50[C],Washington:1964:425-439.
[3]Lin Y,Swift D W.A two-dimensional hybrid simulation of the magnetotail reconnection layer.Journal of Geophysical Research[J].1996,101(A9):19859-19870.
[4]周国成,等.地球近磁尾无碰撞磁重联事件.空间科学学报[J].2003,23(1):25-33.
数值模拟计算 篇5
Altix 450在三维定常N-S层流数值模拟计算中的并行计算效率研究
在共享内存服务器Altix450上,应用Fluent流体计算软件进行三维定常N-S层流数值模拟计算,研究了该系统的加速比,并行效率和可扩展性;在对网格数为500万和1 000万的两次计算对比中发现,当网格数据大的时候,该系统并行优势明显.结果表明,大型的`共享内存服务器在计算流体力学中对解决复杂流场数值模拟越来越重要,符合未来计算流体力学高性能计算的要求.
作 者：南江 高超 郑博睿 NAN Jiang GAO Chao ZHENG Bo-rui 作者单位：西北工业大学,翼型叶栅国防科技重点实验室,陕西,西安,710072刊 名：航空计算技术 ISTIC英文刊名：AERONAUTICAL COMPUTING TECHNIQUE年，卷(期)：40(2)分类号：V211.3关键词：高性能计算 Fluent 数值模拟 并行计算
高桩码头承载特性三维数值模拟 篇6
(上海海事大学 海洋科学与工程学院， 上海 201306)
0 引 言
高桩码头指由桩基及上部结构组成的各种码头，主要有梁板式、桁架式、前板桩、后板桩等结构形式.高桩码头结构轻、受力明确，适宜做成透空式；减弱波浪的效果好，适于软土地基；但结构单薄，耐久性差，构件易损坏并难于修复，对地面超载及装卸工艺变化适应性差.[1-6]
高桩码头处在海陆交接地带，环境脆弱，受力非常复杂.[7-10]目前，对一根单桩的受力和变形都无法进行精确的计算，对高桩码头就更是如此，所以，很有必要加强高桩码头承载性状的研究.目前对高桩码头的研究已取得较多的研究成果[11-14]，但将整个高桩码头作为一个整体，且考虑堆场、岸坡土体、桩基础、上部结构共同作用的研究较少.鉴于此，本文基于有限元软件ABAQUS开展高桩码头三维承载性状的研究.
1 数值仿真模型构建及工况设置
1.1 数值模型构建
ABAQUS是一套功能强大的工程模拟有限元软件，可以分析复杂的固体力学、结构力学系统，能够驾驭非常庞大复杂的问题和模拟高度非线性问题.ABAQUS拥有各种类型的材料模型库，可以模拟典型工程材料的性能，其中包括金属、橡胶、高分子材料、复合材料、钢筋混凝土、可压缩超弹性泡沫材料以及土壤和岩石等地质材料，含有非常丰富的岩土工程本构模型.ABAQUS经过在国内外近30年的广泛应用，其可靠性得到广泛而有效的验证.
根据《港口工程载荷规范》和《高桩码头设计与施工规范》并参考与高桩码头相关的文献建立ABAQUS的几何模型和计算模型.取纵向24 m为研究范围，包括：前承台3个排架，排架间距7 m；后承台7个排架，排架间距3.5 m.前承台为梁板柔性桩台结构，绝大部分为预制构件；前方起重机梁为矩形断面；后方起重机梁是局部加厚的钢筋混凝土面板，该面板为预应力板；横梁为非预应力梁；桩距3.6 m.对岸坡形式进行简化，坡度比为1∶3，简化为单层土进行模拟，接岸处设有混凝土挡土墙.前承台宽14 m，后承台宽15.7 m.
前承台横梁与纵梁均为矩形截面桩，其中：横梁尺寸为1.4 m×0.9 m；纵梁尺寸为0.9 m×0.5 m；前承台直桩、后承台直桩、叉桩均为50 cm×50 cm方桩，叉桩坡度比为3∶1，直桩长23.5 m，采用C40混凝土.
本模型中岸坡土体、桩、梁、板等均采用三维实体单元：桩采用C3D8R实体单元，将桩身划分为均匀的六面体单元；梁板及岸坡土体结构采用C3D4实体单元，将梁板及岸坡土体划分为四面体单元.数值模拟中，桩采用理想弹性模型，土采用理想弹塑性模型.数值计算模型见图1，材料参数见表1.
图1数值计算模型
表1 材料参数
1.2 数值模拟中工况设置
表2 工况设置
基于ABAQUS的三维有限元计算中，设置7种工况(见表2)，主要模拟承台载荷与堆场载荷及其不同组合情况下的桩基受力变形情况.
工况1考虑前、后承台及堆场均布载荷；工况2考虑前承台均布载荷；工况3考虑后承台均布载荷；工况4考虑堆场均布载荷；工况5考虑前、后承台均布载荷；工况6考虑前承台和堆场均布载荷；工况7考虑后承台和堆场均布载荷.
2 码头桩承载性状分析
通过数值分析，获得各工况下的各类位移和内力图.限于篇幅，仅选取工况1，4，6下的位移云图和矢量图，见图2.以码头中的一排桩为代表分析各桩桩身的位移和弯矩，各桩序号见图3.各工况下前、后承台各桩的桩身位移分布见图4和5.
从图2(a)和2(b)可以看出，由于在工况1下整个码头的总位移从堆场、后承台到前承台逐渐减小，后承台桩基础对岸坡土体的遮帘作用使得前承台桩基础位移减小；堆场下土体、承台板的位移方向主要为左下，即向下斜指海侧，而岸坡土体的位移主要为近水平指向海侧，使得岸坡中的桩的位移方向也主要为近水平指向海侧.
(a) 工况1下三维位移场云图(单位:m)
(c) 工况4下三维位移场云图(单位：m)
(e) 工况6下三维位移场云图(单位：m)
(b) 工况1下剖面中的位移矢量图
(d) 工况4下剖面中的位移矢量图
(f) 工况6下剖面中的位移矢量图
图3 中间排桩中的桩序号
图4在工况1～7下前承台桩身位移分布
图5在工况1～7下后承台桩身位移分布
由图2(c)～2(f)发现，工况4和6下码头的位移场云图和位移矢量图与工况1下的大体类似，但也有明显区别，主要体现在图2(c)～2(f)中后承台及其下的桩基础的位移都小于前承台，主要是由于高桩码头只受到堆场的均布载荷，土体位移呈现顺时针绕流的趋势，使得越靠近海测的岸坡其上翘趋势越明显，表现为前承台的位移大于后承台位移.因此，只有堆场载荷时前承台呈现明显的斜向上的移动趋势，对码头前沿港口机械的安装及稳定运行不利，在码头施工及运行时应给予注意.
通过对比工况4和6下的位移矢量图可以看出，土体的运动趋势大体相同，但在工况6下前承台受到向下的均布载荷因而靠近海测的岸坡的上翘趋势减缓，表现为前承台的位移主要集中于水平方向.在工况4下前承台位移是由于在堆场载荷作用下土体挤压使得前承台上翘，表现为前承台既有水平向海侧的位移又有向上的位移；而在工况6下由于前承台受到的竖向力作用抵消其向上的位移，在土的挤压作用下表现出水平向右的位移.
由图4和5可以看出桩身的位移主要有2种情形，在工况1，4，6，7下桩身整体位移较大，其余工况下桩身整体位移则较小.通过比较发现工况1，4，6，7都包含堆场载荷，这说明影响桩身位移的主要因素是堆场载荷.此外由于岸坡土的约束使桩基础7 m以下的入土段的位移较小且集中，上部的悬空段的位移较大且分散.在7个工况下，桩身向海侧的位移最大不超过20 mm.
由弯矩分布图(见图6和7)可以看出，桩身弯矩不会因为堆场载荷是否存在而分为两种不同的趋势，各个工况下的桩身弯矩分布情况较桩身位移集中，只是土体以上的桩身弯矩有所变化，变化范围主要为(-45～20) kN·m.另外，数值模拟结果显示，高桩码头应力最大的部位一般位于叉桩与桩帽的接触部位.
图6在工况1～7下前承台桩身弯矩曲线
图7在工况1～7下后承台桩身弯矩曲线
3 结 论
针对本文的地基情况，对各工况进行基于ABAQUS的数值模拟，得出结论：(1) 各工况下码头的位移各异，但有规律可循.总体来说，整个码头以向海侧的位移为主.(2) 对于有堆场载荷的工况，无论是前承台桩基还是后承台桩基，堆场载荷都是使其产生向海侧位移的主要因素，前、后承台上的竖向载荷对桩基础的影响相对较小；在各工况下桩身向海侧的位移最大不超过20 mm.(3) 各个工况下的桩身弯矩分布较桩身位移集中，桩入土段的弯矩很小，仅土体以上悬空段的桩身弯矩较大，最大值的绝对值在45 kN·m以内.(4) ABAQUS经过长期在岩土工程中的应用，已相当成熟和可靠.本文建立的三维高桩码头有限元模型可用来研究高桩码头的承载性状.
参考文献：
[1]王煜, 郑惠强, 李强. 基于知识工程的煤炭码头堆场分配策略[J]. 上海海事大学学报, 2012, 33(3): 26-30．
[2]陈华林, 李志胜, 郑书禧. 高桩码头结构修复加固[J]. 水运工程, 2012(9): 194-199.
[3]蒋建平, 刘春林, 蒋宏鸣. 遮帘式板桩码头三维地震动响应[J]. 上海海事大学学报, 2013, 34(1): 28-35.
[4]牛兴伟, 沈才华, 孙会想. 高桩码头分缝处工后缝宽有限元分析[J]. 水运工程, 2012(6): 121-124.
[5]顾天意, 梁承姬. 基于矩阵式遗传算法的集装箱码头堆场空间资源分配优化策略[J]. 上海海事大学学报, 2012, 33(2): 40-46.
[6]刘晓平, 凌威, 林积大, 等. 底梁式全直桩高桩码头水平力作用下的特性分析[J]. 中国港湾建设, 2012(2): 13-14.
[7]郑剑, 肖英杰, 白响恩, 等. 基于缆桩表面应力的船舶缆绳载荷测量方法[J]. 上海海事大学学报， 2013, 34(1): 1-4.
[8]王晋, 张华平. 高桩码头排架计算模型简化条件对横梁弯矩的影响[J]. 港工技术, 2012, 49(2): 27-29.
[9]李申, 钟小帅. 高桩码头异常靠泊有限元仿真模拟[J]. 水运工程, 2012(1): 64-68.
[10]廖雄华, 张克绪, 王占生. 岸坡开挖扰动对天津港高桩码头结构安全性影响的数值分析[J]. 中国港湾建设, 2002(2): 33-38.
[11]武清玺, 张旭明. 桩基码头空间结构的内力计算与分析[J]. 水利水电科技进展, 2004, 24(3): 21-23.
[12]王年香, 魏汝龙. 岸坡上桩基码头设计方案的比较分析[J]. 水利水运科学研究, 1995(1): 43-54.
[13]廖雄华, 张克绪. 天津港高桩码头桩基-岸坡土体相互作用的数值分析[J]. 水利学报, 2002(4): 81-87.
数值模拟计算 篇7
在传统的阴极保护工程中,大多采用实际测量或经验估计的方法来获得被保护体表面的电位分布。然而,某些被保护结构如海底管道、海洋平台、深埋的钢桩等实地测量难度很大,对于新项目更不可能事先在现场测试,而经验公式又无法对复杂结构表面的电位分布进行准确的预测,传统方法已难以满足越来越高的安全可靠性以及经济性的要求。应用数值模拟技术代替传统的经验设计模式,利用数值方法求解阴极保护体系的电位和电流分布问题已成为近年来阴极保护领域的热点,并且已经在地下长输管道、近海石油平台等中得到较好的应用,实现了优化设计[1]。
在一般阴极保护体系中,影响阳极发出电流的因素包括阳极接水电阻、阳极极化电阻与阴极极化电阻。但目前在实际阴极保护工程设计中,却主要考虑阳极接水电阻,并将阳极极化电阻(极化电位-1.05 V)和阴极极化电阻(极化电位-0.80 V或-0.85 V)视作定值,而不考虑阳极极化性能的具体影响[2]。对于实际的阴极保护体系,由于阳极极化性能的影响,阴极保护系统的状态并不是一成不变的,而是随着时间发生变化。因此在实际阴极保护工程中,如果阳极质量不佳,其极化性过大而限制了发出电流,往往会使得阴极保护效果不理想,所以在阴极保护工程设计中也应当考虑阳极极化电阻因素[3~5]。
本工作对不同类型阳极的阴极保护达到稳态时的平均电位和电流密度进行了数值模拟,在考虑和不考虑阳极极化电阻两种情况下,通过对比、分析二者的差别,考察阳极极化电阻对阴极保护效果的影响及不同类型阳极的阴极保护效果的优劣,最终判断出阳极极化对阴极保护数值模拟计算的影响。
1 物理模型
1.1 描述方程
国内外对阴极保护体系电位分布的数值研究中,大多采用拉普拉斯方程作为电位分布描述方程,方程形式为k▽2φ=0。式中,φ为电位,k为所研究区域内介质的电导率。要保证拉普拉斯方程的有效性,需要假设所研究区域内导电介质均匀,系统内没有电流的得失,不存在源点或汇点,系统状态不随时间发生变化即处于稳态。当所考察区域内有场源存在时,电位分布方程的形式变为泊松方程:k▽2φ=qs。对于实际的阴极保护体系,由于阴极垢层的形成以及环境参数变化的影响,阴极保护系统的状态随着时间发生变化。本工作采用动态极化曲线来模拟阴极表面的特性,在足够小的时间段内假定阴极保护系统处于稳态,应用拉普拉斯方程进行求解,得到阴极表面电位和电流分布情况。
1.2 试样制备
实施阴极保护的空间区域为1 m×1 m×1 m的正方体;此正方体6个面中,仅底部1个面为金属面(即导电面),其余5个面均为绝缘面;腐蚀环境为海水。
阳极共3种类型:圆柱型阳极直径0.01 m、长0.04m,方条型和长直翼型阳极的长和截面积均与圆柱型阳极相同,因此3种阳极体积也相同。阳极位于正方体区域底面对角线处,呈中心对称放置,即阳极轴线平行于底面一边,阳极中心在底面投影与对角线交点重合;阳极底面距正方体底面为0.01 m。有关物理参数取值如下:阳极开路电位为-1.050 V(vs SCE,下同);海水电阻率为0.25Ω·m;阴极自腐蚀电位为-0.650 V。
1.3 测量系统
室内模拟阴极保护试验自动测量系统见图1。该测量系统主要由钢板、牺牲阳极试样、Ag/Ag X参比电极、精密取样电阻、采样器中继接线板和数据采集/存储器等构成。图中同时标出了每块钢板表面安装的1个牺牲阳极试样A和8个Ag/Ag X参比电极的布设位置及其布线方式,a和b分别为阳极试样和钢板上同一位置的两根引线焊接点,其中a是阳极试样接点,b是阴极钢板接点,即地线。采样器中继接线板上共有12个接线端子,1~8对应于8个Ag/Ag X参比电极,精密取样电阻两端接入接线板输入端9和10,再由其输出端接至采样器上的接线端子。该端电压经过差分放大10倍后作为电压信号存储,并由此计算阳极的发出电流值。
2 结果与讨论
以阴极表观面电阻率分别取0.1,1.0,10.0,100.0Ω·m24种情况为例,分别在不考虑阳极极化电阻和考虑阳极极化电阻为0.002 388 9Ω·m22种情况下,对不同类型阳极在阴极保护达到稳态时的平均电位和电流密度进行数值模拟。
不考虑和考虑阳极极化电阻情况下阴极保护数值模拟见图2。对比是否考虑阳极极化电阻情况下阴极保护达到稳态时的数值模拟结果可以看到:考虑阳极极化电阻情况时数值模拟给出的阴极保护效果要大幅低于不考虑阳极极化电阻情况同等条件下的阴极保护效果。由于阳极极化电阻限制了电流的发出,阻碍了阴极保护效果。因此,是否考虑阳极极化电阻对阴极保护效果包括稳态的平均电位和电流密度的影响关系较大。
阴极保护数值模拟中,物理模型与上文完全相同;使用相同体积不同表面积的牺牲阳极,在阳极极化电阻取不同值的情况下进行数值模拟,对阴极表观面电阻率以其时-空变化方程式为边界条件进行处理。分别取阳极极化面电阻率为0.001,0.003,0.005,0.010Ω·m2,进行不同阳极极化电阻率情况下阴极保数值模拟,利用以上不同阳极极化电阻情况下阴极保护数值模拟结果,作3种阳极的极化电阻对阳极发出电流、阴极表观面电阻率和阴极平均电位的相关性分析,见图3。
从图3可以看到,阳极发出电流、阴极表观面电阻率和阴极平均电位对阳极极化电阻具有明显的相关性,随着阳极极化电阻率的增加其阳极发出电流、阴极表观面电阻率和阴极平均电位均相应下降,并且在同等条件下牺牲阳极接触到的表面积越小,对应值越低。因此根据图3可得到其他阳极极化电阻数值时阳极发出电流、阴极表观面电阻率和阴极平均电位的对应数值。同等体积的3种阳极阴极保护时的效果有明显差异:圆柱型阳极保护效果最差,方条型阳极保护效果居中,长直翼型阳极保护效果最优。
3 结论
(1)在数值模拟计算中,是否考虑阳极极化电阻对包括稳态的平均电位和电流密度在内的阴极保护效果影响比较大。
(2)在实际阴极保护工程中,如果阳极质量不佳,其极化性过大会使得阴极保护效果不理想。因此在阴极保护工程设计中应当考虑阳极极化电阻因素。
(3)同等体积的3种阳极阴极保护时的效果有明显差异:圆柱型阳极保护效果最差,方条型阳极保护效果居中,长直翼型阳极保护效果最优,可为海洋工程中合理选用牺牲阳极提供参考。
参考文献
[1]郝宏娜,李自力,王太源,等.阴极保护数值模拟计算边界条件的确定[J].油气储运,2011,30(7):504~507.
[2]侯德龙,宋月清.铝基牺牲阳极材料的研究与开发[J].稀有金属,2009,33(1):96~100.
[3]陈绍伟.钢在海水中的阴极极化特点及阴极保护系统设计中的新途径[J].腐蚀科学与防护技术,1996,8(1):1 725~1 729.
[4]邱枫,徐乃欣.码头钢管桩阴极保护时的电位分布[J].中国腐蚀与防护学报,1997,17(1):12~18.
喷管内部流场的数值模拟计算与研究 篇8
喷管广泛应用于汽轮机、风洞、引射器等动力和实验装置中。喷管的研究手段有两种:试验和理论分析。前者对人力、物力、财力和时间的消耗巨大, 后者由于计算工具的限制, 使得数值解在工程实际中完全失去了意义。
近年来得到迅速发展的CFD类软件在实现对动力装置内部流场的数值模拟计算方面, 有着极其强大的功能。CFD类软件可以用作在流动基本方程控制下对流体流动的数值模拟, 通过这种数值模拟, 可以得到极其复杂问题的流场内速度、压力、温度、马赫数等的分布, 此外, 与CAD联合, 还可进行结构优化设计等。FLUENT是一个通用的CFD软件包, 用来模拟从不可压缩到高度可压缩范围内的复杂流动。由于采用了多种求解方法和多重网格加速收敛技术, FLUENT能达到最佳的收敛速度和求解精度。
2 喷管内空气流动的控制方程
Navier-Stokes (N-S) 方程组是迄今为止描述流体流动最为完备的控制方程组, 适用于本研究对象的直角坐标系下非定常可压缩的二维轴对称流动的守恒型雷诺平均N-S方程为:
在上述方程中, ρ为密度, vx, v分别为x, y方向的速度矢量, p为压强, μ为粘性系数。
3 FLUENT环境下的数值模拟计算
3.1 喷管模型
本文研究对象为一轴对称Laval喷管, 总长35mm。收敛段长5mm, 喉部直径4mm, 壁型为一半径为5mm的1/4圆弧, 扩张段壁型为一斜率为1:30的斜线。喷管进口温度为295K, 总压为88750Pa。湍流模型选择Spalart-Allmaras湍流模型, 采用耦合隐式求解器。
采用三角形网格, 并对喷管的喉口和边界进行加密。对流动项的离散过程中采用具有二阶精度的二阶迎风格式, 扩散项的离散采用中心差分格式, 对离散方程的求解采用耦合隐式求解方法。边界条件分别为:入口选为压力入口;出口为压力出口;边界采用壁面边界。残差值均选为0.001。图1为Laval喷管的网格分布图。
3.2不同工况下的数值模拟计算结果
模拟中选取了四种工况, 即喷管出口截面下游环境真空度为分别为0.02、0.03、0.05、0.07Mpa。
3.2.1喷管出口质量流量曲线 (见图2)
3.2.2喷管轴线处马赫数分布图 (见图3)
3.2.3喷管壁面压力分布图 (见图4)
3.2.4局部放大速度矢量图 (见图5)
4 结果分析
由图2看出, 随着出口真空度增加, 流量开始增加, 当喷管的喉部压力达到临界压力时, 喷管中的流量达到最大值。再降低背压, 流量m基本保持不变;由图3看出, 当喷管出口截面环境真空度为0.02~0.05Mpa时, 在喉部以后的扩张段的某一位置, 喷管轴线处的马赫数由大于1突然降到小于1, 即在此处产生正激波。用与喷管入口的距离来表示激波产生的位置, 其中喷管入口在试验中的位置为205mm, 可得激波产生的位置见表1。
从激波位置表中可以看到, 喷管中激波的位置随着背压的增加 (真空度减小) 逐渐向喷管喉部移动;
由图4可以看出, 当喷管出口截面真空度为0.07Mp时, 喷管壁面的压力始终下降, 没有出现边界层分离;而喷管出口截面真空度为0.02~0.05Mp时, 在喷管喉口以后出现压力由下降趋势突然变为增加, 就认为此处即为边界层分离位置, 采集分离位置的数据, 得到边界层分离位置数据表。
由表2可以看出, 随着背压的增加 (真空度减小) , 分离位置逐渐向喷管喉口部位移动。从分离处压力一系列数据我们还可以看出, 随着背压的增加 (真空度减小) , 边界层分离处压力也是随之增加的;
此外, 从图5可以看出, 当喷管出口截面真空度为0.03~0.05Mpa时, 喷管内出现局部回流, 且随着背压的增加, 回流的位置越靠近喉部;当喷管出口截面真空度为0.07Mpa时, 在喷管出口截面近壁面处出现局部回流, 当喷管出口截面真空度为0.02Mp时, 无回流产生。
5 结论
5.1 背压与正激波位置的关系表明, 喷管中产生激波的截面随着背压的降低 (真空度增加) 逐渐向喷管出口方向迁移, 因此在设计、制造喷管这一类部件时, 必须充分考虑喷管在不同背压下通道内气体流动特性。
5.2 背压与边界层分离位置的关系表明, 喷管中边界层分离的位置随着背压的降低 (真空度增加) 逐渐向喷管出口分向迁移。为减少边界层分离对喷管流场的影响, 应减小背压, 使边界层分离位置向喷管出口方向移动。
参考文献
[1]韩占忠等, Fluent流体工程仿真计算实例与应用[M].北京:北京理工大学出版社, 2004.
[2]王福军, 计算流体动力学分析[M].北京:清华大学出版社, 2004.
[3]张廷芳, 计算流体力学[M].大连:大连理工大学出版社, 1992.
[4]Laval喷管内激波/湍流边界层干扰的数值模拟[D].鞍山:辽宁科技大学, 2007, 3.
[5]魏庆鼎.湍流边界层分离的实验研究[J].空气动力学学报, 1985, 3:30-38.
数值模拟计算 篇9
燃烧分为层流燃烧和湍流燃烧。层流燃烧在理论中较多提及,实际遇到的燃烧过程几乎都是湍流燃烧。湍流自身空间上的尺度多重性和时间上的高频脉动性以及燃烧室流场中湍流与燃烧两种物理现象的强烈耦合,使得燃烧流动现象异常复杂,并伴随有剧烈的放热化学反应,而且湍流流动的工作机理至今也尚未被充分掌握,这在对湍流燃烧进行数值模拟研究时是不得不面对的难题[1]。目前,随着计算流体力学(CFD)的不断发展,采用CFD方法对湍流燃烧流动进行数值模拟分析,由于成本低、周期短和对多尺度(包括空间和时间尺度)问题实用性好且非常适合对燃烧问题进行研究以及燃烧装置的设计和改造,因此,开展数值模拟的研究工作不仅对于探讨燃烧流动反应机理有重要意义,而且也是一种对试验结果作对比分析的重要手段。
首先对湍流燃烧数值模拟系统体系结构和系统并行计算原理做了阐述,然后在计算流体力学平台OpenFOAM上,基于PaSR模型建立了燃烧流动的数值模拟系统,并对煤气和空气湍流射流扩散燃烧反应进行了数值模拟,实验结果表明系统能够较好地展示湍流燃烧流动反应特性。
1 系统体系结构
图1是湍流燃烧数值模拟系统体系结构图。从图中可以看出,以OpenFOAM作为计算流体力学平台,针对湍流燃烧反应模型设计出湍流燃烧流动解算器,同时使用网格生成软件得到燃烧室的总体网格结构并将其导入到解算器目录中,然后使用网格分区算法将计算网格剖分成若干区域并分配给对应计算节点,最后在MPI[2]并行计算平台使用解算器对燃烧流动反应进行数值求解。
2 并行计算原理
为得到燃烧室燃烧流场细节,需要大量网格,同时化学反应涉及到的反应组分和反应式比完全气体的计算量要多出好几倍,然而这种计算规模用串行计算显然无法满足要求,只有通过并行计算才能充分满足研究需要。
并行求解采用消息传递的编程模式。消息传递并行程序有函数分解和域分解两种形式[3]。函数分解形式也被称为功能分解,即将一个大问题分解成若干个子问题,然后对其分别求解;域分解形式是将一个大问题区域分解成若干个较小问题区域,然后对其分别进行并行求解。燃烧流场数值模拟采用域分解形式进行并行求解,即使用网格分区算法将计算区域划分成若干小区域,然后使用消息传递库进行网格点间的信息交互。网格并行分区算法的关键是区域分解方法,区域分解是指把整个流动区域分割成N个子区域并分配给N个处理器计算,把子区域的初始流场信息、几何信息分别装载入各个子区域对应的处理器内存中,在每一个处理器中启动计算进程来完成流场计算。计算过程中各处理器通过通信完成边界数据交换,同时收集全场数据完成收敛准则判别,并按需要进行保存等其他操作[4]。图2为网格分区划分示意图。
基于消息传递并行编程模型的并行编程语言主要有MPI和PWM两种。MPI吸收了现存许多系统的最突出优点,是当今最为流行的用于并行编程的消息传递库标准, MPI有很多实现,并行编程采用MPI的一个开源实现OpenMPI[6]。
3 平台搭建
3.1 OpenFOAM平台
解算器建立在近几年流行起来的计算流体力学软件OpenFOAM的基础之上。OpenFOAM的全称是Open Field Operation and Manipulation,它的核心是使用C++编写的面向对象的计算流体力学类库,有较高效率的偏微分方程求解模块。OpenFOAM是一款开源软件,这就可以方便用户编写解算器。从程序实现功能角度来看,OpenFOAM软件包括前处理、核心求解器和后处理三大模块[7]。前处理主要用网格设计,然后使用求解器对问题求解,后处理是将计算结果转化成可视化图像以进行结果分析,使用者可根据具体求解情况来设计求解器。
3.2 燃烧模型
本文在进行数值模拟过程中使用的燃烧模型为k-ε双方程模型。该模型对于一般流动情况,如平壁边界层、无力平面射流、管流等流动能给出相当满意的计算结果,而且计算量小。
k-ε双方程模型[8]如下所示
在上述方程中,Gk表示由于平均速度梯度引起的湍动能产生,Gb是用于浮力影响引起的湍动能产生;YM是可压速湍流脉动膨胀对总耗散率的影响。湍流黏性系数
3.3 化学反应动力学机理
在化学反应动力学因素起主导作用的燃烧中,化学反应动力学机理起着决定性作用[9]。本文使用的化学反应机理是OpenFOAM耦合Chemkin[10]后的GRI3.0详细化学反应机理数据库,共涉及煤气和空气燃烧的53组分325个基元反应,其中反应动力学参数可参看文献[11]。其中,对于燃料CmHn,其和氧气的总包反应为
4CmHn+(2m+n)O24mCO+2nH2O (3)
对于CO,其和氧气的总包反应为:
2CO+O22CO2 (4)
3.4 数值计算方法
计算描述燃烧过程的多组分流动守恒方程时,采用基于SIMPLE 方法[12,13]的全隐式有限体积法进行求解。在模拟计算时,对扩散项、对流项、时间项的差分格式分别采用了二阶线性修正的高斯差分格式、二阶迎风、一阶隐式的欧拉差分格式。SIMPLE算法的求解流程示意图如图3所示。
3.5 燃烧反应求解器
构建湍流燃烧反应求解器是根据SIMPLE解算器全隐式有限体积法构建的。计算器包括Make目录、头文件和主程序三部分。程序的主程序为reactingFoam.C,Make目录包括options文件和files文件,Files文件用于指定当前要编译的文件,这里不包含头文件,都是*.C文件。Options文件定义编译选项,用于指定编译用到的头文件位置及其动态库。头文件主要是主程序预编译用到的程序,主要包括hEqn.H, pEqn.H,UEqn.H和YEqn.H等文件。
4 计算结果及分析
4.1 初始条件和边界条件
在OpenFOAM平台上对特定燃烧室网格环境下的煤气和空气湍流射流扩散燃烧反应进行数值模拟。求解的计算区域为轴对称问题,计算区域如图4所示,这是钝体的一部分,即从钝体轴向划分5°的区域,钝体半径为50 mm,燃料入口半径为3.6 mm。钝体网格结构简图如图5所示,全部采用结构化网格离散计算区域,网格数为5 580,以便在OpenFOAM中模拟轴对称问题。反应初始条件如表1所示。
4.2 计算过程及结果分析
在数值模拟阶段,将反应物和伴随流分别从fuelInlet和airInlet以118 m/s和40 m/s的速度注入燃烧室。反应物注入后在燃烧室内进行燃烧反应,经过21 600 s的计算,燃烧反应稳定。
图6、图7和图8分别为燃烧反应稳定后产物二氧化碳、水及速度分布图。从图中可以看出:燃烧反应从喷口到狂涨段都富含燃烧反应产物二氧化碳和水,并在输出口富集。分布云图较好地反映了反应流场在反应稳定时的一些基本特征。
作为燃烧反应的参数之一,温度是与实验结果对比的重要参数。图9是燃烧反应稳定时燃烧室截面中取一条平行于X轴的端点坐标为(0,0)和(0,0.4)的直线,并根据此直线上坐标对应的温度得到的温度分布曲线图,从曲线图中可以得出最高温度超过2 200 K。
Sydney钝体驻定火焰是国际湍流扩散火焰测量与计算研讨会(TNF workshop)的标准系列火焰之一[14,15],常用于与仿真结果做对比。实验温度分布与Sydney钝体驻定火焰的HM1实验结果基本吻合,尤其是对燃烧得到的最高温度所对应的求解结果效果较好。
从图6、图7和图8可以看出,火焰上游没有出现由燃料射流扩散而形成的回流区,这使得模拟出的效果偏向层流燃烧,因此在数值模拟时,k-ε双方程模型中k和ε的设置还需要进一步探讨,以便得到更好的模拟效果。另外,由于在数值模拟时考虑到涉及详细化学反应机理而且求解的规模巨大,为了提高计算效率,将网格规模设计较小,这使得有些区域网格不够密集,尤其是火焰上游区域,这对数值模拟精度也有一定的影响。
4.3 并行效果分析
搭建的并行计算平台为西南科技大学网络与高性能计算实验室IBM刀片服务器。服务器处理器为Intel Xeon(R) E5405 2.00GHZ,内存为8G。图10给出了在并行计算平台上并行计算消耗的时间与参与计算的节点数之间的关系。从图中不难看出,随着参与计算的节点数的增加,计算效率明显提高,几乎呈现倍数减少的关系。
5 结论
本文针对湍流燃烧流动模型设计了解算器并基于并行计算和OpenFOAM平台成功搭建了湍流燃烧流动数值模拟系统,同时又构建了基于PaSR模型的求解器。
通过使用该求解器在特定的燃烧室网格环境下对煤气及空气燃烧流动进行数值模拟并采用GRI3.0详细化学反应机理发现:模拟结果表明系统能够较好地展示湍流燃烧流动的反应特性,符合工程研究需要;另外,使用并行计算技术后,计算速度明显加快,这对于开展大规模高网格数的湍流燃烧流动数值模拟工作有较好的促进作用。
数值模拟计算 篇10
水下系泊监测平台是水下监测系统图像采集仪器与设备的支撑平台,横向布放于安检水域,安装在其上的声纳线阵垂直水面向上,连续发射探测脉冲,实现声纳线阵波束对安检区域的扫描。这类系统需要长期工作于水下环境,能够连续、实时地测量周围水域的多种环境参数,为了获得实时稳定的信号,研究水下系泊监测平台在复杂外载荷作用下的动态响应具有十分重要的意义。
水下系泊监测平台是一种新型的监测装置,承受了波浪力、水流力等很多不确定性载荷,工作环境复杂,其设计尚无现成的设计规范和资料可以参考。这些载荷的随机性都必须在结构设计和分析中予以考虑。国内外对于这种新型监测平台的研究文献很少,但是这种监测平台在力学特性上和海洋系泊平台及水下悬浮隧道有一些相似之处,因此在研究中可以借鉴国内外学者对海洋系泊平台和水下悬浮隧道的研究成果。Chen[1]进行了深海张力腿平台的动力学响应及模糊控制方法研究,Kuiper等[2]研究了自由悬挂的水下流体输送管道的动力稳定性及波浪载荷下的动态响应。李晓平[3]应用Huston提出的缆索有限段方法,考虑流场作用力性质,建立了水下缆索三维有限段模型的动力学方程,但没有考虑缆索受到突变载荷以及浮体和缆索的耦合效应引起的振动。麦继婷[4]研究了波流作用下悬浮隧道的动力学响应,把悬浮隧道的力学模型简化为两端简支的梁,这种模型和水下系泊监测平台的空间受力状态不符合。因此本文在借鉴上述学者的研究成果的基础之上,考虑系泊监测平台浮体细长比小、两端处于无约束为自由状态的特征,研究了监测平台受到冲击的力学模型,借鉴Hopkinson的冲击载荷理论,建立了水下系泊监测平台计算模型和数值模拟模型,分别采用解析求解和数值模拟的方法,利用水力学理论和海洋工程结构力学理论对其进行求解,并通过实例给出了计算结果,分析了水下系泊监测平台的动态响应。
1 解析求解
水下系泊监测平台的组成如图1所示,由锚、系泊缆、被系泊的浮体以及安装在浮体上的声纳、传感器组成。
由文献[4]可知,水下浮体结构的动态特性受到所处水域的流速、水深、浮体放置深度和截面形式等因素的影响,垂直浮体轴向的流体力和绕流阻力最大,沿浮体轴向的流体力和绕流阻力最小,根据监测平台的受力特点和与系泊缆的相互作用可知,影响系泊缆张力的主要运动是监测平台在流速方向的平面运动。图2为水下系泊监测平台的简化示意图。
在水文环境的改变或外来干扰的激励下,浮体会产生不稳定的运动,破坏其初始平衡。由于系泊缆较轻,在外来干扰的作用下,会产生突然的绷紧或松弛,引起浮体和系泊缆之间的碰撞,产生冲击载荷,冲击力可以根据Hopkinson的冲击载荷理论进行计算[5]。同时,浮体还受到系泊缆的拉力、水流方向的拖曳力、流体的阻力、江水的浮力、自身重力的作用。设坐标原点O位于浮体中心的初始位置,忽略系泊缆上的波浪力,根据牛顿力学定律,建立浮体与系泊缆碰撞时的运动方程:
式中,y为垂直于来流方向的位移;z为来流方向的位移;FY为流体垂直升力FL与Y向惯性力FIY的合力;FZ为浮体受到流体的拖曳力Fd与Z向惯性力FIZ的合力;m为浮体的质量;Ac为系泊缆的横截面积;ρc为系泊缆的密度;E为系泊缆的弹性模量;k为系泊缆的刚度系数;C为水动力的阻尼系数。
浮体与系泊缆碰撞结束之后将产生初始位移和初始速度,并继续运动,直到达到平衡位置为止,根据牛顿力学定律,建立浮体与系泊缆碰撞后的运动方程:
1.1系泊缆静张力及静变形
浮体处于静平衡位置时的受力如图3所示。
浮体在静平衡位置受到系泊缆的张力T1、T2、浮力F和浮体重力G的作用。由于系泊缆的安装位置对称,所以有
式中,θ为系泊缆的安装角度。
根据材料力学原理,系泊缆在张力T1、T2作用下的静变形量Δl=T1l/(EAc),其中,l为系泊缆的长度。
1.2浮体运动时系泊缆变形量及张力
1.2.1 系泊缆变形量
式(1)～式(4)右边的第二项为考虑弹性变形时系泊缆所产生的张力,其大小为系泊缆刚度和变形量的乘积。水下监测平台在水流的作用下,两侧系泊缆的变形是不同的,设由于浮体的运动产生的系泊缆的弹性变形为ui,参考图2,并根据位移变形协调条件可知
1.2.2 系泊缆刚度
系泊缆是水下系泊监测平台的主要约束构件,对整个平台的动态性能有着重大的影响。系泊缆作为只受单向拉应力的柔性构件,在其自身重力和水流力的作用下将产生一定垂度。在计算求解的过程中,应考虑垂度的影响,目前对系泊缆的垂度效应的计算主要采用等效弹性模量法。水下监测平台系泊缆处于绷紧状态,在分析过程中,常将其模拟成直杆桁架单元,认为其拉力沿弦向作用,但是由于系泊缆有一定的垂度,其拉力和弦向作用力之间存在一定的差异,根据Ernst建议,采用等效弹性模量法来修正它们之间的偏差[6],具体的计算公式为
式中,q为单位长度系泊缆的质量;lx为系泊缆水平投影长度;T为系泊缆的初始张力。
初始张力可以通过浮体的静平衡位置求解。因此系泊缆的刚度系数可以根据材料力学理论求解:
1.2.3 系泊缆运动过程中的张力
考虑式(6)、式(7),可得系泊缆在Y方向弹性变形产生的力:
同理可得系泊缆在Z方向弹性变形产生的力:
1.3水流作用力
Fd、FIY、FIZ可根据Morison公式计算[7],考虑浮体和水流的相对运动,设江水为均匀定常流体,则有
式中,CD为黏性阻力系数;ρ为江水密度;D为浮体直径;v为浮体布放位置江水的流速;CM为惯性力系数;L为浮体长度。
FL可根据“卡门涡旋”理论求解[7]:
式中,CL为升力系数。
1.4水动力阻尼系数
根据文献[8],水动力阻尼系数为
式中,ω为浮体振动的无阻尼固有频率;AL为浮体的横截面积;U为系泊缆稳态响应的幅值。
1.5监测平台与系泊缆碰撞后的运动方程
考虑到系泊缆是3对,把式(5)～式(16)代入式(1)～式(4)可得浮体与系泊缆碰撞时的运动方程:
浮体与系泊缆碰撞后的运动方程为
式中,H为系泊监测平台所处水域总水深;h为水面到浮体中心位置处深度。
2 长江断面垂线流速分布特性
水下监测平台的布放位置为湖北宜昌的黄陵庙水域,根据该水文站的统计资料分析,该水域的水文断面垂线流速分布符合卡拉乌舍夫椭圆流速分布规律[9]:
式中,v0为水面流速;p为流速分布参数,p=0.6;η=h/H。
由式(21)按积分法计算垂线平均流速为
参考图2,浮体位置处江水的流速为
在本项目的研究中,对黄陵庙水域的水文环境进行了测试,测试结果见表1。
从表1可以看出,测点3的水流量最大,水深最大,平均流速为2.27m/s,选取测点3为水下监测平台布放位置进行分析。
3 数值模拟
3.1模型的建立
为了验证解析解的计算结果,采用Msc/Nastran的流固耦合分析方法进行分析。由于浮体是由钢管焊接而成的圆柱形薄壳体,建模时按其受力特性采用2Dshell板壳单元;系泊缆采用等效弹性模量法,以Rod单元模拟,建立如图4所示的流固耦合有限元分析模型,局部放大如图5所示。流体的附加作用以虚拟质量的方法作用到模型上[10]。
3.2流固耦合面的定义
定义流固耦合面的目的是为了让欧拉网格中定义的材料(水)与拉格朗日网格定义的材料(浮体)发生相互作用。如果不定义耦合关系,即使浮体完全处于水域的网格之中,也不会对水的流动产生任何影响,同时自身也不会受到流体力的作用。要在浮体与水之间建立耦合联系,首先要在浮体上定义一层耦合面,该面是浮体与水之间相互作用力的传递者。对于水,该面充当流场边界,同时,流场的作用力使得外力作用在耦合面上,引起浮体单元的变形,浮筒与水的耦合面的定义如图6所示。
在进行流固耦合分析时,耦合面应当是封闭的,而且必须具有正体积,这就要求所有耦合面单元的法线方向都指向单元外侧,如图7所示。而且封闭的耦合面要位于欧拉网格内,否则,耦合作用不会发生。
4 算例
4.1设计参数
水下系泊监测平台的设计参数如表2所示。
4.2浮体响应量的解析解
把数据代入式(17)、式(18),求解可得碰撞过程中浮体Y、Z方向的位移变化曲线,如图8、图9所示。
从图8、图9可以看出,浮体在冲击结束之后在Y方向产生0.29m的位移,在Z方向产生0.33m的位移。
以碰撞结束之后的位移为新的边界条件求解式(19)、式(20),可以得到浮体在冲击运动结束直到浮体稳定时的响应特性。
从图10、图11可以看出,浮体在冲击结束后在流体的拖曳力、流体的垂直升力、浮力的作用下继续运动,最终在水的阻尼力作用下达到平衡位置,但是在来流方向浮体很快趋于平稳,在垂直于来流方向受水流涡激升力的影响浮体会产生振动,此时Y方向的运动位移达到1.14m,Z方向的位移达到0.43m。
4.3浮体响应量的数值解
求解流固耦合分析模型,可以得到浮体在Y、Z方向的位移的响应量,分别如图12、图13所示。将计算结果进行对比,如表3所示。
从表3可以看出,水下监测平台的Y方向位移响应值的解析解和数值解误差为2.6%,Z方向响应值的误差为7.2%。产生误差的原因是由于解析求解过程中浮体是作为刚性体进行计算的,而数值解的求解过程中以柔性体来进行计算,且采用了浮体的有限元模型,因此流固耦合的数值解的结果更为精确、可靠。
5 结论
(1)水下系泊监测平台浮体在流体的拖曳力、流体的垂直升力、浮力的作用下,会产生较大的垂荡位移和横荡位移,但是在来流方向的横荡位移很快趋于平稳,垂直于来流方向的垂荡位移受水流涡激升力的影响会产生振动。
(2)流固耦合的计算模型和数值分析模型的计算结果一致性较好,但也存在一定误差,产生误差的原因是由于解析求解过程中浮体是作为刚性体进行计算的;而数值解的求解过程采用有限元模型,浮体是以柔性体来进行计算,因此流固耦合的数值解的结果更为精确、可靠。
(3)在进一步研究水下系泊监测平台浮体稳定性的过程中,涡激振动产生的影响是不能忽略的。
参考文献
[1]Chen Cheng-Wu.The Stability of an Oceanic Struc-ture with T-S Fuzzy Models[J].Mathematics andComputers in Simulation,2009,80(2):402-426.
[2]Kuiper G L,Metrikine A V.Dynamic Stability of aSubmerged Free-hanging Riser Conveying[J].Fluid Journal of Sound and Vibration,2005,280(3/5):1051-1065.
[3]李晓平.多体系统动力学建模方法及在水下缆索中的应用研究[D].天津:天津大学,2004.
[4]麦继婷.波流作用下悬浮隧道的响应研究[D].成都:西南交通大学,2005.
[5]杜星文,宋宏伟.圆柱壳冲击动力学及耐撞性设计[M].北京:科学出版社,2004.
[6]Wang Pao-Hsii,Yang Chiung-Guei.ParametricStudies on Cable-stayed Bridges and Structures[J].Computers and Structures,1996,60(2):243-260.
[7]唐友刚,海洋工程结构力学[M].天津:天津大学出版社,2008.
[8]秦银刚,周晓军.张力腿型悬浮隧道涡激响应影响因素分析[J].铁道工程学报,2009,24(1):77-81.Qin Yingang,Zhou Xiaojun.Analysis of the Influ-encing Factors of Vortex-induced Vibration Re-sponse of the Submerged Floating Tunnel Supportedby Tension Leg[J].Journal of Railway EngineeringSociety,2009,24(1):77-81.
[9]姚仕明,卢金友,徐海涛.黄陵庙水文断面垂线流速分布特性研究[J].长江科学院院报,2005,22(4):9-11.Yao Shiming,Lu Jinyou,Xu Haitao.Study of Verti-cal Velocity Distribution at Huanglingmiao Hydro-logic Section[J].Journal of Yangtze River ScientificResearch Institute,2005,22(4):9-11.
【数值模拟计算】相关文章：
动态数值模拟05-07
油藏数值模拟技术05-18
油藏数值模拟学习心得04-26
微气泡减少平板摩擦阻力的数值模拟04-16
高超声速电离绕流的数值模拟04-22
小型轴流CPU风扇设计与数值模拟04-29
框架结构建筑物爆破拆除数值模拟研究09-11
黄土辅助隧道中系统锚杆支护作用数值模拟11-27
半导体断路开关电路——流体耦合数值模拟09-12