有限元数值模拟仿真

2024-08-06

有限元数值模拟仿真(精选9篇)

有限元数值模拟仿真 篇1

0 引言

手持式电动工具内部零件之间常采用过盈配合[1]来形成传递扭矩或承受轴向力或限制周向位移的能力。正确的确定过盈量,在保证工作可靠性的前提下,从经济性的原则出发,合理的运用标准,选择相应的配合,是过盈配合设计的基本工作。

1 模型建立

1.1 问题描述

手持式电动工具内部的铝壳1、钢管和铝壳2相互配合,如图1所示。其中要限制钢管和铝壳1的周向位移,可以通过适当的过盈配合来实现,同时又要保证装配完成的铝壳1能够和铝壳2装配,此处为间隙配合。

设计工程师根据经验选择原始公差配合尺寸[2]如表1所示。

钢管和铝壳1的最大过盈量为:0.079mm,最小过盈量为:0.029mm。

铝壳1和铝壳2的最大间隙量为:0.075mm,最小间隙量为:0.025mm。

在实际装配过程中发现,当完成钢管和铝壳1的装配后,铝壳1和铝壳2配合很紧,装配很难手动完成或根本不能手动完成。初步分析原因应该是钢管和铝壳1过盈装配后铝壳1的外径增大,导致跟铝壳2的间隙配合尺寸减小。

1.2 几何模型

对于轴对称模型,可以采用传统的计算方法,即根据机械设计手册来进行,但是基于本模型结构的不对称性,采用有限元方法[3]进行建模并仿真计算。

仅截取模型中过盈配合面的部分,简化模型中不影响分析精度的小倒角,这样处理有利于减少单元数量和避免生成质量差的网格。模型的配合面尺寸仅采用基本尺寸,不考虑公差尺寸。

1.3 离散模型

单元类型选择三维实体单元C3D10M,为保证计算精度,对接触面处网格细化。约束条件比较简单,只要限制模型不产生刚体位移即可,本模型在螺栓处施加全约束。接触面之间的过盈量采用*Contact Interference设置,模型的尺寸往往会存在数值误差,所以一般在定义接触时设置一个位置误差限度,用来调整从面节点的初始坐标,这一过程会在初始分析步中完成,不会对模型的应变产生影响[4]。铝壳1的材料为铸铝,弹性模量为71000Mpa,泊松比为0.3,屈服强度为160Mpa;钢管的材料为钢,弹性模量为206000Mpa,泊松比为0.3,屈服强度为685Mpa。接触面摩擦系数在整个配合面表面相同,=0.11。

2 强度计算

本模型对最大过盈量(0.079mm)进行计算。计算结果表明,钢管的最大等效应力小于材料的屈服强度685Mpa,大部分接触面上的应力为220Mpa左右,安全系数大,结构不会屈服,说明该设计结构是弹性范围的过盈配合。最大应力点出现在沟槽处,该处是应力集中现象。铝壳1的最大等效应力大于材料的屈服强度160Mpa,主要出现在边缘处,属于应力集中现象,可能会产生局部塑性变形。

3 变形计算

鉴于铝壳1还要与铝壳2进行装配,此处需知道铝壳1外径变化量。图2显示了数据采集方向。

铝壳1在X轴的径向位移变化范围为:0.028mm~0.044mm。

铝壳1在Z轴的径向位移变化范围为:0.052mm~0.060mm。

径向位移较大值在安装端部,故应取大值为0.060mm。此值已远远超过最小间隙量0.025mm,离最大间隙量0.075也很接近,故按原来公差设计安装可能会有问题。

4 模型验证

为了验证有限元分析结果的可靠性,此处选取两个样本进行测量,内径用内径千分尺测量,外径用外径千分尺测量。尺寸测量结果和试验结果如表2所示:

试验结果表明,过盈量接近0.054mm时,压入力约12k N,铝壳1外径增大量接近0.031mm,与铝壳2的间隙量为0.024mm,可完成手动装配。

对过盈量为0.054mm的模型再次进行有限元仿真,可得结果如下:

铝壳1在X轴的径向位移变化范围为:0.018mm~0.028mm。

铝壳1在Z轴的径向位移变化范围为:0.036mm~0.040mm。

取X轴和Z轴径向位移较大值0.040mm,仿真结果相比于试验结果值有所偏大。

采用动态法仿真压入过程,通过几何尺寸来确定过盈量,压入力约为11k N。仿真结果相比于试验结果值有所偏小。

仿真结果和试验结果的偏差原因可能是多方面的:实际零件测量误差,实际零件存在圆度和同轴度的误差等。因此仿真结果还是具有一定的参考价值。

5 设计改进

参考仿真结果和试验结果的对比情况,更改设计尺寸公差见表3。

钢管和铝壳1的最大过盈量为:0.041mm,最小过盈量为:0mm。

铝壳1和铝壳2的最大间隙量为:0.066mm,最小间隙量为:0.025mm。

对过盈量为0.04mm的模型再次进行有限元仿真,可得结果如下:

铝壳1在X轴的径向位移变化范围为:0.014mm~0.022mm。

铝壳1在Z轴的径向位移变化范围为:0.026mm~0.030mm。

取X轴和Z轴径向位移较大值0.030mm,该值大于最小间隙量0.025mm。鉴于仿真结果和试验结果的对比,实际铝壳1的外径增大量会略小于仿真结果值0.030mm。

零件批量生产后零件尺寸处于极限值的概率很小,当模型钢管和铝壳1处于最大过盈量的同时,铝壳1和铝壳2又处于最小间隙量,这相当于每个零件极限值下的概率串联,实际装配过程中出现这样的情况概率是小之又小,故而这个设计可以予以采用。

选取改进设计后的三个样本进行实验验证,其中铝壳1的内径、外径用三坐标测量,钢管的外径用外径千分尺测量。尺寸测量结果和试验结果如表4所示。

试验结果表明,试样过盈量最大为0.0364mm时,压入力约7k N,铝壳1外径增大量约为0.0198mm,小于铝壳1和铝壳2的最小间隙量0.025mm,可完成手动装配。保守估计在极限配合过盈量0.041mm时,铝壳1外径增大量会有少许增大,但是应该还是在铝壳1和铝壳2的最小间隙量0.025mm之内。再者,要使铝壳1内径最小,外径最大,钢管外径最大,铝壳2内径最小同时满足,这个概率非常之小。所以,该设计被最终采纳。

6 结论

在进行过盈配合设计中,应根据设计目标确定最小过盈量,根据材料强度确定最大过盈量,最终选择合适的配合,确定零件的公差尺寸。在设计过程中,还要考虑装配完成的部件是否还要与其他部件进行装配,如果有,则必须同时考虑过盈配合会带来的包容件外径增大及被包容件内径减小的影响,避免出现与其他零件配合失效的问题。

参考文献

[1]黄德武.过盈配合组合轴对称体有限元法应力分析[J].机械设计与制造,1990,(3):2-5.

[2]濮良贵,纪名刚.机械设计[M].北京高等教育出版社,1996:224-225.

[3]石亦平,周玉蓉.ABAQUS有限元分析实例详解[M].机械工业出版社,2006.

[4]赵华,尹辉.各种柱筒与轴的过盈配合分析[J].石油机械,1998,26(2):7-10.

有限元数值模拟仿真 篇2

某轿车侧面碰撞有限元仿真模拟

为了有效地降低交通事故所带来的.严重后果,我国政府将碰撞保护法规作为强制性检测内容.文章采用试验与计算机仿真相结合的方式.从结构和能量变化2个方面,对试验结果和计算机仿真结果进行了比较和分析,结果表明二者符合较好,证明了仿真结果是真实可信的.指出碰撞仿真分析是确保车辆拥有良好碰撞性能的一种重要方法,与实车碰撞试验结果相比,吻合较好.

作 者:周子云 ZHOU Ziyun 作者单位:天津一汽夏利汽车股份有限公司产品开发中心刊 名:汽车工程师英文刊名:TIANJIN AUTO年,卷(期):“”(5)分类号:U4关键词:侧面碰撞 有限元模型 计算机仿真

有限元数值模拟仿真 篇3

【关键词】有限单元法;数值模拟;结构设计

在力学分析开展的过程中,可以使用的方法非常多,但是从归结方面能够将其分为两个类型,也就是解析法与数值法。可是因为结构物的具体形状与其所受荷载情况的两个方面非常复杂,只有少量非常简单的问题可以使用解析法进行求解,由此导致数值法在力学分析中转变成不能够被取代的且普遍使用的方法,同时还获到连续性发展。有限单元法受到电子计算机技术影响比较大,是以其进步为基础而形成的一种具有新颖特点的数值分析方法。其在数学方面上的逻辑非常严谨,在物理的内容方面极其清晰,技术人员比较容易理解与掌握,由此形成其被使用的范围十分普遍,还可以灵活地对多种比较复杂的问题开展处理与求解,尤其是其使用的矩阵形式表达基本公式,在一定程度上可以方便于快速使用计算机编程进行计算。结构设计的过程中使用有限单元法,其主要把结构物当成是由定量个划分的单元构成的整体,其主要使用单元结点形成的位移或者结点力当做基本的未知量进行求解,此内容其也是有限单元法的基本思路内容[1]。

数值模拟技术对我国目前的建筑结构的科学研究和实际工程使用过程中产生着不可被忽略的作用,其在一定程度上可以推进我国建筑的发展。因为建筑结构的具体体积非常大以及价格较高,开展科学研究时非常艰难或者几乎不可以开展相关的建筑结构模型试验活动,对此,造成了大型的模拟仿真技术对结构设计领有着非常重要的影响。

1.有限元数值模拟技术

我国当前经常使用在工程技术范围里的数值模拟方法包括有四种方法,其分为有限、边界和离散三种单元法以及有限差分法,可是根据方法本身的实用性与使用普遍性的两个方面来说,主要使用的还是有限元法。有限元法的基本思想是把问题的求解区域划分为一系列单元,单元间的连接只依靠节点进行连接。而单元节经过选定的函数关系插值获得的数值,便是单元内部点实际的待求量。因为单元的形状非常简单,方便使用平衡关系或者能量关系在节点量间构建有关的方程式,再把每个单元方程聚集在一起,由此便可产生总体代数方程组,在将边界条件进行计入的前提下,可以对方程组进行求解[2]。问题的求解区域当中的单元划分可对计算结果造成影城,其划分范围越细,计算结果也就越准确。

近年来,随着社会和技术发展的需要,形成数值模拟技术经过计算机程序被普遍地使用在工程上。而90年代的到来,世界上越来越多有限元通用软件被使用在的比较大型的工程中,其种类可达几百种,另外,此类软件的功能随着技术的进步不断被完善,产生了可以对在多种条件下工程开展计算的有限元分析程序,同时还出现了具有功能非常强大的可以进行前后处理的程序。因为该类程序具有使用便捷和计算比较精准两个优势,由此当前各类工业的产品设计与性能分析主要以使用其计算得出的结果作为可以依靠的根据。在工程技术范围内使用有限元数值模拟技术对工程开展求解,能够降低设计的成本投入,减少设计与分析的循环时间,影响产品与工程的可靠性,使其呈现增加的趋势。由此可以知道,限元数值模拟技术使用在结构设计中,可以在开展优化设计时减少工程材料的成本与花费,进行工程施工与产品制造之前可以预先察觉形成的问题。

2.有限元求解方法分析

有限元法使用在工程中,其分析计算的思路与作法主要如下所示:

(1)对单元剖分和插值函数的确定。按照构件自身的几何特性、载荷能力以及设计所要求的变形点,构建包含每种单元的计算模型,然后根据单元的性质与精度形成的要求,将可以表达单元里面每个点的位移函数写出来,主要为u(x,y,z),v(x,y,z),w(x,y,z)或者d=S(x,y,z)α,通过这些函数可以计算出相关的数值。

根据节点位置具备的边界条件,写出用α来表示的节点位移qe=Ca,对C-1以及α=C-1qe两者的具体数值求解出来后,再代入到方程d=Sα中,便可以计算出单元体内任意点的插值,其函数式主要如下:

d=SC-1qe=Nqe

(2)对单元特性分析。根据弹性力学的内容可以知道:

ε=Bqe

那么相应的变分函数为:δε=Bδqe

根据物理关系能够得出应变与应力的关系可以用这个函数表示:σ=Dε=DBqe

按照虚位移原理的相关内容,便可得出单元节点力与位移之间的关系方程式为:?e=Keqe,公式中的K主要代表整体结构的刚度矩阵,使用该方程进行具体的计算便可以知道刚度的具体矩阵。

(3)单元组集。将各单元根据节点组的集成和结构大致相似的整体结构,由此得出整体结构的节点力和节点位移两者之间的关系方程组为:?=Kq

(4) 对有限元方程进行求解。可以使用各种计算方法对有限元方程进行求解,获取各节点的具体位移数值。可是开展解题前,需要对结构平衡方程组的边界条件开展处理,然后再计算出节点位移q。

(5)计算应力。计算出每个节点的具体位移q 后,使用ε=Bqe与σ=Dε两个方程式便可以计算出相应的节点应力。

3.有限单元法在工程中的应用

3.1高桩码头模型

某公司在长江中下游的高桩码头建造工程中,工程结构设计中使用限单元法建立了有限元模型。该码头的长度为350m,宽为30m,其面标高7m,排架之间的距离为8m。码头的桩基均使用规格为600mm×600mm的预制混凝土空心方桩。仅仅只有在输油臂地点面板中间增加设置2根直桩与2根叉桩,剩下的标准段内各排架设置8根桩、2根直桩以及4对叉桩。其纵向梁系和横梁在节点位置进行整体现浇,面板主要使用预制叠合板。

使用有限元对结构开展计算,对于结构比较复杂的类型,其可以开展模拟。传统的理论计算方法对情况比较复杂的结构开展理论计算时,经常将结构的实际受力情况开展合适简化,由此计算出来的结果精准度不高,而使用有限元模拟方法可以尽量地将结构实际的受力情况反映出来,可以确保计算结的准确度。另外,使用过去的理论计算通常不可以容易地知道结构各个细部的受力情况,但有限元方法却可以将结构的各点的受力情况如实表达[3]。构造设计构建有限元模型时,必须按照图纸内容对工程进行三维模型的构建,如此能够经过三维模型对图纸开展检查,确认其正确与否。比如上述提到得到某公司在长江中下游的高桩码头建造工程中,假如高桩码头内叉桩的倾斜角度不准确,非常有可能导致与其相邻的直桩或者叉桩出现在空间上相交的现象,此时能够经过有限元三维模型对角度开展有效的调整。

综上所述,在工程范围内使用有限元数值模拟技术,不仅能够处理好工程中出现的各种结构问题,同时还能够有效地处理传热学、流体动力学以及声学等范围存在的问题。因为其在明显增加高产品设计性能和减少设计时候的同时又可以加强产品在市场竞争中的能力,由此导致目前社会上各类工业产品的设计与性能分析均以使用有限元法计算得出的结果作为可靠的根据。

【参考文献】

[1]陈杰,张治民.有限元数值模拟在精密塑性成形中的应用[J].重型机械科技,2007,10(02):189-190.

[2]魏雪丰,莫璇.有限元强度折减法在滑坡稳定性系数计算中的应用[J].科技信息,2010,15(31):156-157.

土中爆炸衬砌有限元数值模拟 篇4

1 有限元建模

对于爆炸形成的地下球形衬砌,考虑爆炸过程中的材料不均匀,将衬砌分成内层、中层和外层三种力学性能有差异的简单混凝土材料组成的结构体,并研究处于一定埋深的土体包容环境下的应力分布和变形情况,考虑到球形衬砌及工作环境具有在空间对称性的特点,将衬砌模拟简化为平面问题,并取半结构进行分析,半结构及材料分布如图1a)、单元离散如图1b)所示。在图1b)中,衬砌外土体在达到一定厚度后,衬砌受此厚度外土体影响不大,因此在模拟分析时外层土体厚度取衬砌半径的一定倍数便可,本文取衬砌外土体厚为衬砌外半径的15倍。

2 材料参数

根据对爆炸衬砌材料的力学性能试验[3],确定爆炸衬砌各层材料用于模拟分析的力学参数见表1。

3 破坏准则

考虑到混凝土拉压力学性能不同,并且破坏与围压有密切关系等特征,选用能反映这些特点的混凝土Hsieh-Ting-Chen四参数破坏准则[4]作为混凝土的破坏判据。

F(Ι1,J2,J3)=aJ2fc2+bJ2fc+cσ1fc+dΙ1fc-1=0

其中,I1,J2,J3分别为第一应力不变量、第二和第三应力偏量不变量;a,b,c,d均为混凝土材料参数,a=2.010 8,b=0.971 4,c=9.141 2,d=0.231 2。

4 有限元模拟

通过分析可知衬砌本身的变形较微小,整体变形主要体现在土体方面,其原因是由于衬砌和周边土体相比刚度和强度相对较大造成的,故而在此对衬砌的变形不做分析和讨论,仅研究衬砌应力状态的分布及变化。研究其在空容(内部无水)、半容量(水位达到衬砌半径高度)、满容下(衬砌正好充满水)的最大主应力S1和最小主应力S3,为明了起见,应力等值线图形仅绘出衬砌部分。通过比较分析知空容和半容量下爆炸衬砌应力分布差别不大,等值线分布见图2,满容下应力等值线分布见图3。

为了解爆炸衬砌在水压作用下的最佳工作情况和应力分布,以及在极限水压下的工作情况,研究了超高水位(即衬砌充满水后上部仍有一定高度的水位)下的应力分布和极限下的应力状态(见图4,图5)。

5 结果分析

通过图2~图5应力等值线分布可知,衬砌工作应力状态最佳时的水位是在一定超容状态下,此时爆破衬砌应力分布状态最好,最大拉应力不到10 kPa,而最大压应力均处于负值。此时爆破衬砌材料受力处于材料的力学性能安全范围。而当衬砌内水位不断变化时,对衬砌工作情况有如下认识:

1)空容下,爆破衬砌在下部外边缘承受了300 kPa左右的拉应力,在外层衬砌混凝土抗拉强度0.7 MPa范围内。半容时,最大和最小主应力变化不大,且应力分布范围一致。当满容时,应力分布范围、模式和空容、半容时一致,但应力分布范围减小,最大主应力和最小主应力的最大值有所降低。这说明内部贮水情况下,随贮水增多,对衬砌受力有利,初砌趋向安全。空容时衬砌处于较不安全状态,但初砌材料未达到破坏的程度。

2)随水位增加到超容过程中,出现最大拉应力的范围逐渐减小,拉应力降低;同时最大压应力也有一定程度提高,但均没有出现正值,即压应力均处于受压状态。在超容过程中,一定范围的水位使衬砌应力分布最为理想:拉应力基本很小,压应力也较小,这时衬砌材料最为安全,此时的水位为衬砌最安全的工作水位。

3)在水位达到一定程度后,衬砌最大拉应力由衬砌下部外边缘向上发展,达到衬砌下部内边缘处。同时最大压应力的最值也在向衬砌内部发展,并且最小主应力也达到正值,此时对混凝土强度安全趋向不利。

4)在上部承受超压水位较大时,可看出衬砌内部下边缘处最大拉应力将接近混凝土的抗拉强度,可能引起破坏。说明此时爆破衬砌达到了临界承载能力。

6 结语

本文利用有限元对在不同水容量和土体包容环境下的爆炸衬砌进行了数值模拟,通过弹塑性分析获得在不同水容量下衬砌的应力状态变化,通过分析确定了衬砌在空容下工作较为不利,但处于安全工作范围内;在内部贮水增多到一定程度时衬砌趋向安全,但当超过一定水位限度,衬砌工作又趋向不利,且在达到限定水位下衬砌将出现破坏。这对爆炸衬砌的工作状态认识和指导衬砌空间在农田水利中的应用具有很强的指导意义。对应用爆炸衬砌进行雨水期时存水,在旱季时农田抗旱,保障农作用物抗旱保收具有重要意义。

摘要:对土中爆炸砂浆形成的爆炸衬砌应用Hsieh-Ting-Chen破坏准则进行了有限元数值模拟,获得在不同水容量情况下衬砌工作时的应力和变形情况,分析了爆破衬砌随水容变化下的安全情况,解决了爆炸衬砌这种新型结构形式的简化数值模拟问题。

关键词:爆炸衬砌,有限元,数值模拟

参考文献

[1]王海亮,马敏,叶朝良,等.爆炸成腔与爆炸衬砌对土体含水率的影响[J].岩土工程学报,2008,30(184):425-428.

[2]郭增强,赵玉成,王海亮.土中球形爆破衬砌承载力的研究[J].铁道建筑,2004(2):43-45.

[3]王海亮.土中爆炸衬砌应用基础研究报告[R].军械工程学院研究院,2005.

有限元数值模拟仿真 篇5

喷丸强化可有效抑制零部件裂纹的产生和发展,提高其抗疲劳和应力腐蚀能力,近年来在模具、汽车、航空等重要领域得到了广泛的应用。然而,喷丸试验成本昂贵,当前已开始使用计算机数值仿真技术进行研究,如建立准静态模型将有限元法应用于喷丸过程模拟[1~4]; 研究丸粒间距与丸粒半径比为1. 0,1. 5,2. 0时残余应力和塑性应变分布[5]; 研究不同喷射角喷丸时的残余应力场[6]等。

本工作利用LS-DYNA软件,建立了按4 -2 -2 -1分层且按此顺序撞击的9丸粒喷丸有限元模型,研究弹丸速度、重复撞击次数以及前人少有研究的搭接率和丸粒分布等对喷丸残余应力大小及分布的影响规律,为探索相关喷丸参数对喷丸残余应力作用机理及影响效果提供有效分析手段。

1 喷丸强化有限元模型的建立

1. 1 模型描述

图1为使用LS -DYNA软件建立的9丸粒有限元模型。其特征如下: 弹丸的相对位置由下向上分为4层,数量分别为4,2,2,1; 9个弹丸以4-2-2-1的顺序分批撞向工件表面; 工件的尺寸为4 mm×4 mm×2 mm,4个侧面采用非反射边界以避免应力波在边界上的反射,底部为全约束; 对工件中心受喷区域进行单元网格细化,以提高残余应力分析精度[7],工件共由4 000个实体单元组成; 铸钢弹丸直径为1 mm,每个弹丸共由3 500个实体单元组成,弹丸弹性模量为210 GPa,密度7. 8 g /cm3,泊松比0. 3; 利用LS -DYNA中的关键字* CONTACT_ATUOMATIC_SURFACE_TO_SURFACE定义工件与弹丸间的接触算法,模拟弹丸对工件的撞击作用。

1. 2 工件材料的本构方程

为了便于将数值模拟结果与文献[8]的试验结果进行对比以检验该数值模型,工件材料选用广泛应用的高强度铝合金Al 2024 -T3,其弹性模量为73. 1 GPa,密度为2. 78 g /cm3,泊松比为0. 3。

由于撞击过程速度较高、应变率高,需要考虑温度对应力流动的影响。因此,采用能较好处理金属材料高应变率、加工硬化及温度效应的Johnson_cook粘塑性本构关系模型。该模型中屈服应力σy与等效塑性应变的关系为

式中A,B,C,n,m———与应变、应变率及温度相关、通过试验获得的数据,其值分别为265 MPa,426 MPa,0. 015,0. 34,1[9]

式中T———任一时刻的瞬时温度,K

Tm———熔解温度,K

Tr———参考温度,292 K

在LS -DYNA中,应用Johnson_cook本构方程需要提供相应的状态方程( EOS) ,用于表示压强与密度间的关系。本工作使用的Gruneisen状态方程是一种绝热熵增的状态方程,见式( 4) :

式中P———压强,GPa

ρ0———初始密度,g /cm3

μ———ρ/ρ0- 1( ρ为当前密度)

c———粒子速度为0时的波速即声速,cm / μs

γ0———材料的Gruneisen参数

α———对γ0和μ的一阶体积修正系数

E———弹性模量,GPa

S1,S2,S3———曲线拟合参数

其参数如下: a = 0. 45,c = 0. 52 cm/μs,γ0= 2. 02,S1=1. 49,S2= 0,S3= 0。

1. 3 模型验证

为验证该喷丸有限元模型的正确性,根据文献[8]进行仿真,速度66. 2 m/s时最大残余压应力为 - 296MPa,发生在工件中心距表面0. 1 mm深的次表层,残余应力层深度为0. 6 mm; 文献[8]中显示试验最大残余压应力值为 - 313 MPa,发生在深度为0. 15 mm的次表层,残余应力层深度0. 55 mm。比较两者结果,最大残余压应力值仿真结果与文献相差5. 4% ,应力层深度仿真结果相差9. 1% ,相关结果吻合较好,故所建模型准确,可用于相关喷丸参数对残余应力场影响规律的数值仿真研究。

2 仿真结果与分析

2. 1 喷丸速度对残余应力分布的影响

喷丸速度是影响喷丸残余应力分布的主要因素之一,弹丸速度变化范围为20 ~ 300 m/s,而国内现有叶轮式喷丸机喷射的弹丸速度一般低于100 m/s。

通过数值仿真手段可以方便地仿真任意速度的喷丸强化过程,预测其残余应力分布状况。图2为喷丸速度与工件中心受喷点残余应力分布的关系。由图2可见: 速度为25 m/s时,最大的残余应力为 - 510 MPa,残余压应力层的深度约为0. 37 mm; 喷丸速度为100m / s时,最大残余应力值达到 - 620. 9 MPa,残余压应力层的深度约为0. 69 mm。这表明随着喷丸速度增大,残余压应力值及对应深度增加,故适当提高喷丸速度可获得较高的残余压应力及更深的残余应力层。但喷丸残余应力层的厚度十分有限,仅为0. 1 ~ 0. 8 mm。在4种喷丸速度下,弹丸撞击工件表面所形成凹坑的最大深度分别为0. 020,0. 041,0. 062,0. 071 mm。可见,速度越大,凹坑深度也越大。凹坑深度从一个侧面反映了工件表面的形貌状况,工件表面的凹坑越深,相应的表面粗糙度越大。

2. 2 喷丸重复撞击次数对残余应力分布的影响

图3为喷丸速度50 m/s,同一落点重复撞击1 ~ 7次时的残余应力分布。

由图3可知: 随着撞击次数的增加,最大残余压应力增大,且发生在深度为0. 1 mm的次表层; 从第5次撞击后,残余应力的增加幅度很小,是因为反复撞击同一位置使得残余应力趋于饱和,第5次撞击后的工件表面残余应力已基本达到理想饱和状态,此后增加幅度不超过10% 。因此,在喷丸强化过程中,应合理安排喷丸撞击次数,以提高喷丸效率。

图4为重复撞击过程中受喷点凹坑深度的变化。由图4可知: 弹丸撞击过程中深度先达到最大,之后存在一定回弹,深度略有减小; 7次撞击所形成的弹坑深度依次为0. 042,0. 057,0. 067,0. 075,0. 082,0. 088,0. 093 mm,随着撞击次数的增加,其弹坑的深度也不断增加,但其深度增加量不断减小并趋于饱和。

2. 3 搭接率对残余应力分布的影响

两弹丸之间重叠区域的大小对残余应力的影响能间接反映覆盖率的影响。丸粒重叠区域的长度与丸粒直径的比值称为搭接率。当搭接率小于1 /2时仿真结果所形成的两弹坑无重叠; 搭接率大于1 /2时,弹坑发生重叠; 随着搭接率的增大,覆盖率也相应增大,搭接率为1( 2弹丸完全重叠) 时覆盖率为100% 。图5是2个直径为1 mm,速度为50 m/s的弹丸以不同搭接率撞击表面时受喷点残余应力分布。由图5可知: 5种搭接率所对应的最大残余应力分别为 - 556. 5,- 564. 6,-598. 6,- 614. 5,- 625. 2 MPa,最大残余应力数值随搭接率的增大而增大,其均发生在深度为0. 1 mm附近的次表层。

2. 4 丸粒分布对残余应力分布的影响

图6为基于有限元模型不同数量弹丸时的丸粒分布。弹丸速度为50 m/s时,喷射4丸撞击时工件受喷面的典型区域为正方形CEGI,点C,E,G,I分别为4个丸粒的中心,4个1 /4圆面积之和即为实际受喷面积,其覆盖率为78. 5% 。采取分层建模、分批顺序撞击,4种丸粒模型对应图6a中的A,B,C,D 4个典型位置处的残余应力分布见图7。图7a为4丸粒( 按4 -0 -0 -0分布) 模型的残余应力分布,C点为直接撞击点,残余应力最大且与单丸粒撞击后的残余应力分布相近,中心A点无残余压应力,说明4个粒子之间相互影响较小。图7b为6丸粒( 按4 -2 -0 -0分布) 模型的残余应力分布,经过2次撞击后,B点残余压应力最大,C点残余压应力由于粒子之间相互影响而减小,A,D点没有撞击到,残余压应力很小。图7c为8丸粒( 按4-2-2-0分布)模型的残余应力分布,B,D点残余应力值最大,A,C点存在一定残余应力,且最大残余应力值在工件表面。图7d为9丸粒( 按4 -2 -2 -1分布) 模型的残余应力分布,整个受喷面积均被覆盖,各处均具有较大残余压应力且分布较为一致,最大残余应力发生在A点。对比可以发现,增加丸粒数提高覆盖率时,各点的残余应力分布趋于一致,最大残余压应力均发生在0. 1 mm的次表层处,残余压应力层的深度大约为0. 5 ~0. 6 mm。

3 结 论

利用LS -DYNA软件建立的含9丸粒的有限元模型可对喷丸残余应力场进行仿真研究。弹丸速度越高,最大残余压应力及应力层深度越大; 同一位置重复打击可以增强强化效果,但达到一定次数后残余应力将趋于饱和; 搭接率过小将会形成强化盲点,增加搭接率则各点的残余应力趋于一致; 通过增加弹丸数量提高其覆盖率后,各点处残余应力分布趋于均匀和一致。

参考文献

[1]Meguid S A,Klair M S.An examination of the relevance of co-indentation studies to incomplete coverage in shot-peening using the finite-element method[J].Journal of Mechanical Working Technology,1985,11(1):87~104.

[2]Meguid S A,Klair M S.Elasto-plastic co-indentation analysis of a bounded solid using finite element method[J].International Journal of Mechanical Sciences,1985,27(3):157~168.

[3]Khabou M T,Castex L,Inglebert G.The effect of material behavior law on the theoretical shot peening results[J].European Journal of Mechanics A.Solids,1990,9(6):537~549.

[4]Li J K,Yao M,Wang D.Mechanical approach to the residual stress field induced by shot peening[J].Materials Science and Engineering:A,1991,147(2):167~173.

[5]Meguid S A,Shagal G,Stranart J C.3D FE analysis of peening of strain-rate sensitive materials using multiple impingement model[J].International Journal of Impact Engineering,2002,27(2):119~134.

[6]Yang F,Chen Z,Meguid S A,et al.Realistic 3D FE modeling of peening residual stresses of strain-rate sensitive materials with oblique incident angles[A].Proceedings of the2013 Annual Conference on Experimental and Applied Mechanics[C].New York:Springer New York,2013:215~220.

[7]杜平安.有限元网格划分的基本原则[J].机械设计与制造,2000(1):36~38.

[8]Miao H Y,Demers D,Larose S,et al.Experimental study of shot peening and stress peen forming[J].Journal of Materials Processing Technology,2010,210(15):2 089~2 102.

有限元数值模拟仿真 篇6

断裂力学是在材料力学的基础上产生与发展起来的。在近30年来,断裂力学的理论与实验研究迅速发展,并且在工程中成功运用,这为岩石断裂力学的发展奠定了良好的基础,并催生了现代岩石力学的一个重要分支——岩石断裂力学。断裂力学研究的重点对象——裂纹,根据裂纹的力学特征,可将其分为Ⅰ型裂纹(张开型)、Ⅱ型裂纹(错开型)和Ⅲ型裂纹(撕开型)三种;根据裂纹的几何特征,可分为穿透裂纹、表面裂纹和深埋裂纹三种;以及由以上两种或两种以上裂纹组成的复合裂纹。在工程中最常见且最危险的裂纹就是Ⅰ型裂纹,即在与裂纹面正交的拉应力σ作用下,裂纹面产生张开位移(位移与裂纹面正交、沿拉应力方向)。裂纹面上的上表面点和下表面点沿y方向的位称分量v不连续。

在断裂力学的研究中,应力强度因子与物体在外力作用下的变形或应力密切相关,在理论和实践上都有着特别重要的地位。近年来,为得到不同裂纹体的应力强度因子,提出了各种实验和测试方法,应用了包括声发射、全息光弹、激光散斑、光弹性和电阻应变片技术在内的几乎所有现代实验手段,为了得到近真实值,这些方法都对裂纹体的裂纹位置、形状以及边界条件的配置加了较多的限制,实验耗资较大,使它们的应用受到了很大的限制。文中主要阐述其有限元数值分析方法对国际岩石力学学会建议的Ⅰ型断裂韧度测试标准实验方法,即带V型切口的三点弯曲圆棒梁试件进行了有限元法线弹性模拟,求得应力强度因子的过程和方法,说明了有限单元法对Ⅰ型裂纹强度因子求解的可行性及实用性。

1 Ⅰ型断裂韧度测试试件参数

采用国际岩石力学学会建议的标准实验方法,即带V型切口的三点弯曲圆棒梁试件,如图1所示。圆棒直径D=50 mm,跨距S=166.5 mm,初始切口深度a0=7.5 mm,泊松比μ=0.31,弹性模量E=0.68 GPa,最大载荷Fm=0.032 kN,Ⅰ型断裂韧度KⅠ=0.036 MN/m3/2(以上数据及图形来自中南大学郭少华的博士毕业论文《岩石类材料压缩断裂的实验与理论研究》)。

2 实体模型建立及网格剖分

根据郭少华博士对石膏进行Ⅰ型断裂韧度测试设计参数(文中仅借用其数据进行有限元分析,进而做出比较。郭少华博士的论文并没有进行有限元分析),建立了实体模型并进行剖分。为了反映裂纹尖端应力场和位移场的奇异性,在网格剖分时,对裂纹尖端进行了细化加密,增加了近2万个节点。以圆棒轴中心为原点,裂纹尖端位于(0,-13.6 mm)。

3 模型KⅠ型裂纹应力强度因子计算

3.1 Ⅰ型对称型裂纹位移形式的应力强度因子与应力和位移之间的关系式

其中,θ为在该有限元模型下,位于y轴的夹角;ρ为距尖端裂 纹的距离;G为剪切弹性模量,;常数,上者为平面应变问题,下者为平面应力问题。

3.2 位移法求解K

有限元法求解应力强度因子有两种方法:位移法和应力法。1)位移法是在求得尖端附近节点位移后,利用式(1)的后两式估算应力强度因子KⅠ,常常根据裂纹表面θ=π上节点的位移v计算各节点的KⅠ,然后绘制KⅠ—ρ曲线,外推知尖端,求得应力强度因子;2)应力法是在Ⅰ型对称的情况下,根据裂纹尖端附近点处的应力,根据式(1)前三式估算KⅠ,作KⅠ—ρ曲线,外推知尖端,求得应力强度因子。常常根据θ=0或90°上的垂直裂纹面的应力估算应力强度因子。本模型通过位移法估算应力强度因子。

在位移分析法中,裂纹尖端Y坐标为-1.36 cm,选择θ=π裂纹表面上的部分节点,如表1所示。

在位移法中,计算节点的选择既不能离裂纹尖端太近(裂纹尖端有奇异性),也不能离裂纹尖端太远(应力应变在有限的范围内是线弹性的),如图2,图3所示(横向0点表示为-1.36 cm)。

由图2和图3可以看出,从-1.371 8~-1.386 0之间4个点应变近似呈线形变化,而在-1.371 8之前由于应力集中的影响,产生奇异性,而在-1.386 0之后,应变开始缓和,变化微小。因此,本模型选择-1.371 8~-1.386 0之间的节点外推尖端处的强度应力因子(如图4所示,回到0点是因为程序的设置,并非实际模拟情况,即最左边的那条线段,不为参考部分,且0点表示为-1.36 cm)。

由图4可看出,-1.371 8和-1.385 1两个节点外推尖端的应力强度因子,得KⅠ=0.034 MN/m3/2。

4 计算结果对比及总结

通过有限单元法模型求出的KⅠ=0.034 MN/m3/2,这与郭少华博士的实验结果相差近0.002 MN/m3/2,几乎吻合,由此可见,通过有限元法求解KⅠ型裂纹强度因子是切实可行的。但是模型建立时,应注意以下几方面:

1)在网格剖分时,一定要细化和加密尖端裂纹处的网格。在此次裂纹尖端处的KⅠ外推计算时,使用的是相差不到0.15 mm,且离间断约0.1 mm,故在裂纹尖端剖分时,至少使节点间距离小于0.1 mm。但是距尖端较远的地方,对KⅠ的求解几乎没影响,应该尽量减少节点,否则,会很大程度地增加计算时间,重要的是容易产生错误,反而降低了计算精度,但也要注意柔度。

2)在求解KⅠ的过程中,因尽量选择多个节点,绘制路径图应变图,逐渐缩小范围,尤其是裂纹尖端附近更是要加密计算节点,找到线形或近似线形变化的部分,然后外推计算。

5结语

通过有限元模型的建立和分析计算,表明通过现有的有限元软件,数值分析模拟Ⅰ型裂纹,计算应力强度因子,具有重要的意义,确实为确定应力强度因子提供了方便、快捷、经济的工具,对断裂力学在工程中应用的推广有着深远的意义。

摘要:叙述了如何运用有限元方法对KⅠ型裂纹断裂韧度测试试验进行数值模拟,并阐述了运用位移法估算裂纹强度因子的过程和方法,重在阐明使用有限元数值模拟计算岩石Ⅰ型裂纹应力强度因子的准确性和可行性。

关键词:有限元,Ⅰ型裂纹,应力强度因子

参考文献

[1]王勖成.有限单元法[M].北京:清华大学出版社,2003.7.

[2]郭少华.岩石类材料压缩断裂的实验与理论研究[D].长沙:中南大学博士毕业论文,2003.3.

有限元数值模拟仿真 篇7

1 三角形环形单元介绍

对空间轴对称问题,采用圆柱坐标系,r表示径向坐标,z表示轴向坐标,任一对称面为rz面。物体内任一点只有径向位移和轴向位移,且仅与坐标r,z有关,因而轴对称问题可简化为二维三角形截面问题,其程序和平面三角形单元类似。离散轴对称体时,取轴对称体内rz平面的一个截面进行三角形离散化。

由于轴对称问题的环向位移恒等于0,径向和轴向位移不等于0,则取位移模式为:

u=α1+α2r+α3zw=α4+α5r+α6z}

(1)

代入节点位移,解出α1~α6,代入上式,得:

u=Νiui+Νjuj+Νmumw=Νiwi+Νjwj+Νmwm}

(2)

则得到[B]矩阵:

[Bi]=12A[bi0fi00cicibi](i,j,m)

(3)

其中,A,ai,bi,ci均与平面问题三角形单元一致。其中:

fi=air+bi+cizr(ijm) (4)

为计算方便,通常取:

为编程方便,其单元刚度矩阵可写为显式:

[krs]=π¯rA32A[brbs+f¯rf¯s+A1(brf¯s+f¯rbs)+A2crcsA1cs(br+f¯r)+A2crbsA1cr(bs+f¯s)+A2brcscrcs+A2brbs](rs=ijm)

(6)

其中,A1=μ1-μ;

A2=1-2μ2(1-μ);

A3=E(1-μ)(1+μ)(1-2μ)

2 有限元程序介绍

由以上三角形环单元求解空间轴对称问题的理论,用FORTRAN语言编程。程序框图如图1所示。

该程序用一维压缩变带宽存储并组集一维总刚,即利用总刚度矩阵的对称、正定、带状、稀疏的特点,只将它的下三角部分逐行从第一个不为零的元素开始直至对角线元素存储在一个一维数组中,这样再辅以单元的局部自由度所对应的结构自由度数组,即单元的地址就可计算总刚半带宽并确定总刚对角线元素在一维总刚内的位置。这样可大大节约存储空间,提高计算效率。方程的求解采用Cholesky法。整个程序结构简单、条理清晰。

3 理论验证

为了验证计算结果的精度,分析了均布荷载作用下厚壁圆筒的应力,并将其与理论值对比。如图2所示,选取一厚壁圆筒受内压P=2×108 Pa作用,材料弹性模量E=2.0×1011 Pa,泊松比μ=0.3;圆筒内径a=30 mm,圆筒外径b=50 mm,圆筒高l=200 mm。

网格划分如图3,图4所示,单元总数:共8×10=80个,节点总数:5×11=55个,在输入程序中输入有关数据,即可进行有限元计算。

部分节点的计算结果见表1。

由表1可知,在圆筒受内压时,圆筒主要受径向应力和环向应力作用,且径向压力为压应力。

在圆筒半径r不变的情况下,其单元径向应力、环向应力计算值都分别大致相等,只有圆筒底部由于约束作用而有偏差。说明应力值大小与z无关,这与理论情况相符。

4 数值结果与理论值比较

在弹性力学中,当厚壁圆筒仅受内压P时,其应力分量的解析式为:

σr=a2pb2-a2(1-b2r2)σθ=a2pb2-a2(1+b2r2)}

(7)

由式(7)计算出的各节点理论值与数值结果作对比,如图5,图6所示,可知该有限元程序计算出的数值是准确的。

图5和图6分别表示z不变时,径向应力与环向应力的分布规律与理论值的比较。

由此可看出在z不变时,其节点径向应力和环向应力的绝对值随r值的增大而减小,说明其应力的最大绝对值均产生于内壁处。

5工程分析

为了进一步验证该方法在水塔中的应用,对合肥地区某3层居民楼顶一圆柱形水塔进行了有限元模拟,塔身弹性模量E=3.0×102 Pa,泊松比μ=0.2。塔内半径a=2.0 m,外半径b=2.15 m,塔高l=1.6 m,可贮水20 m3,网格划分与厚壁圆筒类似,仍用三角形环单元。荷载简化为三角形线性荷载,节点总数:6×21=126个,单元总数:10×20=200。

由水塔顶部的应变和位移可以看出,在线性分布荷载作用下,水塔主要受环向应力。对z=1.6时节点的环向应力和r方向位移的计算值与由实测应变得到的实际应力值和位移值相比较,计算结果与理论及实际值均符合。当z不变时,节点环向应力由内到外逐渐减小,而位移也是逐渐减小的,与实际相符。

6结语

应用三角形环单元有限元法和在其基础上编写的有限元计算程序可有效解决空间轴对称问题。值得一提的是,这类问题的约束应妥善处理,否则很难得到准确的结果。通过对厚壁圆筒计算,验证了该程序计算空间轴对称问题的正确性,在此基础上计算水塔的受力和位移与实测值相差无几;通过改变荷载的加入方式,还可以进一步计算地震作用下水塔的应力与位移,可为水塔的设计和震害控制提供参考。

摘要:为精确求解水塔应力及变形,介绍了用三角形截面的环形单元计算此类空间轴对称问题的有限元解法,用FORTRAN语言编制了有限元程序,通过程序计算结果与理论值的比较验证了该程序计算空间轴对称问题的精确性。并结合实例进行了说明。

关键词:空间轴对称,有限元,柱坐标,三角形环单元

参考文献

[1]李亚智,赵美英,万小朋.有限元法基础与程序设计[M].北京:科学出版社,2004:153-227.

[2]王元汉,李丽娟,李银平.有限元法基础与程序设计[M].广州:华南理工大学出版社,2001:87-92.

[3]王焕定,吴德伦,林家骥.有限单元法及计算程序[M].北京:中国建筑工业出版社,2004:126-133.

[4]徐秉业,刘信声.应用弹塑性力学[M].北京:清华大学出版社,1995:188-190.

有限元数值模拟仿真 篇8

关键词:隧道,有限元方法,数值模拟,强度安全系数

随着社会的发展,如何解决日益拥挤的交通问题成为城市建设中的重头戏。其中,公路隧道等地下工程成为解决这一难题的关键。现行隧道及地下工程结构设计,要求在经济合理的前提下满足安全性、适用性、耐久性,即满足结构的可靠性。因此以可靠度理论为基础的极限状态设计对隧道及地下工程更为必要。隧道衬砌的投资约为隧道全部投资的1/4~1/3[1]。由于各种人为和客观因素,建成的实态衬砌结构与理想设计条件下的结构相差很大,特别表现在因各种施工因素而引起的衬砌厚度的变异性上。隧道衬砌作为主要承担荷载的受力结构,其厚度变化直接影响到隧道的可靠度和使用寿命,厚度的大小也直接影响了施工成本[2],因此,计算隧道衬砌厚度对隧道结构可靠性的影响,评价衬砌强度的安全性,是亟待解决的问题。国外很多学者对隧道开挖的有限元模拟以及衬砌强度安全性做了相关研究[3,4,5,6]。国内学者在隧道衬砌可靠性的研究上也卓有成效。元松等[7]在基于AN-SYS的高速公路隧道设计与开挖过程仿真分析中利用ANSYS有限元分析软件,对高速公路隧道的设计与开挖过程进行了仿真分析,对隧道围岩地层的应力、位移、变形和内力进行了结构计算,从而对隧道设计参数进行了验证。高波等[8]在隧道衬砌结构可靠指标计算方法的研究一文中,采用有限元-响应面法及Monte-Carlo模拟获得衬砌结构内力的特征参数,并利用五种结构可靠度分析法计算明洞与隧道衬砌结构可靠指标,为今后进一步寻求隧道结构设计安全系数与可靠指标之间的对应关系奠定基础。吴德兴等[9]在基于有限元响应面法的隧道衬砌可靠度分析中,结合具体的隧道工程,运用有限元响应面法对确定的最危险截面进行了数值分析,计算出隧道衬砌可靠度指标。计算结果对确定公路隧道长期安全性的可靠度水平具有一定的参考价值。

本文结合土城子~羊头洼港高速公路工程实例,采用有限元数值模拟,对隧道结构的开挖和支护过程进行二维弹塑性分析,计算现有衬砌结构的受力状态,并通过改变模型衬砌厚度,探讨衬砌厚度变化与隧道结构内力的关系,计算衬砌强度安全系数,并以此对隧道衬砌安全性进行评价,探讨衬砌厚度变化对结构强度安全系数的影响,从而为隧道衬砌设计和施工的安全可靠性提供科学依据和技术指导,并为进一步的工程研究积累资料。

1工程概况

隧道所处区域为低谷、沟谷地带,地形起伏较大。隧道所遇地层岩石呈破坏状,风化蚀变严重。隧道入口在K 10+740~K 10+780附近,出口在K 12+230附近。隧道建筑限界:隧道净宽11.25m,有效净宽10.50m;隧道净高7.305m,有效净高5.0m。隧道采用了净宽11.25m,净高7.305m的3心圆弧拱形式。隧道所处区域的地质剖面图如图1所示。

该隧道采用“新奥法”原理设计,隧道以锚喷混凝土作为初期支护,对于软弱围岩以及断层破碎带采取适当的预支护措施,保证围岩的稳定和初期支护的施作。衬砌结构的断面图如图2所示。隧道衬砌作为主要承担荷载的受力结构,其厚度变化直接影响到隧道的可靠度和使用寿命,厚度的大小也直接影响了施工成本,因此本文设置了不同厚度的衬砌,用有限元方法计算隧道衬砌厚度对隧道结构可靠性的影响,评价衬砌强度的安全性。衬砌厚度变化如图3所示。

2 有限元数值模拟方法

根据地质和设计资料建立有限元模型(图4)。该隧道开挖跨度约为11.4 m,本次建模的范围水平方向取80 m,垂直方向取为60 m。所有的有限元分析均由ANSYS软件完成。计算过程中,假定衬砌为小变形弹性梁,采用BEAM3单元来模拟,衬砌离散为多个连续梁单元, 对衬砌进行受力模拟计算分析。围岩和锚杆支护结构均采用PLANE42加以模拟。计算采用弹塑性分析,围岩材料采用D-P模型[10]。隧道断面属深埋Ⅳ级围岩。初期支护采用锚喷支护,锚杆长300 cm ,喷射混凝土厚40 cm。为简化计算,将锚杆支护范围内的围岩提高20 %[11] ,锚杆间距为100×100 cm。简化后,根据实际的地质情况,结合《公路隧道设计规范》围岩级别以及混凝土有关规定选择计算参数。计算模型约束条件为左右两侧水平约束(即X方向),下侧施加垂直约束(即Y方向)。施加荷载为重力荷载,根据ANSYS软件的规定,具体为施加向上的加速度。各个地层及衬砌的物理力学参数根据隧道勘测和设计资料中的实际力学参数确定。围岩以及混凝土的力学参数如表1,表2所示。

注:锚杆加固区的变形模量取10 GPa (考虑锚杆的加固效应)

2.1 隧道开挖模拟

用ANSYS程序提供的单元进行“生死”处理的功能来模拟隧道的分步开挖和支护过程。所谓单元的“死”,并不是将其从模型中删除,而是将该单元的刚度矩阵乘以一个很小的因子(estif),死单元的载荷、质量等效果将设为0;单元的“生”并不是将单元加到模型中去,而是在计算过程中将杀死的单元在适当的载荷步中重新激活,并按照单元当前所代表的材料类型重新设定单元的刚度、质量、载荷等数值。所以,当初次衬砌还未施工的时候,其对应单元处于“死”的状态,当衬砌施作之后,单元被激活,处于“生”的状态,参与模型的运算。围岩采用深埋Ⅳ级围岩。计算时首先计算岩体的自重应力,然后再根据上述方法模拟开挖过程。整个模拟步骤如下:开挖洞室,用杀死单元来模拟;对洞室进行初期支护,用激活单元并改变单元材料属性来模拟[12]。隧道开挖后的有限元模型如图5所示。

2.2 衬砌结构内力及安全系数

计算中通过设置衬砌梁单元厚度,计算了衬砌厚度变化时的结构内力以及弯矩(图6-图9所示)。

为分析衬砌厚度变化对其安全性的影响,根据上述内力结果计算衬砌结构的安全系数,对衬砌的安全性能进行检验。根据《公路隧道设计规范》规定[13],混凝土偏心受压构件按破坏阶段进行强度验算。具体计算方法为根据材料的极限强度,计算出偏心受压构件的极限承载力Nu,与实际内力N相比较,得出截面的抗压(或抗拉)强度安全系数K,检查其是否满足规范要求[Κ][14],即:

Κ=ΝuΝ[Κ] (1)

(1)式中:Nu为偏心受压构件的极限承载力;K为截面的抗压(或抗拉)强度安全系数;[Κ]为规范要求的强度安全系数。

当由抗压强度控制,即e=ΜΝ0.2h(M和N分别为截面厚度为h时的弯矩和轴力)时:

Nu=φαRabh (2)

(2)式中: φ为构件纵向系数,隧道衬砌取1;Ra为混凝土极限抗压强度;b为截面宽度,通常取1 m;h为截面厚度。α为轴力的偏心影响系数,按以下经验公式确定:

α=1-1.5eh (3)

当由抗拉强度控制,即e=ΜΝ0.2h(MN分别为截面厚度为h时的弯矩和轴力)时

Νu=φ1.75R1bh6eh-1 (4)

(4)式中:φ为构件纵向系数,隧道衬砌取1;R1为混凝土极限抗拉强度;b为截面宽度,通常取1 m;h为截面厚度。

采用以上公式,对隧道拱圈单元的安全系数进行了计算。根据以上数据列表统计了变厚度衬砌在拱脚、拱腰、拱顶单元的安全系数,并绘制了单元强度安全系数图。

在对不同厚度的衬砌安全系数的计算过程中出现的最小安全系数为4.751 , 发生在厚度为45 cm的处于拱腰的衬砌断面,即衬砌断面厚度为45 cm时的初期支护的受力是最不利的,但对照中国《公路隧道设计规范》中的规定,初期支护在该地质条件及施工条件下仍然是安全的。

根据上述内力结果计算衬砌结构安全系数。对衬砌的安全性能进行检验。根据《公路隧道设计规范》规定,混凝土偏心受压构件按破坏阶段进行强度验算。具体计算方法为根据材料的极限强度,计算出偏心受压构件的极限承载力,与实际内力相比较,得出截面的抗压(或抗拉)强度安全系数,检查是否满足规范要求[15]。最后绘出安全系数图(内层实线为初期支护中心线,外层虚线为安全系数)。衬砌支护的安全系数如图10所示。衬砌厚度变化下的结构安全系数如表3所示。

3 结 论

通过对二维有限元数值模拟分析,得出以下主要结论:

1)从衬砌结构的轴力、弯矩图可以看出,对于同一厚度的衬砌结构,边墙的轴力比较大,而边墙与仰拱连接处的弯矩和剪力都特别大,有应力集中现象。

2)在衬砌厚度变化的情况下,随着衬砌厚度增加,衬砌各项内力均呈上升趋势,从而造成衬砌强度安全系数增大,可见衬砌厚度是影响其强度安全系数的一个重要因素。

同时可以看出,用有限元数值方法对隧道开挖支护进行模拟能较真实地反映隧道结构的内力和变形情况,也可以用于对二次衬砌结构安全性能进行评价。

有限元数值模拟仿真 篇9

采空区是指由于人类开采地下矿产而使原来的地层中出现的空洞区。采空区地表可能产生连续性或非连续性变形, 并由此带来一系列岩土工程问题。本文拟从有限元数值模拟计算分析的角度出发, 利用Geo Studio有限元模拟软件分析某煤矿采空区对地基稳定性的影响。

1 地质环境概况

该煤矿矿井开拓方式为片盘斜井, 总体巷道是采用一对主、副巷道平行自井口向地下深部呈“之”形斜井延伸。主斜井井口底板标高为+445 m, 井筒倾角25°, 最终水平标高-200 m。井筒总体斜长近2 000 m, 井筒断面面积5.1 m2, 矿用工字钢支护。

1.1 地形地貌

该区位于构造剥蚀低山丘陵河谷地貌区, 区域总体地形是东北高西南低, 高程介于439.50 m~622.00 m之间, 相对高差为182.5 m。

1.2 区域地质构造

该区位于中朝准地台 (Ⅰ) 辽东台隆 (Ⅱ) , 太子河~浑江陷褶断带 (Ⅲ) , 浑江上游凹陷断带 (Ⅳ) 铁厂~八道江复向斜中部。区内断层构造较为发育, 并伴随有岩浆侵入作用。井田内断层可分为K组和R组两组断裂, 见表1。

1.3 地层岩性

该区地层由老至新有:下古生界:O2m;上古生界:C2b, C3t, P1s, P2s, J3l;新生界第四系 (Q) 。

该煤矿井田内煤系地层主要为C3t和P1s, 岩性主要由灰、灰黑色砂岩、粉砂岩、泥岩、黑色页岩及煤层组成, 煤系地层最大厚度138 m, 最小厚度46 m, 平均厚度84 m左右。与上覆地层P2及下伏地层C2b呈平行不整合接触。

根据钻孔揭露资料本井田内含煤11层, 太原组 (C3t) 含煤6层, 分别为六下层、六层、五下层、五层、四下层、四上层;山西组 (P1s) 含煤5层, 分别为三下层、三上层、二层、一下层、一上层。

1.4 水文地质条件

区内碎屑岩类裂隙水主要赋存于O2m, C2b, C3t, P1s, P2s, J3l和Q4的灰岩、泥岩、砂岩、煤、凝灰岩和砂砾石地层中, 主要接受大气降水渗入补给, 向河谷排泄。区内主要发育有2组断裂构造, 断层破碎带两侧构造裂隙发育, 形成构造裂隙含水带。

2 采空区地基变形的有限元计算分析

2.1 概述

本文采用的Geo Studio系列软件是由加拿大岩土工程软件开发商GEO-SLOPE公司开发的面向岩土、采矿、地质工程等领域的一套仿真分析软件, 主要包含了SIOPE/W, SEEP/W, SIGMA/W, TEMP/W, QUAKE/W及CTRAN/W几个模块。本文使用的是SIGMA/W (岩土应力应变场分析模块) , 选用其中的线弹性 (Linear-Elastic) 模型。

2.2 有限元数值模型的建立

本文选用32号勘探线剖面来建模进行分析计算。

1) 计算假定。在模拟计算中做出了以下假设:a.岩体变形为各向同性。b.岩体的初始应力场由自重应力构成。c.不考虑空间效应, 按平面应变问题处理。d.岩土体均按照线弹性模型建模。

2) 参数设定。所选取剖面高程约为-400 m~480 m, 长约1 340 m。其中新老地层7层, 煤层4层, 岩浆岩1处, 断层5处, 共选用材料介质10种, 各力学指标如表2所示。

3) 模型单元网格划分。本课题岩层选用四边形八节点单元或三角形六节点的自由单元进行剖分。对于四边形单元, 尽量减小其相邻边的差值, 对于三角形单元, 尽量避免小角度锐角出现。在应力梯度较高的区域, 如煤层、断层, 采用较密的网格。基本单元大小定位为10 m×10 m的层次, 共划分了6 316个网格单元, 4 604个节点 (见图1) 。两侧的边界各点设定为X方向固定, Y方向可移动;底部各边界点设定为X, Y方向均固定。

4) 模拟施工方案。由于缺少准确的初始地应力资料, 故仅考虑岩体的自重应力, 忽略其构造应力, 在分析第一步首先计算自重应力场, 其次进行模拟开挖。根据煤层的位置及层厚, 把模拟开挖共分为四步。第一步主要开挖三下层、四下层、五层和六层在-150 m以上的部分, 第二步主要开挖三下层、四下层、五层和六层在-150 m~-250 m之间的部分, 第三步主要开挖五层在K2断层和K10断层之间的部分, 第四步开挖五层和六层在-280 m以下的部分。对模拟开挖的计算完成之后, 进行分级加荷模拟地基在建筑物荷载下的变形, 每级荷载设为15 k Pa, 共4级。

2.3 计算结果

1) 分步开挖时的X位移:随着开挖进行, 岩层扰动区显著增大, 紧邻第三步开挖层的上部岩层出现了明显的X位移;第四步开挖后负的最大X位移为11.34 cm, 出现在开挖层五层、六层之间;正的最大X位移6.59 cm, 在第三次开挖的五层上部岩层。

2) 分步开挖时的Y位移:随着开挖区的增大, 开挖区上部岩层下移逐渐明显, 下部岩层回弹量也持续增加;第四步开挖之后最大沉降位移为26.37 cm, 最大回弹位移为3.90 cm。

3) 分步开挖时的最大主应力:随着临空面的出现两侧岩石向内回弹卸荷, 出现小范围的应力释放;在采空地层角部出现了应力集中现象。最大主应力为158.40 MPa。

4) 分步开挖时的最大剪切应力:随着开挖的进行应力集中现象越来越明显。最大剪切应力为88.54 MPa, 出现在三下层与K2断层结合处。

5) 随着开挖区的增大上部的岩层中塑性区也在增大。

6) 施加四级荷载地表最大沉降分别为26.5 cm, 26.6 cm, 26.8 cm和27.1 cm。

2.4 计算结果分析

1) 应力集中现象。在开采地层的角部位置出现了明显的应力集中现象, 如图2, 图3所示。

现在选取了三下层与R13断层结合部位的几个点 (见图4) , 分析其应力状态, 如图5, 图6所示。由此可知, 在被开采地层的角部岩体受到很大的最大主应力和剪切应力, 须加以支护。

2) 沉降和回弹。由于煤层开采而形成了临空面, 附近岩体因卸荷而回弹、沉降。临空面以下的岩体向上回弹隆起, 但位移并不太明显;临空面以上的岩体由于失去了下部的支撑, 应力状态极大改变, 产生了很大的沉降量。

选取临空面两侧的一系列点来观察它们的Y位移, 可见在接近临空面上的点位移达到了最大值 (分别为3.5 cm和-25.6 cm) , 向岩体内部逐渐减小。显然采空地层附近的沉降量直接影响地表的沉降量。

3) 地表沉降。由于开挖而导致的地表沉降很明显, 总体呈中心沉降大, 向两侧逐渐减小的趋势, 由于左侧岩体质量明显好于右侧, 所以沉降分布并不对称。地表各点的Y位移最大沉降为10.4 cm。在铁路路线附近区域的沉降为10.1 cm左右, 对铁路的影响是显著的。

在地面荷载的作用下地基变形明显 (第四级荷载最大沉降量为17.8 cm) , 但增加的这部分沉降主要由Q4覆盖层的变形引起。

4) 塑性区。在开挖完成之后, 塑性区主要分布为左中右三个区。左区和中区边缘均有临空面的应力集中区, 有松动破坏的可能性, 若破坏将可能使塑性区进一步扩大, 需适当支护;而右区所处环境较稳定, 形成原因可能与K10和F5两个断层有关, 推测可能为F5上部被扯断。

施加地面荷载小于120 k Pa时, 塑性区的扩大并不显著, 在右部J3底部与O2接触的部位塑性区有扩大的趋势, 另外井田边界Q4也产生了塑性区。当地面荷载在120 k Pa以上时, 塑性区将有明显的扩展。

3 结语

1) 从模拟的结果来看, 与实际调查得到的结果相近, 说明所建立的平面应变模型、设置的边界条件以及网格划分是可行的。2) 在开挖煤层后临空面附近会出现应力集中现象, 岩石极易在应力作用下破坏。3) 开挖引起的X位移分布不规律;Y位移 (沉降) 在采空部位上部岩层最大, 向两侧逐渐减小, 向上也有减小, 但不明显;在开挖产生临空面后, 下部岩层向上回弹隆起, 但位移很小, 也没有形成塑性区。4) 地表沉降在采空部位的地面投影处的中心最大, 向两侧逐渐减小, 但由于岩性的差异, 沉降量并不对称。地表的最大沉降达到10 cm以上, 对区内的已有建筑和铁路均会造成影响。5) 经过对模拟结果的分析, 我们认为在开挖引起的变形完成后, 建筑地基在120 k Pa的外载作用下变形主要由第四系覆盖物的变形引起, 如果选择合理的地基处理方案和基础方案, 仍然能够满足稳定性的要求。6) 由于开挖而在区内产生了一定规模的塑性区。施加附加引起的塑性区扩大并不显著。

由于所研究问题的复杂性, 在论文撰写过程中仍存在不少不足和问题, 有待于进一步的研究。主要有以下几个方面:

1) 模型的建立仍然有些理想化, 不能进一步细致的反映问题。2) 在建模过程中, 对剖面的选取要求难以很好的满足。3) 对地层细节的刻画、岩土体物理指标的选用有待进一步完善。4) 2D建模的局限性, 平面问题不能考虑空间效应。5) 所考虑地面荷载均为静荷载, 对动荷载的影响未加以考虑。

参考文献

[1]国家煤炭工业局.建筑物、水体、铁路及重要井巷煤柱留设与压煤开采规程[M].北京:煤炭工业出版社, 2000.

[2]戴华阳.岩层与地表移动的矢量预计法[J].煤炭学报, 2002 (5) :473-478.

[3]王金庄.矿山开采沉陷及其损伤防治[M].北京:煤炭工业出版社, 1995.

[4]佴磊, 于清扬.岩土工程数值法[M].长春:吉林大学出版社, 2007.

上一篇:初中数学推理能力培养下一篇:碳排放权会计