有限元数值模拟（精选十篇）
有限元数值模拟 篇1
Freeman研究了页岩气藏中Knudsen扩散对渗流的影响[10];V.Shabro等人研究了不同传质运移机理 (滑脱效应, Knudsen扩散, Lagngmuir解吸附) 在整个生产过中的贡献[11];V.Shabro等人介绍了通过新的数值模拟算法预测有机页岩中产量变化, 并分析了机干酪根中扩散, 滑脱流, Knudsen扩散, Lagngmuir解吸附等机理的作用[1];Vivek Swami等人详细研究了Knudsen扩散, 气体滑脱效应, 气体解吸附效应, 有机干酪根中扩散效应对总产量的影响, 得出产量预测时必须考虑上述因素的结论[12—13];Prince N Azom等人结合双重介质模型, 建立了考虑Knudsen扩散和气体滑脱效应的页岩气数值模型[14];Farzam.Javadpour研究考虑了Knudsen扩散和滑脱流的表观渗透率模型[15]。笔者根据前人研究把页岩储层中的扩散分为:有机干酪中的扩散和微纳米孔隙中的扩散[9]。针对有机干酪根中溶解气, 基于Fick扩散定律, 结合页岩储层特点建立干酪根溶解气扩散数学模型[16], 并利用有限元分析方法对数学模型进行求解, 得出有机干酪根内部浓度变换规律和出口边界平均扩散通量变化规律[17], 为计算和分析页岩气产量变化规律奠定了基础。同时, 可以通过数值计算, 反推出理想模型内部气体有效扩散系数。
1 多孔介质中气体扩散模型
在自然界和实际工程应用中, 多孔介质中气体流动与扩散越来越受到人们的广泛关注和重视。多数研究表明大部分多孔介质具有随机性、无序性和自相似性等特征。早在1855年, 菲克提出了:在单位时间内通过垂直于扩散方向的单位截面积的扩散物质流量与该截面处的浓度梯度成正比, 即Fick第一定律。根据Fick扩散第一定律, 则单位面积多孔介质中气体扩散质量通量可表示为[16]:
式 (1) 中Deff为多孔介质中气体有效扩散系数, m2/s;C为浓度, mol/m3。
与气体自由扩散系数相比, 多孔介质中气体有效扩散系数还受到岩石孔隙结构和孔隙度的影响, 所以气体有效扩散系数可表示为
式 (2) 中Φ为多孔介质孔隙度, τ为迂曲度, 即气体在孔隙中迁移的实际距离与直线距离比值的平方;D为扩散系数。
2 干酪根中溶解气扩散模型
页岩储层有机质含量丰富, 且干酪根中含有大量溶解气。随着孔隙中游离气的采出, 气藏内部压力下降, 当孔隙压力低于临界解吸压力时, 有机干酪根表面吸附气发生解吸附使得吸附气浓度和干酪根内部气体产生浓度差, 页岩气将自发的由高浓度向低浓度扩散, 补充孔隙表面由于解吸附作用减少的吸附气, 维持吸附-解吸附动态平衡, 这样形成了干酪根中溶解气扩散效应。
根据页岩储层特点以及页岩气赋存形式建立干酪根中气体扩散二维物理模型如图1所示。假设纳米孔隙被有机干酪根包围, 灰色区域代表干酪根微粒, 红色区域代表干酪根孔隙之间的孔隙, 中间白色区域代表纳米孔隙 (rn) , 黄色小球代表气体分子 (溶解气溶解在红色区域, 吸附气吸附在白色区域和红色区域交界处) , 干酪根半径 (rk) 约为微纳米孔隙半径 (rn) 的五倍。
针对干酪根中溶解气的存在特点和研究目的, 假设:
(1) 单相流 (只考虑甲烷气体) 。
(2) 干酪根中整个传质运移过程, 只考虑气体扩散, 不存在对流的影响。
(3) 干酪根外表面被硅氧化合物等无机物包围, 视为无流动边界。
(4) 干酪根内部为均质模型, 孔隙均匀分布。
(5) 初始时刻, 干酪根内部气体均匀分布, 浓度为3 mol/m3, 内边界浓度为1 mol/m3。
(6) 干酪根内部基质颗粒视为有机质颗粒且认为基质颗粒是最小单元, 内部无气体分子。
(7) 甲烷气体溶解在干酪根内部, 主要分布在干酪根有机基质颗粒之间的孔隙中。
(8) 假设孔壁吸附气浓度和孔隙中游离气浓度相等, 且在整个模拟过程中保持不变。
假设干酪根中溶解气扩散满足上述条件, 根据Fick扩散第二定律建立干酪根中溶解气扩散数学模型, 扩散方程可表示为[1]。
初始状态下, 溶解气均匀分布在有机干酪根孔隙中, 则给定初始条件
假设有机干酪根外部包围有硅氧化合物, 黏土矿物等无机物, 所以认为有机干酪根和外部没有气体交换, 即给定外边界条件为无流动边界[1]。
有机干酪根内部气体不断向纳米孔隙表面扩散, 所以干酪根内边界条件满足[12]。
由于干酪根内边界和基质孔隙相连, 内边界上气体吸附气浓度和基质孔隙中游离气浓度相同, 内边界条件为:
干酪根扩散物理模型中, 干酪根微粒视为最小基质单元, 基质颗粒内不再含有气体分子, 所以干酪根内部基质颗粒和基质颗粒空间之间也不存在气体交换, 均为无流动边界[13]:
式中km为传质传导系数 (m/s) ;C1为干酪根内边界浓度 (mol/m3) ;Ci为给定浓度值 (mol/m3) ;n为法向方向;D为干酪根中气体扩散系数 (m2/s) ;rk为干酪根半径 (m) ;C为某一时刻干酪根中气体的浓度 (kg/m3) .
3 等效理想几何模型建立
为了便于计算, 根据上述干酪根扩散物理模型建立相应二维均质理想几何模型, 如图2 (a) 所示:红色大圆环代表整个干酪根块, 半径为rk=1×10-3m;中间中等大小的圆代表微纳米孔隙, 孔隙半径为rn=0.2×10-3m, 含有游离气, 表面吸附有吸附气;红色圆环中的小圆代表干酪根块中最小干酪根微粒, 假设大小一致, 半径为0.1×10-3m, 且均匀分布;红色区域为溶解气溶解区域以及气体扩散区域。几何模型采用三角形网格剖分, 如图2 (b) 所示。
4 模型求解
几何模型比较复杂, 计算困难, 故采用有限元分析方法, 借助多物理场耦合模拟软件COMSOL Mul-tiphysics 4.3, 对上述数学模型进行数值求解。模型参数如表1所示:初始时刻干酪根基质孔隙中气体均匀分布, Ci=3 mol/m3, 内孔隙边界处吸附气浓度和微纳米孔隙内部游离气浓度相等, C1=1 mol/m3 (图3.A所示) , 基质孔隙气体浓度和干酪根内边界浓度产生浓度差;由于干酪根外表面由硅氧化合物等无机物包围, 无流量供给, 所以干酪根外边界设置为无流动边界。
模型运行1 s后, 不同时刻干酪根内部气体浓度分布如图3所示, 随着扩散的进行, 浓度梯度由内边界向外边界传播[如图3 (a) ~ (d) 所示], 直到干酪根基质内部气体平均浓度接近内边界浓度[如图3 (d) 所示]。气体扩散过程中, 由于干酪根基质颗粒的不均匀分布以及孔隙连通性影响, 基质内部气体分布不均匀, 如图3 (b) ~ (d) 所示, 上半部分连通性要明显好于下半部分直接导致下半部分扩散速度较慢, 剩余气体浓度分布大于模型上半部分。
不同点处浓度变化曲线如图4所示:当t<0.5 s时, 不同点浓度分布差异明显, 越靠近内边界浓度越大, 但是随着扩散的进行浓度差异不断缩小, 当t>0.5 s时, 各点浓度差异不明显, 当模拟到1 s时, 模型内各点处浓度基本趋于内边界浓度。
内边界平均扩散通量变化曲线如图5所示:当t<0.1 s时, 扩散通量迅速下降, 当0.1 s<t<1 s时, 下降速度变缓, 并逐渐趋于0, 这是由于干酪根内边界浓度 (C1=1 mol/m3) 始终等于孔壁吸附气浓度等于孔隙中游离气浓度且始终保持不变, 然而随着扩散的进行干酪根内部溶解气浓度不断降低, 使得孔壁内外浓度梯度迅速降低, 所以算扩散通量时, 表现出如图5所示的变化规律。
5 多孔介质中气体有效扩散系数的计算
扩散系数是衡量天然气扩散能力大小的重要物理量之一, 由于气体扩散受到孔隙结构和孔隙度影响, 多孔介质中气体有效扩散系数低于自由扩散系数。国内外学者在天然气扩散系数测定方面做了大量探索性工作, 不同学者测定岩石中扩散系数的实验方法不同, 一般情况下, 不同实验方法 (常用实验方法:游离烃浓度法, 水溶烃浓度法, 时滞法) 测量结果差异很大[16]。本文通过修正上述有机干酪根溶解气扩散模型, 利用数值方法计算理想几何模型中气体有效扩散系数。其他条件不变, 假设外边界为供给边界, 可表示为9式 (Cinlet=3 mol/m3) , 内边界为流通边界, 传质传导系数km=5 m/s。
基于有限分析元法求解式 (3) 、式 (6) 、式 (8) 、式 (9) 构成的扩散模型, 由于外边界 (入口边界) 为供给边界, 当扩散进行到一定阶段时整个模型内部浓度分布达到平衡状态并在不改变条件的情况下一直保持下去。模型运行1 s后, 不同时刻各点处浓度分布如图6所示:当0<t<0.1 s时, 内边界附近浓度逐渐降低[如图6 (a) ~ (c) 所示], 当t>0.1 s时, 干酪根内部气体分布逐步达到稳定, 各点处浓度保持动态平衡。x轴正半轴上不同点处浓度变化曲线由图7所示, 内边界 (出口边界) 处浓度为0 mol/m3, 外边界 (入口边界) 处浓度为3 mol/m3, 其他各点浓度由3 mol/m3开始下降, 扩散0.1 s后各点浓度值保持恒定。
内边界处扩散通量变化曲线如图8所示:开始时 (t<0.1 s) 由于内边界附近浓度梯度较大, 所以扩散明显扩散通量较大, 但由于内边界附近浓度下降, 扩散通量迅速下降并最终保持稳定。当气体扩散达到平衡状态时 (平均扩散通量一定, 各点处浓度保持动态平衡) , 可以根据扩散通量定义计算该条件下模型有效扩散系数。
出口端平均扩散通量可表示为式 (10) [17]。由图7可知, 扩散0.1 s以后, 模型内浓度基本达到平衡, 出口端平均扩散通量保持恒定Nave=6.931×10-3mol/ (m2·s) 。
有效扩散系数和平均扩散通量关系可表示为:
出口边界平均浓度表示为
通过式 (12) 计算出口边界处 (内边界) 平均浓度值Cout=2.3×10-3mol/m3, 所以由式 (11) 可计算有效扩散系数为:
式中Nave为出口边界平均通量, mol/ (m2·s) ;Cout为出口边界平均浓度;L0为模型内边界长度 (m) ;L1为模型外边界的长度 (m) ;Cin为外边界 (入口边界) 浓度 (mol/m3) 。
通过数值计算得出该几何模型中气体有效扩散系数 (1.8×10-6m2/s) 小于自由扩散系数 (5×10-6m2/s) 。
6 结论
(1) 基于Fick扩散定律, 结合页岩储层气体赋存特点, 建立了页岩储层有机干酪中溶解气扩散数学模型, 并利用有限元分析软件COMSOL Multiphysics对该模型进行数值求解, 得到干酪根内部浓度分布规律和内边界气体平均摩尔通量变化曲线。
(2) 页岩储层干酪根中溶解气扩散模型满足:外边界被无机质包围为供给边界, 内边界和大孔隙相通为流出边界, 气体均匀分布在干酪根基质颗粒之间, 基质颗粒作为最小单元和颗粒之间孔隙无流体交换。
(3) 在原扩散模型基础上, 改变边界条件, 通过数值模拟方法计算多孔介质中气体有效扩散系数。
有限元数值模拟 篇2
航空发动机薄壁机匣激光焊接有限元数值模拟
对航空发动机燃烧室薄壁机匣建立了激光焊接的数值分析模型.基于Sysweld的`焊接分析功能,用数值方法研究了燃烧室薄壁机匣环缝焊接时引起的焊接变形,模拟了薄壁机匣激光焊接时温度场、应力场的分布以及焊接变形情况.同时对焊接仿真中的热源问题从数学上进行了分析和探讨,在此基础上结合整个薄壁机匣的实际工作状况对其仿真结果进行了定性分析,验证了其分析方法的可行性.
作 者：汪苏 王春侠 Wang Su Wang Chunxia 作者单位：北京航空航天大学机械工程及自动化学院刊 名：航空制造技术 ISTIC英文刊名：AERONAUTICAL MANUFACTURING TECHNOLOGY年，卷(期)：“”(3)分类号：V2关键词：薄壁机匣 激光焊接 温度场 应力场变形
有限元数值模拟 篇3
关键词：机车；吊具；强度；有限元
中图分类号：TP391.72
机车在发运吊装过程中，吊具结构是否具有足够的强度以满足吊装的要求非常重要。吊具在吊装机车过程中主要由发运吊具和车体吊销以及吊车用钢丝绳共同作用实现吊装。机车吊具、吊销的强度决定了吊装过程是否安全。本次分析了新西兰机车发运吊具的强度以满足其发运吊装的要求。
1 机车发运吊具概况
本次分析计算的机车发运吊具主要由横向撑杆和纵向撑杆组装而成。横向撑杆由L50×50×6和L36×36×5的角钢和钢板焊接而成；縱向撑杆由L70×70×8和L36×36×5的角钢和钢板焊接而成，材质均为Q235-A。车体吊销的为圆柱型，吊装时安装在车体底架的吊装位置，材质为45号钢。
机车发运吊具和车体吊销应能达到如下所述的安全系数。（1）在108t载荷作用下，机车发运吊具的安全系数不小于1.8；（2）在27t载荷作用下，车体吊销销体的安全系数不小于1.6。
表1列出了机车发运吊具及吊销材料的基本力学特性：
表1 材料的屈服极限和许用应力
材料部件屈服极限MPa许用应力MPa
Q-235机车发运吊具235131（s=1.8）
45号钢吊销355222（s=1.6）
2 吊具及吊销模型
本文选用大型通用的集CAD/CAE/CAM于一体的软件I-DEAS进行吊具和吊销模型的建立以及有限元强度分析。
2.1 吊具模型
由于吊具的横向和纵向撑杆主要由不同型号的角钢和厚板焊接而成，在有限元仿真模拟时主要采用壳单元进行模拟，所以几何建模时对于角钢建立的几何模型是面，厚板建立的几何模型是体。对撑杆分别建立几何模型然后组装在一起，撑杆的几何模型如图1-4所示。
图1 横向撑杆（一）几何模型 图2 横向撑杆（二）几何模型
图3 纵向撑杆（一）几何模型 图4 纵向撑杆（二）几何模型
由于吊具具有纵向及横向两个对称面，并且所受的载荷也是对称的，因此取1/4结构进行有限元分析。采用壳单元和四面体单元进行结构离散，共划分111314个单元和70349个节点，网格划分如图5所示。
图5 1/4机车发运吊具有限元网格
由其结构特点和力学性质可得：当取1/4结构进行有限元分析计算时，约束施加在吊具结构的对称面上。
2.2 吊销模型
车体吊销安装在车体底架相应的吊装位置，用来吊装机车车辆。所以取出底架侧梁吊装的相应位置和吊销进行接触模拟，具体几何模型及有限元模型如图6、7所示：
图6 车体吊销结构简图 图7 车体吊销和底架装配有限元网格
载荷27t作用在钢丝绳与其接触的位置；吊销销体与底架的接触处设置接触约束；在底架侧梁取出的部分边界设置固定约束。
3 应力计算结果分析
3.1 吊具应力分析
经过有限元分析计算，机车发运吊具节点节点的最大Mises应力为129MPa，发生在端板与纵支承槽钢的交接处，小于材料的许用应力131 MPa，节点的Mises应力分布如图8、9示。
图8机车发运吊具节点Mises应力云图 图9 机车发运吊具局部应力最大点Mises应力云图
3.2 吊销应力分析
车体吊销节点最大Mises应力为220MPa，发生在吊销体与底架侧梁接触处，小于材料的许用应力221 MPa，节点的Mises应力分布如图所示。
图10车体吊销节点Mises应力云图
4 结束语
在机车发运吊具设计时，采用计算机仿真模拟的方法进行强度分析，可以减少设计成本和时间提高效率。由于在实际机车发运吊装中，吊销和吊具在检测时有裂纹出现，所以采用有限元强度分析可以及时发现设计中的不足，并方便改进设计方案，避免和减少在机车吊装试验中出现失误的几率。对于已有吊具，当载重增加时强度是否满足要求，能否保证吊装安全，也可以采用这种方法来实现强度的校核。
本文运用I-DEAS有限元分析软件对新西兰机车发运吊具和吊销进行了有限元强度分析，分析结果表明在吊装时其应力均小于许用应力，强度满足要求，性能可靠。可以为以后的机车发运吊具设计提供参考。
参考文献：
[1]GB/T 20118-2006，一般用途钢丝绳[S].
[2]张质文等.起重机设计手册[M].北京：中国铁道工业出版社，1998.
[3]何翠微.机车整体吊装装置的选择[J].内燃机车，1990（12）：13-16.
[4]杨建新等.I-DEAS入门与提高[M].北京：清华大学出版社，2002.
[5]李伟等.NX Master FEM基础教程[M].北京：清华大学出版社，2005.
作者简介：王玉艳（1977-），女，黑龙江鸡西人，讲师，硕士，研究方向：机车车辆结构强度分析。
作者单位：大连交通大学车辆教研室，辽宁大连 116028
有限元数值模拟 篇4
上海地区广泛分布着具有含水率大、压缩性高、强度低等特性的深厚软土，为了满足工程需求必须对其进行加固处理。真空预压技术作为一种经济可靠、快速简便的软土地基排水固结处理技术在上海地区的应用还不多。为了对软基处理效果进行研究，为工程的设计和施工提供依据，本文以上海浦东某软土地基加固工程为算例，利用ABAQUS建立合理的有限元模型，对真空预压加固软基进行了数值分析，介绍有限元分析的过程和结果，材料的本构模型采用修正剑桥模型。
1 工程及地质情况
场地位于浦东川沙镇，加固区面积60m×80m，砂垫层采用厚度为50cm、含泥量小于5%的中粗砂。密封墙采用水泥土搅拌桩，水泥掺入量为5%，搅拌桩桩长10m。塑料排水板正三角形布置，间距为1.4m，插打深度13.5m。真空预压期间，膜下真空度不小于80kPa。各层土的物理力学性质如表1所示。
2 有限元分析过程
2.1 模型的建立
采用二维平面应变模型，此时根据文献[1]对砂井的渗透系数进行调整。对竖向排水通道的模拟，为了避免网格划分过密以及提高模型运算的速度和效率，将砂墙的间距扩大到4m，根据对称性，考虑到计算的方便，取加固区的一半进行研究。地基的计算宽度为120m，其中40m是加固区，80m是影响区。由于加固区塑料排水板的打设深度为13.5m，考虑到砂垫层以及竖向加固影响范围，按土质条件及真空预压的影响深度，计算深度取40m。
2.2 有限元网格及边界条件
左右两边均为不透水边界，且只有竖向位移无水平位移；底面为无位移边界；网格划分：采用平面应变孔压单元类型(CPE8RP)，8结点二次缩减积分四边形网格，结构化网格(structured)划分。加固区域网格划分较密，往远处逐渐变稀疏。
2.3 荷载及分析步的建立
在定义荷载分析步前先定义一个地应力分析荷载步（即Geostatic步），在加固过程中，整个试验段的加载过程，共分为施加真空荷载、卸除真空荷载两个步骤。根据真空预压加固地基的工程实际情况，将真空荷载分为3个子步，进行分级施加，第一步施加80kPa真空荷载，时间步为一天；第二步持续抽真空，时间步为55天；第三步为卸载真空，时间步10天。在模拟时沿着塑料排板的深度将膜下真空度简化为线性负压设置为80kpa，排水板打设底部的真空度设置为20kPa，中间沿线性递减变化。
3 模拟结果分析
3.1 地基沉降分析
从图中的计算结果我们可以看到，加固区中心位置竖向位移最大，软基的沉降成“碗底”形状，加固区内的沉降较大，并且沉降值随着深度的发展逐渐减小，这与理论的沉降情况相吻合。
3.2 侧向位移结果分析
从图4中可以看出在真空预压期间，整个加固区及加固区外的土体都是向内发生位移的，在加固区边缘处的侧向位移达到最大值，并且影响范围比较广，如图所示，加固区外60m范围内均有一定的影响，在采用真空预压法施工时应该考虑对周围环境的影响。
3.3 孔隙水压力结果分析
孔隙水压力是软土地基处理过程中反映土体中应力变化的重要指标，根据比奥固结方程，孔压的消散与土体的变形是密切相关的，孔压不仅能表示土体固结状态，而且能反映土体的稳定状态。图中反映的是负的超静孔隙水压力，塑料排水版插打深度为13.5m，从图中可以看出，20m深度范围内孔压变化明显，随着深度的增加孔压消散值逐渐减小，在塑料排水板打设深度以下仍有负压存在，即真空预压的影响深度超过排水系统的长度，并且在加固影响区也有一定的影响范围。
4 结论
(1)在砂垫层和排水板单元中施加负的真空压力，排水板中真空度沿深度按线性变化，负的真空压力的这种施加方法与实际情况比较吻合。(2)真空预压处理软土地基，将塑料排水板等效成砂墙地基，采用剑桥模型计算时和实测值比较接近，可以达到理想效果。(3)加固后孔压的变化表明土体有效应力得到增加，有效应力的增加与负压的分布有着重要的关系，同时数值模拟还表明负压可传递至塑料排水板深度以下。(4)利用ABAQUS软件对一般的真空预压加固软土地基工程进行有限元分析计算是可行的，计算结果对工程应用有一定帮助。
参考文献
[1]赵维炳,陈永辉,龚友平.平面应变有限元分析中砂井的处理方法[J].水利学报,1998,(6):53-57.
[2]沈珠江,陆瞬英.软土地基真空排水预压的固结变形分析[J].岩土工程学报,1986,8(3):7-15.
[3]李尧臣,周顺华.真空排水固结法处理软土地基的数值模拟[J].上海铁道大学学报,2000,21(10):96-99.
[4]陈兰云,朱建才.真空预压影响区安全措施有限元分析[J].岩石力学与工程学报,2005,24(2):5712-5715.
[5]徐锴,梅国雄,宰金珉.真空预压土体的变形研究[J].岩土力学,2007,28(11):2471-2474.
隧道全断面施工数值模拟 篇5
根据新奥法原理并结合工程地质条件,采用全断面开挖施工,对隧道施工进行了数值模拟,就其模拟结果进行了详细分析,从而找出隧道施工中围岩及支护结构的应力应交发展变化规律.
作 者：何健 魏碧辉 谭张琴 HE Jian WEI Bi-hui TAN Zhang-qin 作者单位：何健,HE Jian(成都理工大学环境与土木工程学院,四川成都,610059)
魏碧辉,WEI Bi-hui(中铁二十一局集团第三工程有限公司,陕西成阳,71)
谭张琴,TAN Zhang-qin(成都理工大学,四川成都,610059)
刊 名：山西建筑 英文刊名：SHANXI ARCHITECTURE 年，卷(期)：2009 35(3) 分类号：U455 关键词：隧道 全断面施工 数值模拟 围岩 支护结构
岩土工程数值模拟方法的发展 篇6
关键词：岩土工程 数值模拟有限差分有限元边界 元离散 元无界元
1.引言
近几十年来，随着计算机应用的发展，数值计算方法在岩土工程问题分析中迅速得到了广泛应用，大大推动了岩土（体）力学的发展。在岩土（体）力学中所用的数值方法主要有以下几种：有限差分法、有限元法、边界元法、加权余量法、半解析元法、刚体元法、非连续变形分析法、离散元法、无界元法和流形元法等。下面就对这些方法进行简要的介绍和分析。
2.有限差分法
有限差分法是一种比较古老且应用较广的一种数值方法。它的基本思想是将待解决问题的基本方程和边界条件近似地用差分方程来表示，这样就把求解微分方程的问题转化为求解代数方程的问题。亦即它将实际的物理过程在时间和空间上离散，分解成有限数量的有限差分量，近似假设这些差分量足够小，以致在差分量的变化范围内物体的性能和物理过程都是均匀的，并且可以用来描述物理现象的定律，只是在差分量之间发生阶跃式变化。有限差分法的原理是将实际连续的物理过程离散化，近似地置换成一连串的阶跃过程，用函数在一些特定点的有限差商代替微商，建立与原微分方程相应的差分方程，从而将微分方程转化为一组代数方程，通常采用“显式”时间步进方法来求解代数方程组。
3.有限单元法
有限元法将连续的求解域离散为有限数量单元的组合体，解析地模拟或逼近求解区域。由于单元能按各种不同的联结方式组合在一起，且单元本身又可有不同的几何形状，所以可以适应各种复杂几何形状的求解域。它的原理是利用每个单元内假设的近似函数来表示求解区域上待求的未知场函数，单元内的近似函数由未知场函数在各个单元节点上的数值以及插值函数表达。这就使未知场函数的节点值成为新未知量，把一个连续的无限自由度问题变成离散的有限自由度问题。只要解出节点未知量，便可以确定单元组合体上的场函数，随着单元数目的增加，近似解收敛于精确解。按所选未知量的类型，有限元法可分为位移型、平衡型和混合型有限元法。位移型有限元法在计算机上更易实现，且易推广到非线性和动力效应等方面，故比其他类型的有限元法应用广泛。
4.边界元法
边界元法出现在20 世纪60 年代，是一种求解边值问题的数值方法。它是以Betti 互等定理为基础，有直接法与间接法两种。直接边界元法是以互等定理为基础建立起来的，而间接边界元法是以叠加原理为基础建立起来的。边界元法原理是把边值问题归结为求解边界积分方程的问题，在边界上划分单元，求边界积分方程的数值解，进而求出区域内任意点的场变量，故又称为边界积分方程法。边界元法只需对边界进行离散和积分，与有限元法相比，具有降低维数、输入数据较简单、计算工作量少、精度高等优点。比较适合于在无限域或半无限域问题的求解，尤其是等效均质围岩地下工程问题。边界元法的基本解本身就有奇异性，可以比较方便地处理所谓奇异性问题，故目前边界元法得到研究人员的青睐。
5.加权余量法
加权余量法也是一种求解微分方程的数值法，它在流体力学、热传导以及化学工程等方面应用较广。它具有两个方面的优点：①由于加权余量法是直接从控制方程出发去求解问题，理论简单，不需要复杂的数学处理，且它的应用与问题的能量泛函是否存在无关，因而它的应用范围较广，利用加权余量法这一优点去建立有限单元的刚度矩阵，可以大大扩展有限元法的应用范围；②加权余量法的计算程序简单，要求解的代数方程组阶数较低，对计算机内存容量要求不高，计算所需要的原始数据较少，这样就大大减轻了准备工作量。
6.半解析元法
半解析元法是Y. K. Cheung 于1968 年提出来的，同有限元法一样，它也是基于变分原理的。不同点是半解析元法根据结构的类型和特点，利用部分已有的解析结果，选择一定的位移函数，使解沿某些方向直接引入已知解析函数系列，而不再离散为数值计算点，因此自由度和计算量大大降低。这几年半解析法发展很快，种类很多，主要包括有限条法、有限层法、有限厚条法、有限壳条法、样条有限元法以及无限元法等。这类方法适用于求解高维、无限域及动力场等较复杂的问题。
7.无界元法
无界元法是P. Bettess 于1977 年提出来的，用于解决用有限元法求解无限域问题时，人们常会遇到的“计算范围和边界条件不易确定”的问题，是有限元法的推广。其基本思想是适当地选取形函数和位移函数，使得当局部坐标趋近于1 时，整体坐标趋于无穷大而位移为零，从而满足计算范围无限大和无限远处位移为零的条件。它与有限元法等数值方法耦合对于解决岩土（体）力学问题也是一种有效方法。上述介绍的几种数值法都是针对连续介质的，只能获得某一荷载或边界条件下的稳定解。
8.离散单元法
离散单元法随着非连续岩石力学的发展而不断进步，与现有的连续介质力学方法相比，还有以下问题需要研究：
（1）刚体离散单元法是基于非连续岩石力学的，更适合于低应力状态下具有明显发育构造面的坚硬岩体的变形失稳分析。对于软弱破碎、节理裂隙非常发育和高应力状态下的岩体变形失稳分析，则不适合。
（2）岩体介质种类繁多，性质非常复杂。在通常情况下，节理岩体或颗粒体表现为非均质和各向异性，并且常表现有很强的非线性，所处的地质环境不尽相同，这就使得岩土工程计算有很多不确定性因素。离散元的主要计算参数（如阻尼参数、刚度系数），影响到岩土工程稳定过程的正确模拟以及最终结果的可靠性，尤其是离散元计算中的参数选取，没有统一和完善的确定方法。
（3）计算时步的确定。现在的选取原则是出于满足数学方程趋于收敛的条件，与实际工程问题中的“时间”概念如何联系起来，合理地考虑时间效应，是今后需要研究的问题。
（4）迭代运算的时间较长。用计算机进行离散元计算时，CPU 占用时间较多，特别是在考虑岩块变形的情况下，模型划分单元数受到限制，对迭代方法需做进一步的改进。
9.刚体节理元法
刚体节理元法是Asai 在1981 年提出的，它是在Cundall 刚体离散元间夹有Goodman 节理单元的组合单元，但此节理单元有一定厚度而使离散元间不能“叠合”。刚体节理元法也可考虑不含节理单元的情况，即所谓的单一三角形刚体元非连续变形分析法，是石根华博士和古德曼教授于1984 年首次提出的一种新型数值分析方法，至1988 年该方法已形成了一种较为完整的数值计算方法体系。非连续变形分析方法以严格遵循经典力学规则为基础，是一种平行于有限元法的数值计算方法。
10.流形元方法
有限元数值模拟 篇7
TRT焊接机壳是在原有铸造机壳的基础上进行结构改进和创新设计后得出的新产品[1]。与原来的铸造机壳相比,焊接机壳具有多方面的优势。但如此巨大的焊接机壳在各种工况中能否安全可靠,以及在使用过程中的强度、刚度能否满足产品设计要求都是设计人员最关心的问题[2]。因此,虽然国内外学者对这种大型结构件的应力及变形做了大量研究,应用各种有限元方法[3]进行了分析,使得对单个结构件的应力及变形分析变得更迅速和精确,但是如何对同类型的大型结构件进行系统的有限元分析,使生产者和使用者全面了解新产品的力学性能,为新产品的制造和使用提供更可靠的数据参考,仍然是一个亟待解决的问题[4]。
本文通过有限元仿真,获得机壳在刚度和强度方面的特性参数,预知机壳焊接中的变形趋势,分析机壳在使用工况的条件下整个机壳的应力及变形情况,找出应力及变形的最大部位,即最危险的部位,并对结构的调整提出改进性意见。
1 计算模型的建立
1.1结构模型
TRT蜗壳结构如图1所示。蜗壳的最大直径为5.12m,最大厚度为70mm。以中分面法兰为界,蜗壳分为上下两个机壳。壳体内部有大量防变形支撑,如加强筋、支撑柱等。
1.2有限元模型
由于整个TRT蜗壳结构复杂,体积庞大,壳体上存在许多小孔和尖角,造成有限元网格剖分困难。因此,在不影响计算精度的前提下,忽略结构中相对尺寸很小的局部区域,例如许多小孔及倒角等,实现对几何模型的简化[5]。
进行TRT蜗壳结构场模拟分析计算时,壳体采用三维弹塑性实体结构模型。因此,选用三维结构实体单元SOLID45。SOLID45是一种由八节点组成的六面体单元,每个节点具有X、Y、Z方向的3个移动自由度。采用智能和自由网格划分技术对壳体进行有限元网格划分[6],建立与结构尺寸完全相同的三维有限元模型,如图2所示。模型单元总数为301 291,节点总数为86 753。
1.3材料性能参数
在模拟计算中需要给出材料的室温力学性能参数(弹性模量、屈服强度等)。本文计算所针对的TRT蜗壳采用同一种材料Q390C,静力分析以及水压试验分析时,其室温性能参数如表1所示。
1.4边界条件
在TRT蜗壳结构分析中,根据实际运行中的约束条件,在图2中,限制A、B两点Y、Z方向的自由度,限制壳体的竖向以及轴向位移。在C点导向键处限制X方向的自由度,以限制壳体的横向位移。
1.5载荷的施加
在静力结构分析中,所施加的载荷包括壳体自重、承缸重力、静叶调节结构重力以及转子的重力。壳体自重以惯性力的形式加到壳体上,其他构件重力都以壳体配重的形式转化为节点力后,施加到壳体的有限元模型上。
水压试验中,实际测试过程需要在进口和出口处加上挡板,因此在有限元模型上做相应的修改。载荷除了壳体自重外,还包括0.4MPa的静水压力及水的重力。水的重力以节点力的形式加到壳体上。
2 弹塑性有限元模型的计算
TRT蜗壳式焊接机壳实际受力状态非常复杂,其主要载荷包括边界约束力、上下机壳间的摩擦力、自身重力、其他构件重力、内腔压力、热应力等。有限元分析中,机壳总的载荷列阵可表示为
其中,FF为作用于节点的集中力(包括其他构件自重作用力);FV,e、Fs,e、Fσ0,e、Fε0,e分别为作用于单元e的体积力FV(包括自身重力、热应力)、边界分布力Fs(包括边界约束力、上下机壳间的摩擦力、内腔压力)、单元内应力σ0和单元初应变ε0,它们分别为
FV,e=∫VeNTFVdV (2)
Fs,e=∫SeσNTFsdS (3)
Fσ0,e=-∫VeBTσ0dV (4)
Fε0,e=∫VeBTε0dV (5)
式中,N、B、S
机壳的总刚度矩阵K可表示为
式中,D为弹性系数矩阵。
机壳总体有限元方程为
Ka=F (7)
求解式(7)即可获得机壳的综合变形a。通过综合变形a得到单元节点e的位移列向量ae,则单元的应变和应力可通过下式求出:
ε=Bae (8)
σ=D(ε-ε0)+σ0=DBae-Dε0+σ0 (9)
3 计算结果及分析
3.1静力结构分析的计算结果
图3～图7给出了整个TRT蜗壳在静态载荷作用下的应力及变形。从图3可以看出,在室温静力工况下,整体机壳的应力不大,平均在1.2MPa左右,最大等效应力为19.1MPa,远小于材料的屈服强度420MPa,位置如图示3。等效应力主要集中分布在法兰的约束处,这是约束和应力集中的综合效果。
从图4可以看出,机壳的最大合位移为77.6μm,主要位于下机壳支撑柱处,满足变形量小于0.1mm的设计要求。
焊接机壳不同方向位移的最大值如表2所示,各向相对位移值均小于1mm,如图5～图7所示。通过分析可知,在室温静力作用下结构尺寸未发生显著的变化。
μm
3.2整体机壳水压试验计算结果
利用有限元的方法,对TRT蜗壳在承受一定水压作用时,机壳的应力及变形分布情况进行计算,得到的应力和变形如图8～图12所示。从图8可以看出,在水压试验阶段,机壳上的应力大部分布在40～85MPa之间,最大值产生在上机壳中,范围很小,最大值达382MPa。具体位置出现在焊缝处的局部尖角处,属于应力集中。
由以上分析结果可知,水压试验时,机壳应力远未达到材料屈服应力(420MPa),所以机壳变形可以认为是弹性变形。有部分区域应力较大,但分布区域很小,对结构无大的影响。
机壳在生产实际运行过程中,并不加挡板。而水压试验中,整体机壳的最大变形量为1.511mm,最大值出现在挡板上,如图9所示。而机壳的设计要求壳体的变形量不大于5mm。整体机壳受水压作用时,机壳排气腔端板处变形及进气腔的侧板变形较大,且进气腔的侧板有外凸变形,这主要与机壳内所受水压大小不同及侧板的最大变形处的强度有关。X、Y、Z方向的相对变形如图10～图12所示。整体焊接机壳不同方向位移见表3。
mm
可以看出,水压试验中,在水压的作用下,整个壳体的最大合位移是1.511mm,满足变形量小于5mm的设计要求。最大位移是Z向(轴向)的位移,为1.338mm,出现在挡板处。这是由水压所造成的,而焊接机壳本身变形不大。
4 结论
(1)整个焊接机壳受到静力载荷的作用分析结果显示,整体机壳的应力均不大,平均在1.2MPa左右,最大等效应力为19.1MPa,远小于材料的屈服强度420MPa,等效应力最大值主要集中在法兰的约束处。整体机壳的最大合位移为77.6μm,位于下机壳支撑柱处。通过分析得出,整个机壳在室温静力作用下,结构强度和刚度均能满足要求。
(2)通过对整体机壳进行水压试验,得出整体机壳的应力均处于85MPa以下,机壳的应力远未达到材料的屈服强度420MPa,局部应力较大,但分布区域很小,属于应力集中。因此,机壳在水压试验中的强度足够。在水压作用下,壳体的整体合位移最大值为1.511mm,最大位移出现在出口挡板处,机壳上的变形并不大,满足变形量小于5mm的设计要求。因此,机壳在水压试验中的刚度足够。
(3)建议在水压试验中,采用在内部加柱子的方式,进行变形的约束。
参考文献
[1]西安交通大学透平压缩机教研室.离心压缩机原理[M].北京:机械工业出版社,1990.
[2]李雀屏,郑涛,韩清凯.离心式压缩机焊接机壳结构的研究与优化设计[J].装备制造技术,2008(12):78-80.
[3]王瑁成,邵敏.有限单元法基本原理和数值方法[M].北京:清华大学出版社,1997.
[4]孙长辉,刘正先,王斗,等.蜗壳变型线改进离心风机性能的研究[J].流体机械,2007,35(4):1-5.
[5]杜平安.有限元划分的基本原则[J].机械设计与制造,2000(1):34-36.
[6]何磊,田爱梅.低比转速离心泵叶轮造型及有限元网格划分一体化分析[J].流体机械,2006,34(1):27-31.
[7]Zhu X K,Chao Y J.Effect of Temperature-de-pendent Material Properties on Welding Simulation[J].Computers and Structures,2002,80(11):967-976.
有限元数值模拟 篇8
贮液容器在运输和使用过程中,常因受到外部振动、碰撞冲击、跌落冲击而导致破损。容器一旦发生破损就会导致储装的液体泄漏,还可能由于液体具有污染性、腐蚀性、易燃等特点,而导致环境污染,甚至引发人身安全事故。因此,贮液容器在冲击条件下的动态响应问题吸引了国内外学者的兴趣。
梁波[1]针对底部刚性约束支持的薄壁圆柱型贮箱进行了动力特性分析;陈建平等[2]讨论了液体Laplace方程与弹性体Navier方程之间的对应关系,构造了液体晃动的结构弹性体模型;万水等[3]用半解析有限元分析了弹性圆柱贮液器液-固耦合系统的模态特性;李彦民等[4]采用交替解法对地震条件下贮液容器的动态响应问题进行了实例分析和计算;李政等[5]采用ALE有限元法实现了大型车载柔性贮液结构三维动态仿真;孙利民等[6]采用有限元分析软件对制动状态下铁路车载容器的动态响应问题进行了实例分析和计算;Reed等[7]提出当贮液容器长度未超过临界长度时,可以简化为质量弹簧系统进行贮液容器的FEM跌落仿真;Marco等[8]分别采用La-grangian、Eulerian、ALE、SPH 4种算法模拟飞机底板金属水箱跌落碰撞地面的过程,对算法的精度、速度进行了综合比较。综上所述,贮液容器在冲击条件下的动态响应问题的研究主要集中在金属贮液容器的模态分析和容器内液体的动态特性分析方面,且通常将金属容器作刚性化处理,认为容器不存在弹性变形;或是将容器底面视为刚性,侧壁为弹性,整个贮液容器只作微幅运动,而对于工业生产、实验和日常生活中常用的非金属容器动力学分析研究工作尚少。由于非金属材料弹性模量小、变形大,在进行动力学研究时,容器作刚性化处理是不合适的,且还须考虑容器变形与液体间的耦合作用。
本文采用有限元程序ANSYS/LS-DYNA,对PC(聚碳酸酯)饮水桶在自由落体情况下的跌落冲击过程进行数值模拟、分析。
1 贮液容器有限元建模
1.1 问题描述
考虑一圆柱式贮液容器自由跌落过程,如图1所示。其中,L1、L2、L3、r、r0均为容器的外形结构尺寸;b为壁厚;L为贮液高度,它可由容器贮液量和容积结构尺寸计算得出(通常贮液量是用容器的容积百分比表示);h为贮液容器跌落高度。在运输和使用过程中,贮液容器通常以自由落体形式跌落,与地面发生碰撞冲击。一般来说,容器最大的受力和变形出现在容器与地面接触碰撞之时,这也是跌落冲击动力学模拟的重点。
当跌落高度h=0时,即贮液容器静止于地面上,贮液容器受到液体的静压力,容器上最大应力出现在桶底处。由于b≪r,应用材料力学理论,桶体上最大应力可以表示为
undefined
式中,ρ为液体密度;g为重力加速度。
一旦跌落冲击发生,贮液容器上的应力大小与分布就会发生急剧变化。容器结构与液体之间发生耦合作用,即容器结构在液体载荷作用下发生变形或运动,而变形或运动反过来又影响液体的流动。本文主要关注桶体的变形,可将流场部分进行简化,假设水为无黏性、不可压缩的理想流体。
1.2 液固耦合的有限元算法
在流体动力学模拟中,通常采用拉格朗日法和欧拉法两种方法描述流体的运动[9]。拉格朗日法着眼于流体质点,流体质点与其物质坐标始终是一致的。欧拉法基于场的概念,着眼于其上的空间坐标点来观察流经该点的物质质点速度和其他物理量。在贮液容器跌落冲击的数值模拟分析中,为了方便处理碰撞地面时引起的桶体大变形,采用任意拉格朗日-欧拉(ALE)算法。ALE描述方法既不着眼于流体质点,也不着眼于空间点,而是在一个可任意运动的参考点上(即网格点)观察流体运动。
在ANSYS/LS_DYNA程序中的液固耦合算法可以采用以下三种方法:合并流体和结构的界面节点算法、接触算法、欧拉/拉格朗日耦合算法。在欧拉/拉格朗日耦合算法中,结构与流体的几何型以及网格可以重叠在一起,计算中通过一定的约束方法将结构与流体耦合在一起,以实现力学参量的传递。本算例在建立几何模型和有限元网格划分时,结构与流体的几何型以及网格是相互独立的,采用欧拉/拉格朗日耦合算法,对容器结构进行约束,将容器结构的相关参量传递给流体单元,实现真正意义上的液固耦合算法。
1.3 饮水桶有限元模型
本算例中,桶壁厚度b=1.2mm,桶底半径r=137mm,桶顶半径r0=28.5mm,桶体总长L3=486mm。水、空气采用空材料模型和Gruneisen状态方程,桶体采用Bilinear Isotropic模型,地面材料采用Rigid Material模型。由于要分析桶体的变形,因而建模时,必须在水和空气的周围建立一定范围的网格,使得空材料有流动空间。图2所示为饮水桶的有限元模型,整个有限元模型单元总数为40 980,节点总数为44 003。
2 饮水桶跌落冲击数值模拟分析
2.1 数值模拟结果
饮水桶跌落冲击变形和破损一般发生在与地面碰撞接触后,在跌落前的自由落体过程中受力和变形不会变化。所以,在跌落模拟中,为了缩短计算时间,这里假设饮水桶以某一初速度在距离地面1mm处开始跌落。当模拟跌落高度为h时,跌落初速度由undefined确定。
图3所示为饮水桶从1.2m处正跌落姿态时,饮水桶上A、B、C、D点(图2)的应力变化时程曲线。由图3可知:①桶体的等效应力是按照A-B-C-D的顺序出现,即由接触面开始由下而上传递;②桶体各点的应力均呈现由小→大→小的变化趋势,经历时间大约1ms左右,由于数值模拟中没有考虑摩擦等能量耗损,模拟应力并不为0;③桶体各点的应力大小不但与接触面距离有关,还与桶体结构形状相关。从各点最大应力来看,A点在接触面上且又是结构变化处,应力最大;B点离接触面虽然比C点离接触面近,但C点应力比B点应力大,因为C点是结构变化处。
2.2 跌落角度对应力的影响
选取饮水桶口朝左,轴线与地面夹角为0开始,顺时针旋转,每旋转10°进行一次试验,由于饮水桶结构具有对称性,因而旋转180°,进行19次试验,观察其接触瞬时最大等效应力分布情况。
图4~图7分别为轴线与地面夹角为10°、60°、90°、120°时桶体接触地面瞬时的等效应力云图;图8显示了容器跌落角度(姿态)与容器壁受跌落冲击最大等效应力的关系。可以看出:饮水桶跌落角度不同,跌落时冲击地面的接触点也不同;虽然其等效应力均是出现在饮水桶与地面的接触点上,但不同跌落姿态时,桶体受到的等效应力差异较大。
由图8可知:当轴线与地面夹角为120°、130°、140°跌落时,饮水桶与地面接触点均在桶口处,其最大等效应力要远远大于其他情况。饮水桶正跌落90°时,最大等效应力为23.03MPa;与地面夹角140°时,最大等效应力为55.39MPa,是正跌落的2.4倍。由此可知:①在饮水桶设计和耐撞性能试验中,必须充分考虑到跌落角度的多样性和不确定性,不能仅仅考虑某种角度(如正跌落)时的耐撞性能,否则跌落破损就可能时有发生;②桶体底部、顶部、桶口过度处等可能成为跌落接触点,并且容易引起局部应力集中的地方,在结构设计时可采取局部增加厚度、增加过度圆弧半径等措施来增加容器抗跌落冲击耐撞性能。
2.3 跌落高度对应力的影响
以饮水桶壁厚b=1.2mm,贮水量(本文特指桶中水量占整桶水量的百分比)为80%,跌落姿态为底面朝下,与地面夹角为60°,改变跌落高度为0~1.5m,分别进行18次仿真,观察桶体接触地面时等效应力的变化,得到跌落高度与饮水桶跌落地面时的最大等效应力关系曲线如图9所示。
从图9可以看出,随着跌落高度的增大,饮水桶撞击地面时速度增大,饮水桶受到的跌落冲击最大等效应力增大。在正跌落姿态下,取桶内液面高度L=310mm,桶底半径r=137mm,桶壁厚度b=1.2mm,依据式(1)计算求得静态时桶体应力为0.347MPa。与图3中从1.2m处正跌落的应力比较可知,A点等效应力最大值为静态时同一点应力的66.3倍,这时静压力几乎可以忽略不计。这表明饮水桶的结构设计是不能依据液体静压力来进行强度校核的,而必须考虑运输、使用过程中不可避免的跌落冲击破坏等因素。而且,由于跌落高度是不确定的,仅按某一确定高度进行设计和耐撞性试验就会出现要么不可靠,要么太保守的情况。
2.4 容器壁厚对应力的影响
跌落高度h=1.2m,贮水量为80%,跌落姿态为底面朝下,与地面夹角为60°,分别在容器壁厚度为0.4~1.5mm时进行12次跌落仿真,观察桶体接触地面时最大等效应力的变化。图10所示为不同容器壁厚与饮水桶跌落地面时的最大等效应力关系。
从图10可以看出:随着容器壁厚的增大,桶体的跌落冲击等效应力减小,但桶体跌落冲击最大等效应力随壁厚的变化幅度较小,如桶体壁厚由0.4mm增加到0.6mm、1.2mm,其跌落冲击最大等效应力只相应减小了0.3%、4.97%,几乎可以忽略不计;当容器壁厚从0.6mm变化到1.5mm时,桶体跌落冲击最大等效应力直线下降,两者基本呈线性关系,通过回归分析得到两者之间关系为
σ=42-2.64b (2)
由于桶体质量与桶壁厚是正比关系,也即桶体材料成本随重量增大而正比增大,故在容器设计时,必须慎重确定容器壁厚,而不要盲目增加壁厚来提高容器的抗跌落耐撞性。
2.5 装水量对应力的影响
饮水桶壁厚b=1.2mm,跌落姿态为底面朝下,与地面夹角为30°,分别改变贮水量0、20%、40%、60%、80%、100%进行6次仿真,得到桶体接触地面时最大等效应力变化如图11所示。
由图11可知,接触点的最大等效应力随着贮水量的增大而增大。这是因为随着容器内的水量增加,容器总质量增大,冲击力增大,从而使得容器撞击地面时受到的冲击应力必然增大。当贮水量由40%增大到80%,桶体最大等效应力约增大10%,盛满液体的桶体受到的应力比空桶体增大了大约1/3。所以,在容器结构设计和抗跌落耐撞性试验中必须考虑容器贮水量的多少。
3 结论
(1)容器跌落角度(姿态)与容器跌落冲击最大等效应力呈现强差异性。在设计中可考虑在转角等直接与地面相接触点增加过度圆弧半径和容器局部壁厚,以便增强容器抗跌落耐撞性能。
(2)从容器壁厚的变化来看,在一定跌落冲击参数前提下,容器跌落冲击等效应力随壁厚增大而减小。但必须根据强度要求,谨慎设计确定桶体壁厚,而不是盲目增加壁厚来提高耐撞性。
(3)由于跌落冲击参数的随机不确定性对容器抗跌落耐撞性能有显著影响,因而在容器设计和耐撞性试验中确定某一个单一状态指标可能是不合适的。另外,贮液量对容器的等效应力影响也很大,在容器的抗跌落耐撞性设计中应加以考虑。
参考文献
[1]梁波.弹性贮液圆柱壳的动力特性分析[J].国防科技大学学报,1990,12(2):30-35.
[2]陈建平,周儒荣,韦忠?.用结构有限元程序进行液体晃动动力分析[J].中国机械工程,2003,14(19):1631-1633.
[3]万水,朱德懋.圆柱贮液器固液耦合模态分析[J].工程力学,2007,17(3):87-92.
[4]李彦民,徐刚,任文敏,等.储液容器流固耦合动力响应分析计算[J].工程力学,2002,19(4):29-32.
[5]李政,金先龙,申杰,等.车载柔性储液结构动态仿真方法[J].上海交通大学学报,2007,41(1):19-22.
[6]孙利民,赵勇,张庆华.列车制动状态储油罐瞬态响应分析[J].郑州大学学报,2006,27(2):32-35.
[7]Reed P E,Breedveld G,Lim B C.Simulation of theDrop Impact Test for Moulded Thermoplasticcon-tainers[J].International Journal of Impact Engi-neering,2000,24:133-153.
[8]Marco A,Luigi M L.Castelletti,Maurizio Tirelli.Fluid-structure Interaction of Water Filled Tanksduring the Impact with the Ground[J].Internation-al Journal of Impact Engineering,2005,31:235-254.
有限元数值模拟 篇9
数值模拟技术相对类比法、实验室模拟以及现场观察法具有适用性和通用性强的特点, 在防治煤与瓦斯突出中也有应用。配合数值模拟和现场测试, 有限元可以成为研究煤与瓦斯突出规律的有利工具。本文采用数值模拟理论分析与现场试验论证相结合的方法探讨有限元数值模拟在巷道掘进中的应用。
1现场实验地点概况
综掘开挖工艺直接关系到巷道支护和采掘速度, 决定巷道前方集中应力带的形状和范围, 对防治煤与瓦斯突出具有重要意义。现场试验地点选择在平煤集团八矿己15煤层中的己15-12130采面。平顶山煤业集团八矿是个高瓦斯矿, 己15-12130采面位于己二下延采区下部, 上部为已开采的己15-12110采面, 下部为未开采的二水平己二采区, 东临己二下延回风下山, 西临己四采区边界, 标高为-490～-538m, 垂深为574～622m, 方位126°39′, 煤层厚度3.0～3.5m之间, 平均厚3.2m, 倾角12°, 整个走势为东高西低。
己15-12130采面走向长680m, 采长171m, 采高3.2m。煤炭储量46.3万吨。煤层顶板为深灰色砂质泥岩与浅灰色细砂岩, 粉砂岩互层, 厚14.2m, 节理发育, 底部夹0.05～0.15的煤线;底板为砂质泥岩, 厚5～8.5米。己15-12110机巷距己15-12130风巷中心20m, 在己15-12110机巷掘进过程中出现一落差0.8m, 倾角52°的正断层, 己15-12130机巷外段出现一落差0.8m, 倾角32°的正断层。这两条断层对整个采面影响不大, 己15-12130风、机巷掘进过程中不会出现大的地质构造。己15-12130采面开采平面图如图1所示。
2数值模拟方案设计
数值模拟一般分为以下几个步骤: (1) 建立模型; (2) 确定模型受力条件; (3) 模拟; (4) 根据模拟结果进行应力、位移等分析; (5) 根据现场结果验证模拟。
2.1计算模型
完全准确地模拟各种岩体的变形既不可能也不必要。数值模拟方法作为一项重要研究手段, 在人为考虑到各种因素的基础上, 简化模型, 抓住重点, 可以取得事半功倍的效果。在本次有限元数值模拟分析过程中作如下假设和说明:
(1) 由于主要研究综掘开挖过程中顶板的应力变化, 同时直接顶相对上覆岩层的尺寸很小, 可以忽略不计。同时将煤层的顶底板按同一介质处理, 煤层则按照另一介质处理。即在计算模型中, 将顶、底板岩体、矿体当作不同材料物性参数考虑, 但反应在模型中二者之间是紧密连接的连续单元, 在结合面处有相同的节点。
(2) 上覆岩层、煤层和底板为均质、各向同性的, 不考虑岩层和煤层的结构面、裂隙和软弱夹层的影响。
(3) 在计算中取开挖煤层的厚度和倾角保持不变, 为计算和分析方便, 煤层取水平, 煤层厚度取3.5～4.0m。
(4) 己15-12130工作面标高为-490～-538m, 垂深为574～622m, 为工程施工稳妥安全起见, 计算中取最不利荷载和最不利开采情况, 取上覆岩层厚度为650m。
计算模型如图2所示, 其中 (a) 图为计算模型, (b) 图为实际巷道模型。
2.2物理模型及几何边界条件和约束条件
在三维坐标系统中, 以底板基点为坐标原点, 基岩底面为XOY平面, 煤层走向为Y正方向, 煤层倾斜方向为X方向, 垂直向上为Z轴正方向。计算模型范围取:
为使开挖部分不影响边界, X方向取20m, Y方向取15m, Z方向取60m。开挖掘进工作面为己15-12130机巷的实际几何尺寸。
模型底面为位移三向固定约束, 模型侧面为水平垂直约束。顶部受岩石自重应力作用, 如图3所示。
2.3网格划分
划分网格是建立有限元模型的一个重要环节, 需要的工作量较大, 所划分的网格形式对计算精度和计算规模将产生直接影响。网格参数包括单元格长度、网格数量和网格疏密程度。确定网格参数中最重要的是单元格长度的确定, 其取值的大小将直接影响网格划分的数量和网格的疏密。单元格长度值越大, 网格数量越少、间隙越疏, 会忽略一些危险点的计算, 不能准确全面地反映应力应变情况;如果单元格长度值过小, 计算点会太密, 计算时间太长, 效率太低, 从而影响有限元分析的结果。[2]在有限元模型基础上生成的网格划分如图4所示。网格划分模型后, 可以得到, 四面体单元为17316个。
2.4计算模型边界受力的确定
根据对平顶山煤业集团八矿地应力场特征分析, 该区主要受自重应力场作用, 模型顶面作用等效垂直压力, 模拟上覆岩层的压力, 取矿石上方岩石平均体重为2.7t/m3。
undefined
式中:σ2—垂直压力, MPa;
ρ—岩体体重, t/m3;
g—重力加速度, N/kg;
h—深度, m。
3掘进现场试验开挖方案设计
开挖方案在考察现场的基础上, 研究机械掘进过程的技术规程以及听取操作人员的要求后, 将开挖分成两部分进行。第一部分, 主要是以扫底和一次开挖完0.5m为主, 第二部分以开挖四块煤体为主, 确定了四套开挖方案。四套开挖方案的开挖顺序如图5所示。
第一套方案从第一块开挖 (从上帮开始) , 然后第二块, 第三块, 最后再挖下帮。即1→2→3→4, 如图5中的 (a) 图所示。第二套方案开挖从第二块开挖, 然后第一块, 第三块, 最后再挖下帮 (第四块) , 即2→1→3→4, 如图5中的 (b) 图所示。第三套开挖方案从第三块开挖, 然后第四块, 第二块, 最后再挖上帮 (第一块) , 即3→4→2→1, 如图5中的 (c) 图所示。第四套开挖方案从第四块开挖, 然后第三块, 第二块, 最后再挖上帮 (第一块) , 即4→3→2→1, 如图5中的 (d) 图所示。
4模拟结果及其分析
开挖方案的第一部分开挖范围小, 对整体结果也没有影响。因此, 为了讨论方便, 将第一部分的开挖情况忽略, 主要讨论第二部分, 即讨论第一至第四步。四套开挖方案的数值模拟分析是相类似的, 这里以第一套开挖方案为例进行讨论。
4.1第一套方案开挖后的垂直应力变化
第一套开挖方案数值模拟垂直应力变化如图6所示。图6中的 (a) ～ (d) 图表示了第一套开挖方案从第一步到第四步开挖后的巷道面上的垂直应力变化, 从图6的 (a) 和 (b) 图可以看出, 第一步到第二步, 巷道周围出现了压应力 (负号) , 在未开挖的区域, 出现煤体从受拉过渡到受压的情况, 此时煤体很容易发生破碎。总体看应力发生转移, 在已经开挖的断面上出现较大的压应力。第三步至第四步, 断面出现以压应力为主, 在靠近两帮和底板附近出现煤体受拉状态。
4.2第一套方案开挖后的水平应力变化
第一套开挖方案数值模拟水平应力变化如图7所示。图7中的 (a) ～ (d) 图表示了第一套开挖方案四步开挖后, 煤体的水平应力变化。从图可以看出, 随着一步步的开挖, 在煤体边脚处, 煤体受到一定的拉应力, 在开挖第四步后, 开挖后的底板受到最大拉应力在3.5MPa左右。总的来看水平应力相对垂直应力值要小。但每开挖一步, 其水平应力值在增加。
4.3第一套方案开挖后的位移变化
第一套开挖方案数值模拟垂直位移和水平位移变化如图8所示。图8中的 (a) 和 (b) 分别表示了在开挖完四步后, 水平位移和垂直位移的大小, 根据模拟发现, 在水平模拟开挖过程中, 煤体的下帮出现位移增大的趋势, 在靠近下帮的边脚、煤体发生较大的位移, 此处更容易出现掉顶、垮塌现象, 因而需要进行有针对性的支护。
在垂直方向上, 煤体的位移随着开挖的顺序进行移动, 先开挖的煤体有较大的位移, 同时在底板出现较小底鼓现象。
4.4模拟结果讨论及评价
通过数值模拟不同开挖顺序, 得出的巷道前方应力峰值的变化规律还是比较清晰的。为讨论方便, 将各方案的集中应力系数变化和峰值波及范围汇总, 见表1。
从防治突出的观点分析, 集中应力的峰值离煤壁越近, 突出的可能越大;应力峰值波及范围大, 形状平缓, 塑性变形较大, 排放瓦斯的范围较大, 突出可能性相对较小。从四套方案的模拟结果看出, 在开挖开始阶段, 产生的最高应力基本接近, 峰值的变化第一套方案和第二套方案有所降低, 第三、第四有所增加, 特别是第三套增加最大, 这样对开挖后续工作不利, 同时存在较大隐患。从煤壁产生的应力变化看, 第一和第三较优。
当开挖顺序由中间向两边开挖时, 峰值变化较大, 特别是第三套方案, 明显高于由上帮到下帮和由下帮到上帮。而从左往右开挖对煤体扰动和破坏相对较小。从上述可以看出, 第一种方案相对要优于其他几种方案。同时从应力集中系数和后续开挖的安全性上分析, 第一种方案为最好的方案。
5现场试验
现场试验要求在掘进机扫完底后, 暂停作业5分钟, 保证掘进头的瓦斯涌出趋于正常。这时瓦斯检查员在固定位置的巷道断面进行四点瓦斯浓度的测定 (如图9所示) , 瓦斯探头记录员开始记录瓦斯浓度, 每隔1分钟记录一次瓦斯浓度。
掘进机按照四套开挖方案开始切割, 掘进机切割深度0.5m, 宽度1.2m, 在切割第一步过程中, 瓦斯检查员进行上述四点处的瓦斯浓度测定, 同时瓦斯探头处的记录人员每隔1分钟记录瓦斯探头的浓度, 同时注明工艺过程。在切割第二步、第三步、第四步时, 重复上述瓦斯浓度的测定工作, 直至掘进机切割工作面完毕。每一步割煤完毕后, 测定瓦斯浓度10分钟。
在四种开挖方案跟踪测定过程中, 根据实际测定数据和风量计算掘进过程中瓦斯涌出量的结果如表2所示。
从表2可以看出, 第一套方案要比其他方案的瓦斯涌出量小得多。
6结论
数值模拟分析技术在巷道、围岩、岩石方面有着广泛的应用, 在防治煤与瓦斯突出方面也可以找到切入点, 对煤矿安全工作是一个有益的启示。选择合理的开挖方式可以有效地减小应力集中, 减少瓦斯涌出量和浓度, 降低风险。
从数值模拟理论分析和现场试验看出, 第一套方案和第四套方案与平时挖掘过程中断面瓦斯浓度 (挖掘0.2m条件下) 基本接近, 而且在挖掘过程中瓦斯涌出量相对较小。而第二套和第三套方案挖掘过程中, 浓度明显增加, 瓦斯涌出量也较大。原因在于第一、四套方案对煤体前方扰动较小, 煤体落煤后涌出较小所致。综合应力和瓦斯突出两方面因素, 第一套挖掘方案造成的瓦斯涌出量和浓度较小, 而且产生应力集中较小, 是最优的开挖方案, 数值模拟理论分析和现场试验论证的结果得到了较好的吻合。
参考文献
[1]石必明, 刘泽功.保护层开采上覆煤层变形特性数值模拟.煤炭学报, 2008, 33 (1) :18.SHI Bi-ming, LIU Ze-gong.Numerical simulation of the upper coal and rock deformation characteristic caused by mining protecting stratum.JOURNAL OF CHINA COAL SOCIETY, 2008, 33 (1) :18.
[2]饶运章, 柴炜, 黄本文.ANSYS数值计算在紫金山金铜矿矿柱稳定性分析中的应用.黄金, 2008, 29 (7) :22.RAO Yun-zhang, CHAI Wei, HUANG Ben-wen.ANSYS numerical calculation of ore pillar stability in Zijinshan gold-copper mine.GOLD, 2008, 29 (7) :22.
[3]朱秀娟.有限元分析网格划分的关键技巧.机械工程与自动化, 2009, 1 (总第152期) :185.ZHU Xiu-juan.Pivotal Skills in Grid Division in FEA.MECHANICAL ENGINEERING&AUTOMATION, 2009, 1 (Series No.152) :185.
[4]卢成, 程五一.巷道应力拱成拱及防冒顶技术分析.中国安全生产科学技术, 2008, 4 (6) .LU Cheng, CHENG Wu-yi.Analysis of stress arch action and preventing roof falling in laneway engineering.Journal OF SAFETY SCIENCE AND TECHNOLOGY, 2008, 4 (6) .
小振幅波有限差分方法数值模拟 篇10
小振幅波波动模型具有简单实用的特点, 同时又是考虑方程非线性, 能量耗散和波浪随机性等高级模型的基础, 因此具有理论和实用价值。重力作用下的水波运动, 通常可以将它视为不可压缩, 无粘, 无旋的理想流体势运动。当波高相对于波长和静水深度很小, 可以用线性化理论处理时叫作小振幅波[1]。国内外学者对小振幅波进行了大量研究并将其运用到工程领域。Airy[2]首先提出小振幅波理论, 将自由表面边界条件略去非线性项, 简化为线性的自由表面条件。G.Dean[3]详细介绍求解线性势流方程组解析解的全过程, 并对驻波, 行进波进行模拟。Stokes[4]首次将摄动法用于非线性波研究, 将有限振幅的非线性波视为无限多个线性波的Fourier组合, 求得有限振幅非线性波的二阶级数解析解。T.Ramadan[5]在线性浅水波基础上研究海啸的产生和传播。Ben-Menahem[6]利用小振幅波理论计算移动源的二维辐射方向图。Hayir[7]将线性浅水波应用到海洋深度对海啸振幅影响的研究。Das[8]构建自由面势波, 研究在忽略和考虑两种不同情况下, 表面张力对深水自由表面的作用。钱铁[9]使用区域变换方法求解综合考虑折射, 反射, 绕射影响下小振幅波数学模型。张双志[10]根据流体力学基本原理, 对无限深水中的小振幅波进行了总结和讨论。沈纪苹[11]分析了Airy线性波和二阶Stokes波对水中的梁所受水动力及非线性振动特征。
基于以上分析, 对线性势流方程组采用CrankNicolson隐式有限差分方法进行离散。对小振幅波进行数值模拟, 数值结果表明数值解很好地吻合解析解。从小振幅波不同时刻的波形发现, 波在波节处静止不动。
1 小振幅波数学模型
本节介绍定义在二维水槽的小振幅波数学模型, 假设水槽中的流体为不可压缩, 无粘, 无旋的理想流体。笛卡尔坐标系 (xoz) 建立在自由面上 (如图1) 。
其中静止液面中点为原点, x轴正方向水平向右, z轴正方向垂直向上, 其中原点到水槽左右两侧边界为b1, b2, 到底部边界为-d, 自由面到水槽底部的距离表示为h=d+η, η表示静止水面到自由面的波高。
下面给出小振幅波数学模型[3]:
其中φ表示速度势, t表示时间, p表示液体静压力, p0表示自由面压力, ρ表示液体密度, g为重力加速度。
速度与速度势的关系
2 区域坐标变换和方程组无量纲化
2.1 坐标变换通过坐标变换公式
将如图1的物理区域转换为如图2的正方形计算区域, 用平行于坐标轴x, z的线段将物理区域分割成 (M+1) × (N+1) 个小区域, 其中
2.2 线性势流方程组无量纲化
下面对速度势, 自由面波高, 压力, 时间, 速度分量进行无量纲化处理, 选用静水深度d和重力速度g为基量纲。
将上述公式和坐标变换公式 (7) 代入方程 (1) - (6) 进行无量纲化, 可得
其中ΦT表示Φ对无量纲量T的导数, 其他导数类似, 系数c1, c2, …, c6表示如下
3 Crank-Nicolson有限差分方法
在本节中, 将介绍求解线性势流方程组的数值算法, 在交错网格下的矩形单元网格 (如图3) 中定义速度势Φ和速度分量U, W, 其中速度势Φi, jn定义在矩形网格单元中心 (Xi, Zj) 处, 速度分量Ui, jn定义在网格单元中心 (Xi, Zj) 左0.5ΔZ网格处, Wi, jn定义在网格单元中心 (Xi, Zj) 下0.5ΔZ网格处, n表示时间层nΔT。在时间层T= (n+1/2) ΔT和空间位置 (Xi, Zj) 上对方程 (8) (9) (10) 建立Crank-Nicolson隐式格式。
3.1
线性势流方程组离散
得到的差分方程组用超松弛迭代法 (SOR) 求解, 其中0<ω<2为松弛因子, 是Φi, jn的系数, 偏导数的离散, 时间方向用一阶向前差分离散, 空间方向上, 一阶偏导数用两点二阶中心差分格式离散, 二阶偏导数用二阶中心差分格式离散, 边界处用三点二阶差分格式离散。差分格式中时间步长, 空间步长和速度分量的稳定性条件限制为 (Δt<min (Δxmin/|ui, j|, Δzmin/|wi, j|) 。截断误差为O (Δt+Δx2+Δz2) , 当Δx, Δz和Δt趋于零时, 截断误差趋于零。
3.2 迭代算法步骤
设速度势Φ和波高H的迭代收敛性条件为是迭代次数, 在相同的时间层T=nΔT, T= (n+1) ΔT和空间点 (Xi, Zj) 上, 线性势流方程组Crank-Nicolson隐式格式有限差分方法的迭代算法为:
(1) 在正方形计算区域上定义初始U, W, Φ和H的值;
(2) 由自由面动力学边界条件方程 (15) 和固面边界条件 (11) , 计算水槽区域四个边界上的Φ;
(3) 由方程 (14) 计算出区域内部的速度势Φ, 检查Φ是否满足收敛标准, 如果满足, 则进入下一步 (4) , 否则转到第 (2) 步;
(4) 由方程 (16) 计算自由面波高H, 检查H是否满足收敛标准, 如果满足, 则进入下一步 (5) , 否则转到第 (2) 步;
(5) 如果Φ和H都满足收敛标准, 则转入下一时间。
4 数值结果及分析
4.1 数值算法验证
本节将验证数值解的准确性以及与网格的关系。计算出自由表面波高, 给出整个自由表面10秒内两种网格下数值解与解析解的最大误差.设水槽宽度和静水深度分别为 (b2-b1) =1.0m和d=0.5, 水槽自然圆频率, 其中kn=nπ/ (b2-b1) , n=1, 2, …。初始波Y=Acos (knx) , A为振幅, M×N表示网格数, 时间步Δt=0.001sec。Er1, Er2分别表示自由面左端点和右端点波高在10秒内|ηanal-ηnum|的最大误差。振幅A=0.00012m, n=1。在不同网格数下数值解与解析解[3]相比较, 10sec内自由面左端点, 中间点和右端点的最大误差 (如表1) , 从表1中可以看出数值解很好地吻合解析解, 误差随网格的增大而减小。
4.2 数值算例
上节验证了本文建立的CrankNicolson有限差分方法能有效的数值模拟小振幅波。本节取网格大小 (M, N) = (50, 30) 。
图4, 图5可知, 振幅A=0.0015, 0.0025时有限差分数值解与解析解[3]吻合很好。图6, 图7, 图8, 图9给出了不同时刻小振幅波的波形, 结果发现小振幅波表现出线性的特征以及在波节处静止不动。
5 结论
【有限元数值模拟】相关文章：
航空发动机薄壁机匣激光焊接有限元数值模拟05-26
离散裂缝网络模型与有限元方法在页岩气藏数值模拟中的研究10-09
ANSYS有限元模拟07-12
有限空间作业模拟考试05-14
压边圈参数设置对有限元模拟结果的影响09-11
模拟有限责任公司设立程序是怎样的08-15
有限元心得06-23
有限元模型06-09
有限元法08-12