数值试验范文

数值试验范文（精选10篇）

数值试验 第1篇

在零件磨损数值仿真模型中,要使仿真结果准确就必须建立准确的几何过渡关系,而磨损量的变化规律是一个重要的影响因素,根据磨损机理的不同,磨损量随载荷的变化也不尽相同。

建立准确的几何过渡关系必须考察变量之间的各种关系,而这种关系又是多种多样的,其大致可分为:①确定性关系:即变量之间的关系可以用精确的函数关系来表达;②非确定性关系:也称为相关关系,即在实际的物理研究过程中,有些问题很难由经典物理理论推导出物理量的函数表达式或推导出的表达式十分复杂,为研究需要,就希望能得到这些物理量之间的函数关系的近似表达,即经验公式。回归分析就是研究相关关系的数学工具,它可以采用曲线拟合的方式,用由试验数据得到的物理量之间的近似函数表达式y=f(x)来描述。

本文对经离子渗碳处理的不锈钢的磨损量随载荷变化的规律采用最小二乘法的回归方式来进行数值拟合以验证试验结果,以此为摩擦学系统特性的预测提供一定的帮助。

1试验设计方案

在MM-200型磨损试验机上对不同载荷情况下与对磨滚轮GCr15配对的试样(抛光)进行耐磨性试验,下试样转速n=400r/min,磨损试验用湿摩擦方式。试验开始时,在每个下试样磨损面上滴一滴润滑油(20#机油),转动下试样,使下试样有效工作磨损面均匀敷上一层油膜,并在试验过程中,用自制滴油装置,以1滴/min的量滴入下试样。用光电天平测量在一定载荷、速度下试样磨损后的质量,其精度为0.1mg。

2试验结果及数据处理

2.1 试验结果

试样磨损量试验数据见表1,磨损量曲线见图1。

2.2 试验数据的处理

从图1可看出,试样磨损量的变化趋势比较接近指数函数和幂函数,这时可以借助数学方法将其分别转化为其对应的线性回归曲线y=f(x),从中选出最佳拟合函数来描述试样磨损量的变化趋势。

2.2.1 方法一

把试验数据回归成y=aebx(a>0,b>0)的指数形式,其中,x为载荷值,y为试样磨损量,进行回归计算。

将y=aebx取对数,改写为:lny=lna+bx。并令Y=lny,n=lna,则有:Y=n+bx,其中Y与x为线性关系。

对样本(x1,y1),(x2,y2),(x3,y3),(x4,y4),(x5,y5)用最小二乘法估计a,b值。经计算得出:

undefined。

其中:undefined;undefined。具有undefined;a=en。表2为方法一的回归分析表,将数值代入计算得到:

undefined。

所以,回归表达式为:

y=0.43e0.001 52x。

下面进行回归显著性检验:

回归平方和(S回):S回=b2Sxx=0.231。

离差平方和(SYY):undefinedundefined2=0.273,其中N为样本数。

残差平方和(Qe):Qe=SYY-S回=0.042。

F检验法:undefined。

F≥Fα=0.05(1,3)=10.13,其中α为检验水平。

2.2.2 方法二

把试验数据回归成y=axb(a>0,b>0)的幂函数形式,其中x为载荷值,y为试样磨损量,进行回归计算。

将y=axb取对数,改写为:lny=lna+blnx,令Y=lny,X=lnx,n=lna,则有:Y=n+bX,Y与X为线性关系,用线性回归法确定a,b。表3为方法二的回归分析表,解法同方法一,得到:

b=0.37,n=-3.42,a=0.033。

所以,回归表达式为:

Y=0.033x0.37 。

下面进行回归显著性检验:

回归平方和(S回):S回=b2Sxx=0.221。

离差平方和(SYY):undefinedundefined2=0.273。

残差平方和(Qe):Qe=SYY-S回=0.052。

F检验法:undefined。

F≥Fα=0.05(1,3)=10.13。

通过以上两种方法可以看出,指数函数的F大于幂函数的,其回归分析的置信概率也稍高,从而可知指数函数为最佳拟合函数;从试验数据曲线与采用最小二乘法对试验数据进行数值拟合的拟合曲线的对比曲线图(见图2)看出,指数函数曲线更逼近试验数据曲线。所以,该组摩擦学的磨损量与施加载荷服从指数关系:y=0.43e0.001 52x。

3结论

在MM-200型磨损试验机上用最小二乘法拟合经离子渗碳处理的不锈钢的磨损量随载荷变化的规律,并通过回归分析得到了最佳拟合函数(指数):y=0.43e0.001 52x。

摘要：数值仿真技术因可以显著提高摩擦学研究的效率和摩擦学系统特性预测的准确性而成为当今摩擦学研究领域的重要方法之一,其研究重点是用一定的数学方法建立合适的仿真数学模型。对通过离子渗碳处理的不锈钢的摩擦量随载荷变化的规律进行数值拟合以验证试验结果;在MM-200型磨损试验机上用最小二乘法拟合其磨损量随载荷变化的规律,并得到了最佳拟合函数(指数函数)。

关键词：不锈钢,磨损量,载荷,指数函数,数值仿真

参考文献

[1]汪选国,严新平,李涛生,等.磨损数值仿真技术的研究进展[J].摩擦学学报,2004,4(2):188-192.

[2]陈魁.应用概率统计[M].北京:清华大学出版社,2002.

[3]江亲瑜.零件磨损过程及寿命预测的数值仿真[J].润滑与密封,1997,29(6):29-30.

[4]江亲瑜.数值仿真技术及其在磨损研究中的应用[J].大连铁道学院学报,1997,18(2):18-22.

[5]江亲瑜,李宝良,孙晓云.磨损数值仿真模型中材料磨损率试验研究[J].大连铁道学院学报,2000,21(1):29-32.

数值试验 第2篇

一次特大暴雨天气的数值模拟试验研究

目的 检验各非绝热物理过程参数化方案对模式输出结果的影响,进一步了解和认识MM5模式中不同积云对流参数化方案的性能、特点和适用范围,为合理选择使用模式中的物理过程参数化方案提供依据,提高模式对类似暴雨天气过程的预报能力.方法 运用相互比较和试验分析的方法在不同模式分辨率下,对MM5V3模式中的不同物理过程参数化方案设计18种组合方案,进行比较试验和分析.结果 主要雨带位置对参数化方案并不是十分敏感,但模拟的暴雨中心强度、范围、雨带走向随参数化方案的不同有较大变化,甚至出现虚假的.强降水中心.各组合方案模拟的天气尺度系统水平环流结构差别较小,模拟的中尺度系统结构存在着较大差别,这种差别对定点、定量降水的准确预报产生影响.结论 在实际预报中可根据初始气象场的分型特点选择较合适的组合方案启动数值预报模式,从而达到提高模式对类似暴雨天气过程的预报能力.

作 者：王繁强 戴进 WANG Fan-qiang DAI Jin  作者单位：王繁强,WANG Fan-qiang(陕西省气象科学研究所;陕西省气候中心,陕西,西安,710015)

戴进,DAI Jin(陕西省气象科学研究所)

刊 名：西北大学学报（自然科学版）  ISTIC PKU英文刊名：JOURNAL OF NORTHWEST UNIVERSITY(NATURAL SCIENCE EDITION) 年，卷(期)： 37(3) 分类号：P435 关键词：特大暴雨   数值模拟   天气预报  

数值试验 第3篇

摘要：建立了自然暴晒下汽车乘员舱内的温度场流场模型，考虑车身传热、车内空气流动和人体散热的影响，对乘员舱内传热过程进行模拟计算。通过暴晒升温试验和制冷降温试验分析制冷模式关闭和开启2种工况下乘员舱温度动态变化规律。仿真结果表明，在制冷工况下前排的冷却效果比后排好，但垂直方向上局部温度差异较大，人体表面存在较大温差。试验结果与仿真结果较吻合，并发现前、后排的热不均匀性在冷却风作用下有不同程度的改善，在制冷初始阶段（10 min内）乘员舱的温度变化较为剧烈，对人体热舒适性影响大。研究成果对优化空调设计和提高汽车的乘员热舒适性有重要意义。

关键词：汽车热舒适性；太阳辐射；空调系统；人体散热；车内温度场

中图分类号：U463。83；TB61文献标识码：A

舒适环保是消费者购买汽车的一个重要的考虑因素，汽车乘员舱内的热舒适问题已经成为广大消费者、研究学者和汽车厂商关注的热点和焦点问题。车内小空间热环境和乘员热舒适性密切相关，车内的空气温度、相对湿度、气体流速和平均辐射温度等环境因素直接作用于人体体表，影响乘员的主观热舒适感觉。汽车空调是调节车内温度、湿度、空气清洁度和空气流动性的主要手段，也是提高乘员乘坐舒适性的主要途径。据统计，在城市行车中，85%的旅程在18 km以内，行车耗时约15～30 min＼[1＼]。文献＼[2＼]和文献＼[3＼]的研究表明，在不考虑人体散热的情况下，空调制冷模式开启后的10 min和15 min内，车内温度变化相对剧烈，之后车厢热环境变化较小。可见，人们从外部环境进入车内后短时间内的热舒适感觉对整车的乘坐热舒适性评价是相当重要的。汽车空调系统应在尽可能短的时间内将乘员舱内的热环境调节至人体的舒适范围，还应避免因温度变化过快引起人体的热不舒适感。汽车空调的设计不仅要着眼于缩短使车内热环境达到稳定状态的时间和保证达到稳态后乘员的热舒适性，还应保证在升温或降温的动态过程中乘员的热舒适性。因此，研究自然暴晒中空调开启后乘员舱内热环境的变化规律，对提高整车的乘坐热舒适性、增强汽车的市场竞争力有重要意义。

太阳辐射、空调系统和人体散热是影响乘员舱热环境的主要因素。国内外学者就其中一个或某几个因素对乘员舱热环境的影响进行了研究。文献＼[4＼]将车身壁面温度设为定温，对比不同局部制冷策略下驾驶员的热感觉和热舒适性。文献＼[5＼]通过数值方法研究了空调出风的流量、速度和温度等流量参数和出风口面积、形状等几何参数对驾乘人员局部和整体热舒适性的影响。文献＼[6＼]比较有/无搭载人员两种情况下，空调冷却风射出角度对车内温度场和流场的影响。文献＼[7＼]通过试验分析湿度对车内温度及乘员热舒适性的影响。国内学者也针对轿车、客车等不同车型开展相关研究＼[8-10＼]。已有学者根据乘员热舒适性和车内温度场的要求来改进空调系统[11-12]。如芦克龙等＼[13＼]将重型货车空调系统和乘员舱作为一个整体，加入驾驶员模型，根据乘员的热舒适性对空调的风道设计进行改进，达到改善驾驶员热舒适性的目的。黄木生＼[14＼]构建了适用于评价微型车乘员舱乘员热舒适度的准则，在某款微型汽车的开发中得到应用。

本文结合数值仿真和物理试验对太阳辐射、空调系统和人体散热共同作用下乘员舱内热环境和人体体表温度的热响应过程进行研究，建立了乘员舱内的温度场流场模型，包括乘员舱传热模型、车内空气流模型和人体散热模型，仿真得到乘员舱的稳态温度场和人体表面温度分布；在空调关闭的情况下对试验车辆进行1 h的空车自然暴晒试验，紧接着试验人员进入副驾驶座开始30 min的空调制冷试验，研究成果为设计高效的空调控制策略提供科学依据和理论指导，实现减小乘员舱热负荷、提升乘员的热舒适感觉的目的。研究成果对构建节能、舒适、安全的车内乘坐环境具有深远意义。

1乘员舱内温度场流场模型

1。1物理模型

以某轿车车厢简化模型为研究对象，见图1。该车厢几何模型长2 610 mm，宽1 715 mm，高1 780 mm；副驾驶座上安装人体模型；空调系统的4个出风口对称布置在前排的仪表板上，中间2个、两侧各一个；空调回风口位于副驾驶座下方。

1。2数学模型

计算流体力学理论以3个流动基本方程为基础，遵守质量守恒定律、动量守恒定律和能量守恒定律。乘员舱内空气的流动和热传递满足质量守恒方程、动量守恒方程和能量守恒方程。直角坐标系下这3个基本方程可分别表示为：

ρt+ρuixi=0； （1）

（ρui）t+（ρuiuj）xi=-Pxi+（τij）xj+Fi； （2）

（ρT）t+（ρuiT）xi=xi（kcpTxi）+ST。（3）

式中：ρ为密度；t为时间；x为坐标分量；i，j=1，2，3，表示直角坐标中相互垂直的坐标方向；ui为速度矢量；P为流体单元受到的压力；τij为流体单元表面上的粘性应力分量；Fi为流体单元上的体力；T为温度；k为流体的传热系数；cp为比热容；ST为粘性耗散项。

人体因新陈代谢会向外散发热量，本文将乘员简化成散热恒定的热源，热流密度为150 W/m2，不考虑服装热阻。

1。3边界条件

边界条件的设置对数值仿真的结果有很大影响，应尽量贴合实际。参考已有研究成果，并考虑汽车实际的使用工况，设置的边界条件具体如下：

1）4个空调出风口采用速度入口边界条件，仪表板两侧的出风风速为7 m/s、温度为290 K，湍动能k=0。184，湍动耗散率ε=1。43；仪表板中间的2个出风口风速为10 m/s、温度为290 K，湍动能k=0。175，湍动耗散率ε=1。33；

2）空调回风口为压力出口，P=101 325 Pa；

3）车厢地板和车身定义为无滑移的壁面边界条件；

4）车窗、天花板、仪表板及座椅参与车厢内的对流、传导及辐射热交换过程，设为外部辐射边界条件；

5）选用DO模型计算太阳辐射强度，认为车内空气是不可压缩的。

2乘员舱内热环境仿真分析

仿真得到乘员舱内稳态温度场，见图2。人体模型的表面温度大多处于306～310 K；两侧手臂由于受空调冷风直吹，温度较低，温度为298～300 K；两侧肩膀和大腿根部的温度最高，约318 K。汽车内表面的最高温度出现在后窗台的边角位置，约322 K；仪表板靠近前风挡玻璃的位置也呈现出较高的温度。截取4个截面，研究乘员舱内的空气温度分布。具体位置是：截面Ⅰ，前排座椅侧前150 mm垂直剖面；截面Ⅱ，后排座椅侧前150 mm垂直剖面；截面Ⅲ，副驾驶员位置纵向中剖面；截面Ⅳ，副驾驶员胸部水平剖面，见图3。

图4为各截面温度分布情况，4个截面中温度波动都保持在4 K以内，乘员舱内车内稳态温度场的温度均匀性较好。截面I和截面Ⅲ均显示前、后乘员舱的脚部温度都比胸部和头部温度高出2 K；截面I显示，两侧车门附近、高度相当于乘员小臂位置的空气温度低于其他位置的空气温度；截面Ⅱ中，各个位置的温度差异较小；截面Ⅳ显示乘员胸部水平面上车厢内各处的温度值，前仪表板温度较高，其他位置的温度基本相同。

太阳辐射、空调送风方向和车内结构设计影响乘员舱内不同位置的冷却效果。驾驶座和副驾驶座的乘员外侧小臂位置受到冷却风的直吹作用，冷却效果明显，温度比肩部和脚部低。太阳光通过车窗玻璃直接照射在仪表板和后窗台上，二者吸收的热量较多；仪表板表面位于空调出风口上方，冷却风只能通过自由扩散达到该区域，对流冷却强度明显不足；空调出风口全都安装在前排的仪表板上；冷却风受到人体和座椅等固体区域的阻碍，到达后排时速度减小且因吸热而导致温度上升，对后排的冷却效果较差，在后窗台贴近后车窗的狭小空间内形成气流死角，冷却效果较差。冷却风流到后排时，方向性减弱而向各个方向发散，因而各个位置的冷却程度没有明显的差别。

确测

量分布在乘员舱各处的空气温度，在乘员的头部、胸部和脚部位置都布置了温度传感器，前后排各3个测量点；车内零部件表面温度的7个温度测量点分布在仪表板、天花板、前风挡玻璃、驾驶员侧的坐垫和靠背、副驾驶员侧的坐垫和靠背，具体的测

点位置见图6。数据采集器每5 s采集一次温度数据。

3。1试验结果与仿真结果对比验证

各零部件表面温度和人体体表温度的试验结果与数值仿真结果对比如表1和表2所示。除了座椅外，其他位置的试验值和仿真值误差均控制在10%之内，整体上仿真计算结果与试验结果吻合较好，可认为符合工程计算的要求，验证了前述计算模型选择和边界条件设置的合理性和准确性。

（a）试验车辆

（b）部分测点布置

图5 试验现场

Fig。5Test site

（a） 主视图

（b） 俯视图

图6温度测点分布

Fig。6Temperature sensors distribution

表1零部件表面温度

Tab。1Parts surface temperature

名称

试验结果

/℃

仿真结果

/℃

误差

/%

仪表板

30。6

32。1

-4。6

前风挡玻璃

25。2

23。4

7。1

座椅

28。1

24。9

11。1

天花板

27。4

25。9

5。5

表2人体表面皮肤温度

Tab。2Human body skin temperature

名称

试验结果

/℃

仿真结果

/℃

误差

/%

头部

33。3

32。1

3。6

背部

36。6

34。4

6。0

胸部

32。6

34。4

-5。2

腹部

33。1

31。2

2。7

胳膊

29。3

27。6

5。8

大腿

34。5

35。7

-3。4

小腿

32。7

34。1

-5。8

3。2试验结果分析

3。2。1暴晒升温阶段的试验结果分析

利用Matlab对采集到的温度数据进行曲线拟合。暴晒升温阶段车内没有搭载人员、空调关闭，乘员舱温度主要受太阳辐射的影响。

车内空气平均温度和零部件表面平均温度的变化如图7所示，该组数据通过对测点温度进行算术平均得到。在太阳下暴晒1 h后，车内空气平均温度和零部件表面平均温度分别达到了44 ℃和49 ℃。各零部件表面温度变化如图8所示。由于仪表板受到通过车窗玻璃的太阳光直射，所以其温度始终高于前挡风玻璃和天花板，最高温接近65 ℃；前挡风玻璃透射太阳光，吸热量较小，最高温约50 ℃；无法接受太阳直射的天花板最高温只有46 ℃。驾驶员侧坐垫和靠背的温度都高于副驾驶员侧对应位置的温度，温升速率也快于副驾驶侧。太阳光主要从驾驶员侧车窗玻璃进入车厢内部，导致正驾驶和副驾驶两侧的座椅表面温度差异较大。图9和图10分别为乘员舱前排和后排的垂直方向上的空气温度变化。前排和后排垂直方向的温度场分布基本相同，温度从低处往高处逐渐增加，脚部、胸部和头部温度分别约为37 ℃，47 ℃，19 ℃。随着暴晒时间的增长，垂直方向上的温差也逐渐拉大。暴晒初始，头部位置和脚部位置的温度差约5 ℃，而在暴晒后期，温差增大到12 ℃左右。造成这种差异的主要原因是车门阻挡部分太阳辐射，使脚部位置空气的受光量较小，所以温度上升幅度较小。在暴晒结束时，仪表板表面的温度最高（65 ℃），前、后排脚部位置的温度最低（约37 ℃），说明车内热环境整体温度较高、区域温度差异性十分明显。如果不能在短时间内有效地降低乘员舱内的温差和高温状态，将会对前后排的乘坐热舒适性造成极大影响。

3。2。2制冷降温阶段的试验结果分析

在制冷降温阶段乘员舱热环境受太阳辐射、冷却风和人体散热3个因素的综合影响。

如图11所示，在整个降温过程零部件表面平均温度始终高于空气的平均温度。在制冷的初始阶段（10 min内），零部件和车内空气的平均温度均出现大幅度下降，分别降低10 ℃和17 ℃，零部件的表面温度变化曲线相对平缓。在随后的20 min降温过程中，零部件表面温度和车内空气温度变化较小，可认为车内的热环境达到平衡状态。

4结论

结合数值仿真和物理试验对汽车乘员舱内温度场变化规律进行研究，建立乘员舱内的温度场流场模型，计算制冷工况下乘员舱内的稳态温度分布，通过试验验证仿真结果，并通过对比空调制冷关闭/开启2种情况下汽车乘员舱热环境的动态变化规律，分析太阳辐射、人体散热和空调系统对乘员舱温度场稳态分布和动态变化规律的影响。

太阳辐射对乘员舱热环境的影响在空调系统关闭时较为突出，太阳光照射的时间越长，乘员舱垂直温差越大，暴晒升温过程中，头部和脚部位置的温差从5 ℃增加到12 ℃，会给人体造成强烈的瞬时热冲击，引起乘员的热不舒适感。人体会阻碍空调风的流动，但人体散热对乘员舱温度分布的影响可以忽略。

空调系统的出风口位置和送风方向对乘员舱热环境和乘员体表温度有决定性的影响。空调出风口全布置在前排仪表板上，对前排的冷却效果优于后排，能不同程度地降低前、后排垂直温差，对后排空间的冷却程度较一致；受到冷却风直吹的人体部位呈现较低温度，乘员体表出现明显温差；在制冷的初始阶段（<10 min），乘员舱内温度变化剧烈，这段时间内乘员可能因冷却效果不佳、冷却时间过长或降温过快而产生热不舒适感。

参考文献

[1]ALAHMER A， MAYYAS A， MAYYAS A A， et al。 Vehicular thermal comfort models， a comprehensive review＼[J＼]。 Applied Thermal Engineering， 2011， 31（6/7）：995-1002。

[2]HAN T， CHEN K。 Assessment of various environmental thermal loads on passenger compartment soak and cooldown analyses＼[R＼]。SAE Technical Papers，2009-01-1148。

[3]崔津楠。汽车乘员舱空调制冷CFD瞬态分析＼[J＼]。计算机辅助工程，2012，21（6）：54-57。

CUI Jinnan。 Transient CFD analysis on airconditioner cooling of automobile occupant cabin＼[J＼]。Computer Aided Engineering，2012，21（6）：54-57。（In Chinese）

[4]KAUSHIK S， CHEN K， HAN T， et al。 Microcooling/heating strategy for energy efficient HVAC system＼[J＼] 。 SAE International Journal of Materials and Manufacturing，2011，4（1）：853-863。

[5]CURRLE J， MAU J。 Numerical study of the influence of air vent area and air mass flux on the thermal comfort of car occupants＼[C＼]//SAE Technical Paper。2000-01-0980。

[6]VIVEK P， ABDUL N， NAGPURWALA Q H。 Numerical studies on the effect of cooling vent setting and solar radiation on air flow and temperature distribution in a passenger car＼[R＼]。 SAE Technical Paper，2009-28-0048。

[7]ALAHMER A， ABDELHAMID M， OMAR M。 Design for thermal sensation and comfort states in vehicles cabins＼[J＼]。 Applied Thermal Engineering， 2012，36（1）：126-140。

[8]肖红林，李洪亮，王远，等。考虑人体散热的轿车座舱内流场数值仿真研究＼[J＼]。计算机仿真，2011，28（6）：320-325。

XIAO Honglin，LI Hongliang，WANG Yuan，et al。 Investigation on numerical simulation of flow field in car cabin with demic heat dissipation＼[J＼]。 Computer Simulation， 2011，28（6）：320-325。 （In Chinese）

[9]熊可嘉，王伟，张万平。考虑太阳辐射的轿车客舱流场与温度场的数值计算＼[J＼]。制冷与空调，2008，22（6）：12-15。

XIONG Kejia，WANG Wei，ZHANG Wanping。 Automotive passenger compartment HVAC simulation with solar radiation＼[J＼]。 Refrigeration and Air Conditioning，2008，22（6）：12-15。（In Chinese）

[10]张炳力，薛铁龙。考虑太阳辐射的乘员舱热舒适性仿真分析＼[C＼]//2013中国汽车工程学会年会论文集。北京：北京理工大学出版社，2013：1344-1348。

ZHANG Bingli，XUE Tielong。 Simulation and analysis of thermal comfortable inside a passenger compartment with solar radiation＼[C＼]//2013 SAEChina Congress Proceedings。 Beijing：Beijing Institute of Technology Press，2013：1344-1348。（In Chinese）

[11]张文阁。基于人体热舒适性的汽车空调优化设计＼[D＼]。长沙：湖南大学机械与运载工程学院，2012。

ZHANG Wenge。 Optimal design of vehicle air condition based on thermal comfort of human body＼[D＼]。Changsha： College of Mechanical and Vehicle Engineering，Hunan University，2012。（In Chinese）

[12]吴志武。皮卡车空调和太阳辐射对乘客舱内热流场的影响研究＼[D＼]。广州：华南理工大学机械与汽车工程学院，2011。

WU Zhiwu。 Research on the influence of the air condition and solar radiation on indoor thermal flow field for a pickup＼[D＼]。 Guangzhou： School of Mechanical and Automotive Engineering，South China University of Technology，2011。（In Chinese）

[13]芦克龙，谷正气，贾新建，等。某重型货车空调系统对乘员舱热舒适性影响的分析与改进＼[J＼]。汽车工程，2011，33（2）：162-166。

LU Kelong，GU Zhengqi，JIA Xinjian，et al。 Analysis and improvement of the effects of airconditioning system on the thermal comfort of passenger compartment in a heavy truck＼[J＼]。 Automotive Engineering，2011，33（2）：162-166。（In Chinese）

[14]黄木生。基于CFD的微型车乘员舱乘员热舒适度分析＼[D＼]。长沙：湖南大学机械与运载工程学院，2013。

模具刃口的激光淬火试验与数值模拟 第4篇

激光淬火应用于模具的局部强化，可显著提高表面硬度、改善刃口等部位的耐磨性能，提高模具的使用寿命;可实现低成本材料替代高成本材料，改进模具制造和使用方式;克服传统模具因加工和热处理带来的诸多弊病，提高模具整体质量，大幅度降低模具制造成本，简化制造工艺，缩短制造周期，增强模具企业的竞争力[3]。从已发表的论文来看，绝大部分学者采用激光器对金属表面进行激光扫描，对硬化区进行X射线衍射谱分析、光学显微分析和硬度测试[4,5,6,7,8],这种方法成本高，对变换淬火工艺参数或材料不具备通用性。王秀凤[9]等学者曾采用试验与数值模拟相结合的方法，对激光淬火工艺进行了初步研究，从而减少了试验次数，提高了效率。

本文以模具常用材料45钢为例，采用CO2激光器对试件边缘(模拟模具刃口)进行激光淬火，实验测定淬火层硬度、深度和宽度，并通过有限元分析软件模拟了激光淬火过程，得到了淬火过程中的温度场、应力场及位移场等信息，经验证，模拟结果与试验相吻合。结论表明，通过改变淬火的工艺参数，模拟不同情况下淬硬层深度变化，可为实际应用提供参考依据。

1 试验

试验在中科院力学所表面改性实验室进行。试验设备为2kw连续输出式CO2激光器。试件材料为45钢，尺寸为50mm25mm25mm。为了对模具刃口的强化效果进行研究，试验时将试件平放在工作台上，沿试件的棱边(此处为Y方向)进行激光连续扫描，扫描前将扫描区域涂黑，以提高材料对激光的吸收率，如图1所示。激光器输出功率为900W，光斑直径为5mm，扫描速度为13mm/s。

扫描完毕，在光斑作用处沿X-Z切取截面取试样，做金相试验。在制备金相试样时，为防止磨坏扫描区，在其上表面粘贴薄金属片。试样打磨、抛光后，用3%的硝酸酒精腐蚀，通过NEOPHOT21金相显微镜观察其组织如图2所示，初步测定淬硬层最大深度为0.53mm，半宽度为2.65mm，如图3所示，图中刻度尺的最小刻度为0.1mm。

采用硬度测量仪测量试件相变区不同深度处的硬度值，如图4所示，根据测量值绘制硬度随深度变化的曲线如图5所示，图中测量点的垂直间距均为0.05mm。

图5表明，在距试件表面0.45mm以内为相变硬化区，硬度值在60-64HRC范围内变化，远大于常规淬火所达到的最高硬度值(55HRC);0.45～0.55mm深度处为过渡区，硬度值明显下降，大于0.55mm处则为基体材料。因此有效淬硬层的深度可取为0.45mm，比由金相图读取的数值(0.53mm)略小。

由图5可以看出，淬火区的硬度很高且硬化层硬度基本无梯度变化，十分有利于提高材料的寿命。而常规淬火硬化层的硬度从表面到内部有明显的下降梯度，表面硬度最高，耐磨性最好，但随着零件正常工作时的相对硬度，表面将逐渐被磨去，其相对运动接触面的硬度值随之逐渐降低，磨损随之加剧，最终导致零件因磨损量过大而失效。

2 数值模拟

2.1 模型的建立

利用大型非线性有限元软件MSC.Marc模拟激光淬火过程。模型尺寸为50mm25mm25mm，其工艺参数和边界条件与试验相同，环境温度为27.5℃;材料的力学性能参数和热物性能参数均是温度的函数，取自手册[10,11];通过Marc软件的子程序接口，利用FORTRAN编写FLUX用户子程序来模拟移动的圆形激光光斑的自动加载。考虑到激光束辐照后在试件作用区产生强烈的温度梯度，故采用三维八节点六面体减缩积分单元以达到较高的精度;为提高计算效率，单元划分时在激光扫描区网格加密，最小单元尺寸为0.5mm0.5mm0.5mm，远离扫描区网格逐渐稀疏，如图6所示。

2.2 温度场模拟结果

通过热-机耦合模拟了45钢的激光淬火过程。45钢的相变临界温度为727℃，且随加热速度的增加而升高，可由试验值估取。选取模型在任一光斑作用中心处的横截面，如图7所示，光斑中心处表面最高温度为1065.9℃;利用最小刻度为0.1mm的刻度尺测量750℃处硬化层深度为0.44mm、宽度为4.12mm，如图8所示，与试验测得值较吻合，认为750℃处作为相变的临界温度较为合理。

任取激光扫描路径上的一点，其温度随时间变化的曲线如图9所示。

图9可以看出，激光淬火时，其加热和冷却速度极快，加热相变区的温度梯度很大，在极短的时间内生成的大量细小的奥氏体形核，在快速冷却时来不及长大和扩散，致使相变区晶粒细化。所以激光淬火时只要表面温度低于熔点，不使其熔化即可。为了增加硬化层深度，在表面不熔化的前提下，可以加大能量注入，如增加激光输出功率或降低扫描速度。如图10所示等温线图是在其它参数不变的条件下，采用功率1000W，速度8mm/s进行模拟计算的结果，对比图8可以看出淬火层深度明显增加，约为0.7mm，宽度基本不变，最高温度上升至1189℃。

2.3 应力场模拟结果

激光扫描完毕后，在板料中存在残余应力，它是由热应力和相变应力共同作用的结果。通过模拟计算可以方便地得到残余应力的大小及分布，如图11所示。在激光扫描的部位及其垂直下方的区域产生压应力，压应力的周围尤其是靠近试件侧面的部位产生较大的残余拉应力，此外试件在远离扫描区的底面棱角处产生压应力。总体来看是压应力和拉应力相互交织在一起的比较复杂的分布，众所周知残余压应力可以提高材料的可靠性和使用寿命，而残余拉应力则将导致裂纹的产生及扩展。所以，激光淬火在选定工艺参数的时候，要注意控制残余应力的大小及分布。

2.4 位移场模拟结果

激光淬火过程中，由于热应力和相变应力的存在，扫描完后试件将产生塑性变形。扫描完后试件的位移分布如图12(a)所示;在试件上表面变形最为明显的区域取一点，画出该点位移随时间变化曲线，如图12(b)所示，可以看出，激光经过扫描点时，在其表面产生较大的塑性变形;激光扫描过后，试件快速冷却，变形大幅回复，最终的塑性变形却很小，位于激光扫描表面的最大位移为0.009756mm，相对于试件的几何尺寸可忽略不计，由此可见，激光淬火后硬化部分的热影响区较常规淬火方法要小很多，产生的热应力也小的多，在实际应用时可将模具刃口处加工到最终尺寸，激光淬火后就可以投入使用，无需后序精加工。

3 结论

数值模拟与试验的相互印证表明，通过建立较准确的数值模型，采用热-机耦合计算，获得全面的激光淬火过程信息。通过温度场预测了硬化层的深度和宽度，试件在激光淬火后表面硬度可以达到模具钢的使用要求;通过应力场得知在试件中存在复杂交织的残余拉、压应力，为达到预期的性能要求可以选取合适的激光工艺参数;通过位移场得知在热应力和相变应力的共同作用下，试件最终产生很小的塑性变形，可以忽略不计。

研究表明，激光淬火可显著提高模具表面的硬度、细化晶粒，增加模具的耐磨性，提高模具的使用寿命;激光淬火的可控性可提高产品淬火质量的稳定性，其高效率可缩短产品的制造周期。同时，材料性能的大幅提高使低成本材料替代高成本材料成为可能。

参考文献

[1]张国顺.现代激光制造技术[M].北京:化学工业出版社,2006:164-181

[2]关振忠.激光加工工艺手册[M].北京:中国计量出版社,1998:133-134

[3]冯荣元.模具激光强化技术[J].模具制造.2006,11:76-80

[4]李同道,王勇,赫庆坤,等.45钢多次激光相变硬化组织与性能研究[J].中国表面工程.2007,20(2):33-36

[5]王大承.HT300链轮齿表激光强化技术研究[J].应用激光.2001,21(6):403-406

[6]刘其斌,魏远翔,李海,等.T10钢宽带激光相变硬化的组织与性能研究[J].现代机械.2003.6:85-87

[7]周健,温宗胤,李新华.激光相变硬化在Cr12汽车模具材料表面强化中的应用研究[J].模具工业.2007,33(4):67-70

[8]吴钢,刘宝,宋光明.激光工艺参数对表面硬度影响规律的研究[J].机械工人·热加工.2005,1:59-61

[9]王秀凤,吕晓东,陈光南,等.激光强化温度场的数值模拟与校验[J].激光技术.2004,28(4):162-165

[10]谭真,郭广文.工程合金热物性[M].北京:冶金工业出版社.1994.9.

非定常喷流有限差分方法数值试验 第5篇

非定常喷流有限差分方法数值试验

以欠膨胀自由喷流初期流动为例,采用Euler方程和Beam-Warming/TVD有限差分格式,对比分析了Beam-Warming AF方法,隐式亚迭代方法和简化Runge-Kutta五步格式的非定常流场描述能力,结果表明:(1)隐式近似因式分解方法基本上可以描述非定常流动现象;(2)隐式亚迭代一阶时间精度格式会导致流场结构的变化,其精度可能是不足的.;(3)隐式亚迭代二阶时间精度与简化Runge-Kutta五步格式的计算结果一致,可以认为是计算非定常问题的适当方法.

作 者：马汉东 袁宏生 作者单位：北京空气动力研究所,北京,100074刊 名：计算物理 ISTIC EI PKU英文刊名：CHINESE JOURNAL OF COMPUTATIONAL PHYSICS年，卷(期)：18(4)分类号：V211关键词：非定常流 有限差分方法 数值试验 自由喷流

中尺度风资源数值模拟试验分析 第6篇

风能是一种可再生、无污染、能量大、前景广的能源, 大力发展清洁能源是世界各国的战略选择。风能资源评估是风能开发利用的关键环节, 它是制定风能规划、风电场选址、风电功率预测、风机选型、机组布置和发电量计算的重要基础。而风电场中由于测风塔分布数量有限, 仅依靠单纯测风塔资料很难对风能资源空间变化格局做出准确评价, 因此选用模型模拟风能资源分布已成为一种重要而被广泛接受的研究方式[1,2,3]。

1 风场数值模拟实施方案

假设风电场位于近海区域, 风场内没有测风塔, 附近有一座海上测风塔和若干海岸上陆地测风塔, 且它们离风场均较远, 用其对本风场进行风资源评价存在较大误差。因此, 本阶段采用数值模拟方法, 对此类型风场区域风资源进行垂直和水平双向计算, 并通过引入实测数据, 对其进行检验和校正, 能较大提高拟建区域的风能资源评估精度, 为风机选型、机位布置和发电量计算打下坚实基础, 提高项目收益, 降低项目风险。针对此类风场, 所采用的数值模拟方法为WRF (Weather Research and Forecast) 模式。

1.1 WRF数值模式简介

WRF (Weather Research and Forecast) 模式[4]系统是由美国多家研究部门及大学的科学家共同参与开发研究的新一代中尺度预报模式和资料同化系统。WRF模式是一个完全可压非静力模式, 控制方程组都写为通量形式。垂直坐标采用地形跟随静力气压垂直坐标, 网格形式采用Arakawa C格点, 时间积分采用时间分裂显示方案来解动力方程组, 即模式中垂直高频波的求解采用隐式方案, 其它波动则采用显示方案。显示时间积分方案提供了2阶和3阶Runge-Kutta时间积分。

WRF模式集成了多年来在中小尺度天气动力学方面的研究成果。WRF模式采用Arakawa-C网格、三阶Runge-Kutta时间积分、高阶平流方案、标量守恒差分方案等技术, 可使模式耗散项更小、精度更高, 从而使模式的最佳水平分辨率可以达到1 km~10 km。针对中小尺度天气特点, WRF模式对辐射过程、边界层参数化过程、对流参数化过程、次网格湍流扩散过程及微物理过程等进行了改进和优化处理, 发展出了丰富而成套的物理过程选项, 模式应用技术人员可根据当地的地形、气候特点加以选择, 从而取得较好的数值预报和数值模拟效果。

在软件设计方面, WRF模式应用了继承式软件设计、多级并行分解算法、选择式软件管理工具、中间软件包 (连接信息交换、输入/输出及其它服务程序的外部软件包) 结构, 具有单重、多重及移动嵌套网格功能。

因此, WRF模式特别适合于为用户提供高精度和高分辨率数值模拟结果。近年来, 该模式在世界各地不仅被广泛应用于区域数值天气预报业务, 还被广泛应用于大气环境影响评估与预测、风能资源评估与预报及城市或区域发展规划等领域中。

1.2 模拟方案设计

根据用户需求和当前数值模式的最新研究结果, 在尽可能提高水平分辨率的基础上, 我们采用了三重嵌套方案, 模式最细网格水平分辨选择了相对经济、安全和可靠的4 km最高水平分辨率, 垂直方向分为32层, 并适度增加了边界层内的垂直层数。用户关心区置于最细网格区域的中心, 并适度放大最细网格区以避免侧边界带来的影响。

在考虑了关注区的气候特点和模式分辨率的基础上, 我们对模式所采用的物理过程加以选择。模式中所采用的物理过程包括:显式云降水过程、积云参数化过程、浅对流参数过程、陆面过程、边界层物理过程、辐射物理过程、土壤温度模式等。

模式初始化的大尺度资料采用NCEP发布的水平分辨率1010的全球格点分析资料。首先用WPS模块将地形、植被、土壤等静态资料及大尺度背景资料分析到模式区域, 形成模式初猜场, 然后通过obsgrid模块把气象探空、船舶和地面观测资料分析融入初猜场, 形成模式的初始场。粗细网格区域的时间积分步长分别为180 s、60 s和20 s, 模式积分时间为36 h, 每小时输出一次积分结果。

1.3 资料后处理

由于模式生成资料无法被用户直接使用, 需对模式资料进行后处理, 使其成为用户可直接使用的资料。资料后处理一般有分析和诊断两种方法。CALMET模式是一种诊断模式, 有人使用CALMET模式对WRF模式资料进行加工。这种方法的优点主要是可以引入更高分辨率的地形资料, 并通过风的最小化无辐散迭代的办法, 对描述由于地形所产生的局地风阻挡、下坡风等有一定好处。但是在本项目工作中, 我们发现在不引入新观测资料情况下, 无辐射迭代在某些点产生虚假源汇, 从而造成新误差源。考虑到本项目关心区位于海上, 地形简单, 因而, 本项目中我们采用客观分析的办法, 把WRF模式数据分析到1 km1 km水平格点的相应高度层上, 形成最终用户数据。具体分析方法如下。

1.3.1 水平差值方法

水平差值方法采用Cressman客观分析方法, 对于任意要素H, 有:

式 (1) 中, H为任意气象要素;n为观测台站数;s为从1到n自然数;Ws为权重系数;hs为第s个观测台站观测的气象要素值。

变量h为WRF模式的相应变量, W为权重系数, W定义为:

式 (2) 中, R为最大扫描半径, m, 取为1.5个WRF模式格距;r为WRF模式点到分析点的距离, m;n为扫描半径内的所有格点数。

1.3.2 垂直差值方法

近地面层风速随高度的变化有指数律风廓线和对数律风廓线两种。研究表明, 指数律风廓线使用范围更广泛, 且适用高度也较高。本项目采用指数差值方法。在某一高度z上, 风速v为:

式 (3) 中, v为在高度z上的风速, m/s;v1和v2分别为z高度的上、下层风速, m/s;z1和z2是上、下层的对应高度, m。

2 模拟成果分析

通过以上数值模拟方法对本区域风资源进行模拟, 网格分辨率为1 km1 km, 最终得到风场区域每1 km2的逐小时风资源数据。风速包括10 m、30 m、50 m、70 m、90 m和100 m共6个高度。不同高度风速变化梯度见图1。

由图1可知, 离海边越远, 海上风速越大, 且风速变化梯度较小, 这主要是因为海面粗糙度较小, 风经过海面时, 风速衰减并不明显, 模拟结果符合风速的一般变化规律。

根据数值模拟结果, 找出距离1#和2#测风塔最近的数据点, 其编号为2818和5188。二者同期数据统计结果如表1, 由表1可知, 数值模拟结果与实测数据的100 m同期风速相差分别为0.2 m/s和0.17 m/s;考虑到测风仪器本身误差, 可判断二者风速基本相当。以上结果验证了本次模拟的准确性。

3 结语

介绍了WRF (Weather Research and Forecast) 数值模式模拟方法, 并采用WRF模式对近海风场进行模拟分析。模拟结果较好地反映了风速变化特征, 并通过对模拟数据与实测数据进行对比分析, 两者误差较小, 模拟成果满足设计要求。因此, 利用该数值模式可对近海风电场风资源进行模拟分析, 可有效解决近海风电场风资源资料分辨率不足的问题, 为近海风资源开发利用提供有力的科学依据。

摘要：应用中尺度数值模式模拟地区的风资源状况是一种比较先进的方法。该模拟成果对大范围区域风能资源的宏观评估和风电场宏观选址具有很好的参考价值。着重介绍了WRF (Weather Research and Forecast) 数值模式模拟方法, 并采用该模拟方法对近海风电场的风资源进行了模拟。

关键词：数值模式,WRF,风电场,风资源

参考文献

[1]程麟生.中尺度大气数值模式发展现状和应用前景[J].高原气象, 1999, 18 (3) :350-360.

[2]张华, 孙科, 田玲, 等.应用WRF模型模拟分析风力发电场风速[J].天津大学学报, 2012, 45 (12) :1116-1120.

[3]张云海.用于风电场选址的风能资源评估软件[J].气象科技, 2004, 32 (1) :1-4.

数值试验 第7篇

颗粒材料广泛存在于自然界中,颗粒的形状千差万别,颗粒间的接触特性、约束特性和力学特性差异显著。为了从细观角度描述颗粒间的相互作用,需要对颗粒进行简化处理。本实验采用离散元法将单个颗粒材料简化为圆形,模拟了颗粒体之间具有一定约束作用的双轴压缩试验,得到了各指标之间的关系。

在利用离散元法对颗粒材料进行数值模拟的研究中,周忠等利用离散元法模拟了宕渣碾压过程,解释了宕渣碾压效果与碾压遍数和摊铺厚度的关系,确定了最佳摊铺厚度(30cm)和碾压遍数(5～6遍)[1]。蒋明镜等将考虑颗粒抗转动作用的颗粒接触模型导入PFC2D中,数值模拟了砂土的双轴试验,研究了颗粒抗转动作用对砂土力学行为的影响,结果表明,颗粒抗转作用越强,砂土强度越高,当抗转系数为0.4时,松砂出现软化及剪胀现象[2]。王桂萱等利用颗粒离散元法对二维颗粒体进行了双轴压缩的数值模拟,将离散系统中的变形与连续体中相应的变形相联系,用于计算颗粒离散元法中的平均应力,比较了不同颗粒形状、颗粒间摩擦系数对颗粒体微观-宏观力学性能的影响,并分析了不同摩擦系数下颗粒间接触力的法向分量和切向分量对剪应力的贡献[3]。张翀等采用二维离散单元法,建立了4种不同形状的颗粒,研究了颗粒形状对模拟双轴试验的影响以及4种不同颗粒试样宏观特性与颗粒细观参数的关系,结果表明,在其它细观参数相同的情况下,其中3种异形颗粒显著提高了颗粒试样的剪切强度[4]。周健等采用颗粒流程序数值模拟了砂土试样的双轴试验,通过比较模拟结果和室内试验实测结果发现颗粒流方法能较好地模拟室内试验,并给出了不同颗粒单元参数时砂土试样的宏观性质[5]。叶勇等通过模拟单轴压缩、拉伸及双轴压缩试验,分析了岩石的离散单元法中各参数的关系,模拟结果表明,由杨氏模量、泊松比与颗粒接触刚度的关系可以反向确定接触刚度,法向接触刚度与法向连接强度的比值对岩石离散单元模型分析的影响很小,法向与切向连接强度的比值影响裂纹的破坏方式,而围压控制双轴试验时裂纹的类型[6]。

本实验利用二维离散元法以圆形刚性体表征颗粒材料,并对颗粒间赋予粘结,研究了不同粒径和不同围压作用下粘结颗粒材料双轴试验的响应结果,包括粒径对材料弹性模量、泊松比、强度以及裂缝初始应力的影响,组种试件出现裂缝时围压对裂缝初始应力的影响,并以荷载作用后试件内部颗粒间的接触力和位移矢量表征试件的变形和颗粒运动状态。

1 双轴试验数值模拟

1.1 双轴试验环境模型

在利用离散元法数值模拟双轴试验的过程中,首先在宽度w、高度h的矩形范围内生成试件,其中包括4个墙体,上下墙体作为承载板施加加载速率(vp),左右墙体用来施加围压,本实验模拟了6种围压(0MPa、10MPa、20MPa、30MPa、40MPa、50MPa),双轴试验模型的参数及其取值见表1。试验终止标准:作用于左右墙体的平均力与峰值力的比值大于某一指定值时模拟即结束,加载模型见图1。

1.2 模型的构建与参数的选取

以5组不同粒径范围的颗粒填充4个墙体所围成的矩形区域,粒径分别为0.60～1.18mm、1.18～2.36mm、2.36～4.75mm、4.75～9.50mm及9.50～13.20mm。由于颗粒材料具有不规则的边界,颗粒间存在点-点接触、点-面接触和面-面接触,在离散元模型中把颗粒简化为圆形,因此本实验中颗粒之间采用点-点接触的形式,并在接触点上赋予约束。颗粒模型参数见表2,其中颗粒刚度比为法向与切向刚度的比值,生成的5组粒径试件见图2,颗粒个数分别为2653、663、171、42、16。

2 结果分析与讨论

达到试验终止标准后,通过提取20MPa围压条件下的检测结果可以得到弹性模量、泊松比、强度和裂缝初始应力的变化规律,以及颗粒间的接触力分布和位移矢量分布,由此得到粒径对弹性模量、泊松比、强度的影响,以及粒径、围压与裂缝初始应力的关系。

2.1 弹性模量变化规律分析

循环结束后提取5组试件的弹性模量和泊松比,经拟合分析得出粒径与弹性模量和泊松比的关系曲线,如图3、图4所示。由图3可知,随着粒径范围的增大,弹性模量呈减小的趋势,且4.75～9.50mm粒径对应的强度值与拟合曲线的偏离最大,说明该粒径范围的影响较大。粒径与泊松比的拟合曲线具有较大的波动性,由图4仅可看出4.75～9.50mm粒径试件的泊松比出现突变,这与该粒径对强度的影响是相对应的。

2.2 强度变化规律分析

图5为各粒径范围对强度的影响。由图5可以看出,随着粒径的增大,强度呈减小的趋势,且2.36～4.75mm和4.75～9.50mm粒径对应的强度值与拟合曲线的偏离较大,说明该粒径范围的影响较大。这主要是由于粒径较大的颗粒比表面积较小,颗粒间的接触点较少,因此荷载作用下在较少的接触点处产生错动,降低了强度,与混合料中以4.75mm筛孔作为关键筛孔的表述一致。

2.3 裂缝初始应力变化规律分析

试件在承受荷载的过程中出现第一条裂缝时对应的应力为裂缝初始应力。分析了引起裂缝初始应力的2个因素粒径和围压,如图6、图7所示。

由图6可以得出,粒径的增大使裂缝初始应力减小,其变化趋势与粒径对弹性模量、强度的影响一致。

由图7可以看出,对于粒径为0.60～1.18mm和1.18～2.36mm的两组试件,围压对裂缝初始应力的影响具有一致性,两条曲线基本重叠,并且同一粒径下的试件裂缝初始应力随围压的增大基本呈直线增大;对于粒径为4.75～9.50mm和9.50～13.20mm的试件,围压对裂缝初始应力的影响具有与上述两种粒径试件相似的特点,同一粒径下裂缝初始应力与围压基本呈直线增长趋势。粒径为2.36～4.75mm的试件的曲线波动幅度较大,未表现出一定的规律性,其余4种试件其粒径越小同一围压下的裂缝初始应力越大。

2.4 加载后的接触力和位移矢量分析

加载结束后可以得到试件内部颗粒间的接触力分布和位移矢量图,如图8、图9所示。由于试件受压力的作用颗粒间会出现法向力和切向力,颗粒在荷载作用后会产生位移,包括大小和方向的变化。

图8为5组试件内部颗粒的接触力,其中黑色表示压力,其余为拉力,线形越粗表示接触力越大。由图8可以看出,粒径越小颗粒越多的试件接触力分布越均匀,最大接触力越小,粒径越大的试件与之相反,说明较大粒径颗粒形成试件的骨架,抵抗荷载作用主要依赖较大粒径颗粒,粒径较小颗粒对试件的主要贡献为填充空隙。

图9为颗粒间的位移矢量,其中箭头的指向为颗粒运动方向,箭头的长短表示颗粒运动强弱。由图9可以看出,上半部分的颗粒基本向右下角运动,下半部分的颗粒基本向左上角运动,说明试件受到剪切作用而被破坏;粒径较大的试件最大颗粒位移亦较大,粒径0.60～1.18mm试件的最大位移为0.19mm,粒径9.50～13.20mm试件的最大位移为0.27mm,相差0.08mm,并且由于粒径增大颗粒数量减少,试件的位移矢量分布更加稀疏,印证了荷载作用下试件的变形主要源于大粒径颗粒的运动。

3 结论

(1)随着粒径的增大,试件的弹性模量减小,4.75～9.50mm粒径对弹性模量和泊松比的影响最大。

(2)随着粒径的增大,强度减小,且2.36～4.75mm和4.75～9.50mm粒径对强度的影响较大,与混合料级配设计中以4.75mm筛孔作为关键筛孔的表述一致。

(3)围压对各粒径试件的裂缝初始应力的影响基本一致,即围压增大,裂缝初始应力增大。

(4)较大粒径颗粒,其接触力和位移较大,较大粒径颗粒形成骨架,较小粒径颗粒主要起填充作用,并且试件的变形主要取决于大粒径颗粒的运动。

摘要：利用离散元法数值模拟了5组试件的双轴试验,得到了粒径对弹性模量、泊松比、强度和裂缝初始应力的影响,分析了20MPa围压条件下粒径与弹性模量、泊松比、强度和裂缝初始应力的关系,以及6种围压与裂缝初始应力的关系,给出了5组试件荷载作用后颗粒间的接触力和颗粒位移矢量,从细观角度阐释了荷载作用下颗粒的运动模式。模拟结果表明,随着粒径的增大,弹性模量、强度和裂缝初始应力均呈减小的趋势;随着围压的增大,各粒径范围试件的裂缝初始应力呈增加的趋势,并且粒径越大的试件最大接触力和最大位移越大。

关键词：道路工程,粘结颗粒材料,双轴试验,离散元方法,数值模拟

参考文献

[1] Zhou Zhong(周忠),Qin Changguo(秦长国),Xu Yongfu(徐永福),et al.Numerical simulation of roller compactionprocess for slag embankment(宕渣路堤碾压过程的数值模拟)[J].J China Foreign Highway(中外公路),2008,28(6):45

[2] Jiang Mingjing(蒋明镜),Li Xiumei(李秀梅),Sun Yugang(孙渝刚),et al.Discrete element simulation of biaxial com-pression test considering rolling resistance(考虑颗粒抗转动的砂土双轴试验离散元模拟)[J].Rock Soil Mechan(岩土力学),2009,30(S2):514

[3] Wang Guixuan(王桂萱),Qin Jianmin(秦建敏).Numericalsimulation of mechanical behaviors of ellipsoid granular ma-terials by discrete element method(颗粒材料力学性质的离散元数值模拟)[J].J Dalian University(大连大学学报),2007,28(6):514

[4] Zhang Chong(张翀),Shu Ganping(舒赣平).Effect of parti-cle shape on biaxial tests simulated by particle flow code(颗粒形状对颗粒流模拟双轴压缩试验的影响研究)[J].Chi-nese J Geotechn Eng(岩土工程学报),2009,31(8):1281

[5] Zhou Jian(周健),Chi Yuwei(池毓蔚),Chi Yong(池永),etal.Simulation of biaxial test on sand by particle flow code(砂土双轴试验的颗粒流模拟)[J].Chinese J Geotechn Eng(岩土工程学报),2000,22(6):701

数值试验 第8篇

自20 世纪90 年代以来, 锚杆支护以其显著的技术和经济优越性在煤矿巷道控制中获得广泛应用, 是巷道支护的一场革命。国内外很多专家对锚杆支护理论、设计方法、监测技术进行了大量研究, 取得了很多有益的成果。但是, 在锚杆支护巷道中还存在局部冒顶的现象。在众多的锚杆支护巷道冒顶事故中, 其中有一类是锚固系统的失效。巷道顶板锚固系统如图1 所示[1,2,3]。由图1 可知, 锚固剂受孔壁围岩对其的摩擦阻力F、孔底围岩对其的拉力T, 锚杆尾部对其的拉力t, 锚杆壁对锚固剂的摩擦力f。锚固剂在复杂的受力系统中主要受F和f两个力的作用。正因为两个力的存在, 使巷道顶板锚固系统安全工作。而锚固长度及其完整性是两个力保持良好工作的必要条件。巷道锚杆支护系统是隐蔽的, 施工后锚固长度及环形锚固剂的完整程度很难直观得出。因此, 在地面土木与岩土工程中, 采用弹性波进行锚杆无损检测方法已被推荐为规范方法[4]。

国内外对于煤矿巷道锚杆质量无损检测进行了大量研究, 取得了一些重要的成果[5]。但是, 针对于巷道顶板锚固系统进行的无损检测有待进一步研究, 本文基于巷道顶板锚固系统进行不同锚固长度情况下锚杆波导特性, 寻求与现场更接近情况下锚固长度的无损检测方法。

1 煤矿巷道顶板锚杆波导理论分析

1. 1 顶板锚杆波动方程

实际巷道顶板所使用锚杆构件的直径 φ ( φ20 -22mm) 一般远小于其长度L ( 2000 - 2400mm) , 即L> > φ, 故振动时锚杆可采用一维杆件的弹性波理论分析。应力反射波法以一维杆件的弹性波理论为基础[6], 假定锚杆的受激振动在弹性限度内、锚杆材料均匀或分段均匀且各向同性, 锚杆受激振动时, 截面保持为平面。根据运动平衡方程、几何方程及物理方程可得到波在锚杆锚固系统中的波动方程为:

式中, t为时间, u为质点纵向位移, A为截面面积, R为锚固介质与锚杆界面之间单位长度上的粘结摩擦力, C为应力波波速, 其与介质密度 ρ 和弹性模量E有关, 可近似表示为[7]

1. 2 树脂锚杆锚固质量检测原理

土木和岩土工程中锚杆质量无损检测的内容为锚杆长度和注浆密实度。在煤矿巷道中, 锚固体由锚杆、树脂锚固剂和围岩组成。当在树脂锚杆中传播的应力波波长 λ > 10d ( d为树脂锚杆直径) 且 λ< < L ( L为锚杆长度) , 可将树脂锚杆简化为嵌入围岩的一维匀质变截面杆件, 锚固长度和树脂锚固剂密实度的变化表现为杆件截面面积的变化, 锚杆长度表现为材质的变化。无论锚杆长度和树脂锚固剂密实度的改变, 均表现为广义波阻抗的变化。

当在锚杆端头作用激振力时, 就会在杆端产生应力波, 应力波沿杆轴向传播, 当杆中截面面积或材料性质发生变化时, 入射波将在该截面上发生反射和透射。其反射波和透射波幅值的大小与截面积和波阻抗相对变化的程度有关。由于反射波携带锚杆体内的信息, 通过对反射波内所含的信息进行分析, 对锚杆的锚固质量进行分析评价。

根据声波反射原理, 由于波阻抗变化, 应力波遇到波阻抗界面时, 既有反射, 又有透射。依据反射波的传播时间计算锚杆长度 ( 杆底反射波) 和锚固段长度 ( 锚固界面反射波) , 计算公式如下:

式中, L1为未锚固段长度, C为应力波波速, △t为初波和锚固界面反射波第一次到达的时间, 式 ( 3) 为锚杆锚固体系应力波时域分析的重要依据[8,9,10]。

2煤矿巷道顶板锚杆锚固长度检测数值试验

2. 1 模拟方案确定

对于煤矿锚杆支护巷道顶板, 锚杆直径相对于岩层尺寸来说, 周围岩层是无限大的。根据弹性力学理论, 由单孔周围的切向应力分布衰减情况可知, 它有一个剧烈影响范围, 一般以超过原岩应力的5% 为界。令此影响半径为Ri, 根据相关理论推导得出[6]:, r为圆孔半径。两孔不受影响合理间距为B ≥ 2Ri。以巷道顶板锚固孔为例, 影响范围。因此, 为避免边界条件的影响, 数值模型将锚固介质的半径范围取为100mm。另外, 在建模过程中, 将围岩的最外层一圈作为无限元边界吸收本应该透射到围岩中的能量, 避免在边界处形成波的反射[11,12]。

分别建立自由状态和端锚、加长锚和全锚情况下计算模型。顶板锚杆长2. 4m, 直径20mm, 钻孔直径28mm, 锚固剂厚度4mm, 围岩半径100mm。

2. 2 模型建立

锚固体系简化为由锚杆、树脂锚固剂和围岩3种介质组成的3 层柱状结构, 为了保证计算收敛, 时间步长必须足够小, 同时还应满足动力荷载的波长小于模型中单元的最小边长, 这样才能保证计算结果具有足够的精度。

由于低应变动力检测中所用激振能量很小, 在模拟过程中锚杆、锚固剂和围岩处于弹性状态, 界面之间没有产生相对滑移, 另外, 振动过程中界面之间并未产生显著滑动, 可以认为变形是连续的, 为了简化计算, 锚杆系统均采用线弹性本构模型, 且假定界面间无相对滑动[13,14]。

自由状态和锚固状态下的锚杆模型如图2、3所示。

2. 3 模型物理力学参数选取及边界条件

参考现场工程经验, 顶板锚固系统锚固力一般为120K, 锚杆、树脂及围岩介质的物理力学参数如表1。

在动力问题中, 模型周围边界条件的选取是一个十分重要的内容, 因为边界上会存在波的反射, 对动力分析的结果产生影响。理论上讲, 分析模型的范围设置的越大, 分析结果就越好, 但较大的模型会导致很大的计算代价。为了减少边界对波的反射的影响, 在动态分析中, 无限元可以给有限单元提供“静态”边界, 传递来自有限单元的能量, 吸收效果是很好的, 在垂直入射的情况下, 反射能量少于2% 。因此, 围岩边界采用无限元单元来模拟。

由于无限元单元要求单元节点编号外侧为大数编号, 通过sweep方式划分网格, sweep方向为径向无限远处, 这样可以避免节点排序问题。把围岩的最外层划分出来一层圆环包围, 将圆环做成无限元单元, 会吸收本应该透射到围岩中的能量, 不会在边界处形成波的反射, 这也是本文中边界条件设定的特色所在。

2. 4 激励输入及采样间隔

常规振源为一瞬态锤击脉冲信号, 脉冲信号的主要指标是振幅A和脉冲宽度 τ ( 时间/ms) 。本文中激振力采用半正弦波近似描述, 如图4 所示, 取以下形式:

在瞬态分析求解中, 时间步长的选择非常重要, 对于应力波传播问题, 时间步长的选取应当以应力波在单元之间传播时能捕捉到波为准则。经过实践证明, 对集中质量有限元, 时间步长应满足[13,14,15,16]

其中, h为单元网格尺寸, Vp为应力波传播速度。

2. 5 单元尺寸

在二维和三维有限元模型中, 研究应力波的传播问题, 除一维模型的频散和截止频率外, 随着模型维数的增加还将引入新的问题。因此, 采用数值方法对应力波传播与其衰减特征进行模拟分析时, 应使单元网格划分足够小, 保证能够精确描述应力波的传播现象。在三维模型中, 单元网格尺寸的基本选取原则是 Δλ /20, 在选取原则下采用有限元模型取代连续介质模型分析动力问题导致的误差通常可忽略不计。

3 模拟结果分析

在锚固状态下, 为了使接触关系平稳的建立且计算更容易收敛, 在建立模型过程中设置一个时间长度为1s的接触分析步。计算完成后, 把模拟过程中所记录的锚杆端头振动的相对速度数据提取出来, 结合CAE软件的查询值和XY数据管理器功能精确地提取首波、反射波中振幅最大的点的时间值, 得到不同情况下锚杆波导特性曲线。

3. 1 自由状态时锚杆激振响应特性

当锚杆长度为2. 4m, 自由状态下锚杆激振响应曲线如图5 所示。由图5 可知, 在自由锚杆的底端有多次十分明显的反射, 波形十分清晰, 反射波与入射波的方向相同, 随着反射次数的增加, 反射波的幅值有降低的现象, 说明随着波的传播, 能量有很小程度的衰减。由CAE软件的查询值和XY数据管理器功能可知, 反射波的第一个波峰与初波的波峰之间的位置相差0. 93 10- 4s, 锚杆长度为2. 4m, 则锚杆杆体的纵波波速为2. 4 2 / ( 0. 93 10- 4) =5161 ( m / s) , 与由式 ( 2) 计算的的结果5189m / s误差相差不大。

3. 2 不同锚固方式时锚杆激振响应特性

3. 2. 1 端锚锚固时锚杆激振响应特性

当锚杆长度为2. 4m, 用1 支Z2330 树脂锚固剂锚固0. 4m时, 树脂锚杆激振响应曲线如图6 所示。

由CAE软件的查询值功能和XY数据管理器功能可知, 初波的波峰位置在1. 00002s处, 当应力波从锚杆传递到初始锚固界面时, 广义波阻抗增大, 反射波的第一个波谷位置在1. 000805s处, 与初波的波峰之间位置相差7. 85 10- 4s, 则未锚固段长度为: 5161 7. 85 10- 4/2 = 2. 0256925 ( m) , 锚固长度为2. 4 - 2. 0256925 = 0. 3743075 ( m) 。模型中实际锚固的长度为0. 4m, 可见两者比较接近。

3. 2. 2 加长锚固时锚杆激振响应特性

当锚杆长度为2. 4m, 用1 支Z2360 树脂锚固剂锚固0. 8m时, 树脂锚杆激振响应曲线如图7 所示。

由CAE软件的查询值功能和XY数据管理器功能可知, 初波的波峰位置在1. 00002s处, 当应力波从锚杆传递到初始锚固界面时, 广义波阻抗增大, 反射波的第一个波谷位置在1. 00066s处, 与初波的波峰之间位置相差6. 4 10- 4s, 则未锚固段长度为: 5161 6. 4 10- 4/2 = 1. 65152 ( m) , 锚固长度为2. 4 - 1. 65152 = 0. 74848 ( m) 。模型中实际锚固的长度为0. 8m, 检测数据与实际数据两者比较接近。

当锚杆长度为2. 4m, 同时用1 支Z2330 和1 支Z2360 树脂锚固剂锚固1. 2m时, 树脂锚杆激振响应曲线如图8 所示。

由CAE软件的查询值功能和XY数据管理器功能可知, 初波的波峰位置在1. 00002s处, 当应力波从锚杆传递到初始锚固界面时, 广义波阻抗增大, 反射波的第一个波谷位置在1. 00051s处, 与初波的波峰之间位置相差4. 9 10- 4s, 则未锚固段长度为:5161 4. 9 10- 4/2 = 1. 264445 ( m) , 锚固长度为2. 4 - 1. 264445 = 1. 135555 ( m) 。模型中实际锚固的长度为1. 2m, 可见两者比较接近。

当锚杆长度为2. 4m, 用2 支Z2360 树脂锚固剂锚固1. 6m时, 锚固树脂锚杆激振响应曲线如图9所示。

由CAE软件的查询值功能和XY数据管理器功能可知, 初波的波峰位置在1. 00002s处, 当应力波从锚杆传递到初始锚固界面时, 广义波阻抗增大, 反射波十分清晰, 第一个波谷位置在1. 00036s处, 与初波的波峰之间位置相差3. 4 10- 4s, 则未锚固段长度为: 5161 3. 4 10- 4/2 = 0. 87737 ( m) , 锚固长度为2. 4 - 0. 87737 = 1. 52263 ( m) 。模型中实际锚固的长度为1. 6m, 可见两者比较接近。

当锚杆长度为2. 4m, 同时用2 支Z2360 和1 支Z2330 树脂锚固剂锚固2m时, 树脂锚杆激振响应曲线如图10 所示。

由CAE软件的查询值功能和XY数据管理器功能可知, 初波的波峰位置在1. 00002s处, 当应力波从锚杆传递到初始锚固界面时, 广义波阻抗增大, 反射波的第一个波谷位置在1. 000185s处, 与初波的波峰之间位置相差1. 65 10- 4s, 则未锚固段长度为: 5161 1. 65 10- 4/2 = 0. 4257825 ( m) , 锚固长度为2. 4 - 0. 4257825 = 1. 9742175 ( m) 。模型中实际锚固的长度为2m, 可见两者比较接近。从图中还可以看出, 波在锚固体中传播时, 速度响应曲线衰减很快, 初始锚固界面处的反射波比较清晰, 之后的反射波并不明显。

3. 2. 3 全锚锚固时的锚杆响应曲线

当锚杆长度为2. 4m, 同时用2 支Z2360 和1 支Z2350 树脂锚固剂锚固2. 2m时, 树脂锚杆激振响应曲线如图11 所示。

由CAE软件的查询值功能和XY数据管理器功能可知, 初波的波峰位置在1. 00002s处, 当应力波从锚杆传递到初始锚固界面时, 广义波阻抗增大, 反射波的第一个波谷位置在1. 000131s处, 与初波的波峰之间位置相差1. 11 10- 4s, 则未锚固段长度为: 5161 1. 11 10- 4/2 = 0. 2864355 ( m) , 锚固长度为2. 4 - 0. 2864355 = 2. 1135645 ( m) 。模型中实际锚固的长度为2. 2m, 可见两者比较接近。从图中还可以看出, 当波在锚固体中传播时, 速度响应曲线衰减很快, 除第一次反射波比较明显外, 之后的反射波不明显。

3. 2. 4 误差对比分析

不同锚固方式时模拟检测与实际锚固长度的误差对比分析如表2 所示。

由表2 可知, 数值模拟检测锚固长度与实际锚固长度虽有一定的误差, 但误差均在5% 左右, 检测长度较实际锚固长度偏小, 在实际锚固长度无损检测过程中, 为使检测的锚固长度更接近实际, 可引入一定的修正系数消除误差。

4 结论

1) 研究表明, 当锚固长度小于2 m时, 锚杆开始锚固位置处广义波阻抗增大, 反射波十分清晰, 与入射波的方向相反。随着锚固长度的增加, 当锚固长度超过2 m时, 波在锚固段传播时, 速度响应曲线衰减很快, 之后的反射波不明显, 不容易辨别。

2) 端锚、加长锚固和全锚锚固方式下检测结果和实际情况误差大约在1% - 6% , 检测长度较实际锚固长度偏小, 在实际锚固长度无损检测过程中, 为使检测的锚固长度更接近实际, 可引入一定的修正系数消除误差。

3) 对不同锚固方式时树脂锚杆锚固长度进行数值试验, 检测与实际锚固长度数据比较接近, 表明采用应力波无损检测方法检验锚杆锚固长度具有较大的可靠性。

摘要：针对煤岩巷道锚杆支护锚固失效引起的冒顶问题, 建立了煤岩巷道顶板锚杆树脂锚固模型, 采用数值模拟的方法对不同锚固方式下树脂锚杆的波导特性进行研究, 得到了不同锚固长度与树脂锚杆激振响应特性曲线。通过对比分析, 研究得出无论何种锚固方式, 检测结果与实际锚固长度的误差大约在1%-6%, 同时可采用修正补偿的方法降低检测误差。研究结果得出:采用应力反射波法检测树脂锚杆的锚固长度是一种行之有效的无损检测方法, 可实现树脂锚杆锚固长度的便携、准确检测提供技术支持。

泥石流启动过程试验与数值模拟研究 第9篇

1 试验模型设计

1. 1 试验模型及材料

太行山区泥石流发生区域地质构造发育,呈现出坡度陡、沟谷比降大,沟源底部呈扇形或漏斗形的特点,其所在区域内易于松散物及洪水汇集。根据现场调查分析发现泥石流引发的灾害坡度大多数在10° ~ 25°范围之内,为了分析地质灾害破坏最严重的启动机理,本文选取坡度为25°。试验模型沟长为200 cm,宽深比为5∶ 4,在主沟断面上布置高为20cm的玻璃板,其目的是拦截水流,保证过往水流平顺。整个模型结构由供水系统、沉淀区域以及沟槽结构三部分组成( 图1) 。为了保证试验的准确性, 取来自太行山地区的土样进行级配试 验 ( 图2 ) 可知,堆积土体的级配呈现不连续性,为不均匀的砾石土 。

1. 2 泥石流启动试验过程与试验结果

按照图2中的颗粒级配配置试验土样,分层堆积并压实,土层厚度为10 cm,配置土体后的试验模型如图3所示,按照0. 72 cm2/ s的初始单宽流量进行试验,试验结果如图4所示。

由图4( a) 可看出径流历时420 s时,坡趾位置细粒颗粒被水流携带流走,随着土体饱和程度的提高,坡趾位置处的水流在土体表面形成径流,当土体界面处的渗透压力和浮脱力超过了渗流断面土体的滑动摩擦力,坡趾的端部出现滑动流失[2],此时由于坡趾位置土体流走,相邻土体的部分将在饱和和非饱和位置处产生滑移面并发生溯源侵蚀。滑移面的产生诱发土体不断流失且逐渐向上部土体发展[6],对泥石流的发生起到了铺垫作用。当径流时间达到720 s时,由于顶部土体的快速填充造成中上部土体发生局部雍高,土体表面出现径流现象,而表层土体在水流的作用下发生塌陷、密实使得含水量不断增大,泥石流启动出现了一个暂停时期[图4( b) ],随着土体“饱水状态”的发展,土体的基质吸力作用不断弱化,水体达到了一个“供过于求”的极限状态,由于周围土体的带动以及沉降和重力作用的积聚诱发了土体内部的应力释放,土体发生崩塌,形成冲沟现象[图4( c) ]。

2 数值试验

2. 1 离散元理论与本构模型

离散单元法是研究散体材料等非连续介质的一种数值计算方法,专门用于模拟固体力学大变形问题,其理论基础是基于牛顿第二定律与力-位移定律对模型进行循环计算。计算循环过程如图5所示。

颗粒离散元法不用事先定义材料的宏观本构关系,而是由颗粒接触点的接触本构模型来模拟材料的宏观力学行为。为了模拟不同的物理问题,接触模型主要包含三部分: 1接触刚度模型; 2滑移模型; 3链接模型。接触刚度模型是离散单元法的最基本的接触模型,通过法向和切向接触刚度和规定了颗粒间接触力与相对位移间遵循线性或非线性关系; 滑移模型通过接触点处的摩擦系数,来确定颗粒接触处的切向抗剪强度,用来判断颗粒是否在切向发生滑动; 链接模型通过设置颗粒间的法向、切向链接强度来约束颗粒在接触点处的相对运动,只有当法向应力或切向应力超过各自的链接强度值时,相应链接产生断裂,颗粒才会产生相对的运动。根据三种接触模型的特性,选择链接模型来模拟泥石流的启动过程,但在颗粒离散单元法中链接模型提供了两种最基本的接触模型: BCM( 颗粒接触模型) 和BPM( 颗粒平行链接模型) 。两种接触模型各自具有其特点: 对于颗粒接触模型当接触范围很小时链接产生并且接触力通过接触链接传递; 然而平行链接模型其链接主要存在于颗粒或块体之间的接触范围以内,通过链接同时传递力和力矩[7]。目前BCM和BPM两种模型的链接现在均有被采用[8,9]。BCM接触链接模型其链接的断裂对材料的宏观刚度的变化影响较小; 相反在BPM模型中,只要平行链接发生断裂,接触刚度会立即降低,这样不仅会影响邻近的颗粒,而且对材料的宏观刚度会产生较大的变化,基于此现利用BPM构建土体颗粒的接触本构模型,模拟泥石流的启动过程[10]。

2. 2 数值模型的细观参数

颗粒流模拟泥石流的启动过程,计算模型需要输入的是碎屑物的细观参数,而室内试验仅能获得碎屑物的宏观指标( 如抗剪强度、模量) ,因此,需要建立材料的宏观与细观参数间的关系。但由于目前还没有成熟的理论能够把岩土介质的细观参数和宏观物理参数直接联系起来,这就要求在细观参数的选取利用“试凑法”,拟定细观参数后不断进行调整,进行反复的数值试验,直至所得的试样的宏观参数与物理试验的结果相近,最后确定最适合的细观参数。在试验中,利用双轴压缩试验来确定摩擦系数fc、法向接触刚度kn、切向接触刚度ks等细观参数。采用PFC2D进行碎屑物的双轴试验模拟,并以室内三轴试验为校核目标,最终获得可用于后续模拟的碎屑物细观参数。模型的颗粒细观参数见表1。

2. 3 泥石流过程模拟

按照室内物理模拟试验平台实比例建立数值模拟模型,即水槽长为2. 0 m,高0. 4 m,试验选取水槽坡度与室内试验一致为25°,采用自重沉降法生成颗粒球,即首先用wall、generate及property等命令分别建立模型区域、在较大范围内生成颗粒球以及赋予颗粒参数,然后施加重力场,颗粒沉积达到平衡状态,最后删除模型外的颗粒,建立的模型及颗粒情况如图6所示。

泥石流的主要的控制作用并不在于颗粒间的作用剪力[10],而是由以下两方面决定: 水流冲刷能量与碎屑物势能之和要大于碎屑物摩擦碰撞损耗的能量; 水流冲刷力能使碎屑物颗粒运动。泥石流的启动从前缘开始( 图7) 这与室内试验结果一致,碎屑体端部对泥石流的启动起到了一个铺垫作用,中部土体才是泥石流启动的内因。在水流的推动下,碎屑堆积物的后缘部分开始向前运动,在中部堆积,当发展到一定高度就不再发生变化,呈整体向前推动,这个高度是由龙头颗粒消耗的能量与水流冲刷及颗粒碰撞的能量相平衡决定的。这一现象在室内模型试验中得到了验证,龙头的端部几乎看不见水流,固相的体积浓度明显高于后续泥石流。

随着后续泥石流流速增大,出现明显的表层冲沟,以下切侵蚀为主,兼有坡面片蚀作用。发展到一定程度后,冲沟两侧碎屑物垮塌,被水流带走形成新的泥石流。通过对泥石流的启动变化过程发现,坡面的侵蚀机理是由表及里发展的,直至颗粒被全部冲出,没有整体破坏的现象,见图8。

在泥石流的运动过程中,泥石流碎屑在坡面上高速运动,颗粒间的接触力的规律会产生很大变化( 图9) 。图9( a) 是泥石流尚未启动时的接触力链,接触力的大小与线条的宽度成正比,可以看出接触力的大小与碎屑的堆积深度成正比,力链的分布有很强的规律性; 图9( b) 是泥石流高速运动时的力链分布,可以看出颗粒与边坡的接触力呈现跳跃性的变化,反映出泥石流碎屑的高速运动过程中,颗粒接触力的运动变化过程,同时在力链的顶部可以明显看到由于颗粒相互碰撞而产生的接触力突然增大的情况,颗粒速度分布如图10所示,颗粒在滑动过程中,土体颗粒高速运动。综上可知,PFC2D模拟洪水冲刷形成泥石流的启动、运动及堆积过程与室内试验结果吻合较好。

3 结论

( 1) 泥石流启动过程经历了坡趾溯源侵蚀、滑移面初步形成、相对平静、坍塌破坏四个阶段,在泥石流临近破坏阶段,土体的基质吸力作用不断弱化,水体达到了一个“供过于求”的极限状态,沉降和重力作用的积聚诱发了土体内部的应力释放,土体发生崩塌,形成冲沟现象,对于滑动面的形成对泥石流的启动过程起到了铺垫作用。

数值试验 第10篇

随着计算机技术的进步和基于有限元的计算流体动力学方法的发展,用有限元软件来实现风洞试验模拟的技术也不断地在发展。基于有限元的数值风洞模拟应用中,解决二维翼型的升阻力虽已用一些CFD方法所数值模拟过,但结果并不理想,与此同时随着有限元法应用领域以及有限元方法的迅速发展,ANSYS通过质量、动量和能量三个守恒性质来计算流体的速度分量、压力以及温度,使得其结果更接近风洞试验结果。

1 ANSYS数值模拟风洞试验

1.1 ANSYS软件简介

ANSYS公司由JohnSwanson博士创立于1970年,总部位于美国宾夕法尼亚州的匹兹堡,ANSYS软件是融合结构、热、流体、电磁、声学于一体的大型通用有限元分析软件,可广泛用于核工业、铁道、石油化工、航空航天、机械制造、能源、交通(汽车)、国防军工、电子、土木工程、造船、生物医学、轻工地矿、水利、日用家电等一般工业及科学研究。在其所有的产品系列和工作平台上均兼容。ANSYS的多物理场耦合功能,允许在同一模型上进行各式各样的耦合计算成本,如:热-结构耦合、磁-结构耦合以及电-磁-流体-热耦合。

2.2 ANSYS流体流动分析功能在风洞数值模拟中的应用

ANSYS中的计算流体力学分析功能FLOTRAN可以用于分析二维及三维流场问题,使用ANSYS中用于计算流体动力学分析的FLUID 141和FLU-ID 142单元,可以进行传热或绝热、层流或紊流、压缩或不可压缩等问题的研究。根据不同的分析对象选择不同的分析的方法。

ANSYS流体分析功能已明确了可解决作用于气动翼型上的升力与阻力,并通过后处理模块可以得到流场的速度、压力分布以及翼型附面层内速度与压力分布,并可以得到这些量的相关分布曲线。图1为翼型流场网格与流场速度分布图。

2 翼型三分力系数曲线求解模块的二次开发

利用ANSYS流体流动分析中的紊流分析与不可压缩分析,ANSYS中的二方程紊流模型可计及在平均流动下的紊流速度波动的影响。假设该流体的密度在流动过程中保持不变,或者当流体压缩时只消耗很少的能量,即该流体是不可压缩的,不可压缩流的温度方程将忽略流体动能的变化和粘性耗散。在此状态下,在利用参数化设计语言APDL(ANSYS Parametric Design Language)和用户界面设计语言UIDL(User Programmable Features)对ANSYS二次开发,实现翼型三分力系数曲线的求解。

通过二次开发,将ANSYS中的三分力系数的求解过程制成宏,并在Toolbar上设置启动菜单,见图3,图4为翼型参数的输入方式菜单图。

3 算 例

3.1 模型建立

现以标模翼型为例,利用三分力系数曲线求解模块来求解三分力系数曲线。此次数值模拟为了与风洞试验的数值作对比,在ANSYS中所建模型与NF3风洞试验中的尺寸完全一致,并对边界条件、网格点、紊流模型等参数的设置使得雷诺数为风洞试验值即6106,风速87.63 m/s。有限元标模翼型如图5。

4.2 结果分析

由风洞试验、ANSYS计算得到的三分力系数如下表1所示,三分力系数曲线如图6所示。

图6可见,ANSYS与风洞试验得到的曲线的形状基本一致,如果用影响因子Δl(%)来表示数值模拟与风洞试验之间的差值:

$\Delta l=\frac{A{\rm N}SYS\text{计}\text{算}\text{值}-\text{风}\text{洞}\text{试}\text{验}\text{值}}{\text{风}\text{洞}\text{试}\text{验}\text{值}}$。

对于CL值,Δlmax为26.64%;对于CD值,Δlmax为24.16%;对于CM值,Δlmax为11.30%。由以往其他基于有限差分法的CFD软件数值模拟与风洞试验数据的差值相比较的经验,以上各分力系数的Δlmax相比其他CFD软件所得分力系数较小,即在与其他数值模拟软件相比,由ANSYS得到的三分力系数更接近试验结果。根据流体动力学可知,如果在ANSYS计算中建立1:1模型,必然可以得到更精确的结果,在此算例中,网格采用了非结构化,附面层内网格也达到了一般的计算要求,收敛检测量以压力值110-8作为参考,这对计算机有着比较高的要求。

图7 为从翼型攻角为4 °时翼型周围的压力与速度分布。

4 结 论

ANSYS的流体流动分析功能使得数值模拟翼型风洞试验的功能更加具体化,数值模拟结果有较高的精确度。与此同时,ANSYS所具备的通用后处理器POST1与时间历程后处理器POST26使分析结果更加直观、生动,有效地展示了流场中速度与压力等量的分布。对ANSYS的二次开发,使得其模拟风洞更加方便。

参考文献

[1]王庆五,左.,胡仁喜,等.ANSYS10.0机械设计高级应用实例.北京:机械工业出版社,2006;94—123,308—344

[2]商跃进.有限元原理与ANSYS应用指南(第二版).北京:清华大学出版社,2006:224—256,321—344

[3]李庆龄.ANSYS中网格划分方法研究.上海电机学院学报,2006;14(3):28—30

[4]高兴军,赵恒华.大型通用有限元分析软件ANSYS简介.辽宁石油化工大学学报,2004;24(3):94—98

[5]吴鹏,曾红,韩迈.基于ANSYS的二次开发技术的实现方法.辽宁工学院学报,2004;24(5):25—29

[6]王瑞,陈海霞,王广峰.ANSYS有限元网格划分浅析.天津工业大学学报,2005;(1):8—11