数值有限元方法（共12篇）
数值有限元方法 篇1
0 引言
晶体在应力超过屈服点时会发生塑性变形,其根本原因是在外力作用下,材料的内部临界切应力超过了极限值,导致晶体内部的滑移系开动,位错开始滑移进而产生宏观的塑性变形。以上晶体材料基于滑移系的变形塑性理论早在20世纪30 年代就已经被提出,但并没有被马上应用到有限元模型中,连续体力学的有限元模拟仍然采用各向同性的材料模型。Taylor[1]、Hill[2,3]和Kroner等[4]早期开始着手研究描述晶体中各向异性的塑性,该理论着眼于同质的单一晶粒,仅可以反映在简单加载条件下应力应变的相互作用。为了进一步研究单晶和多晶中复杂的应力应变相互作用以及不同边界条件的影响,晶体塑性有限元方法(Crystal plastic finite element method,CPFEM)被发展起来。
晶体塑性理论可以追溯到Taylor[1]在1938 年提出的晶体塑性的位错运动理论,之后便与有限元方法相结合形成了晶体塑性有限元方法。第一个CPFEM模拟由Peirce[5]于1982年发表。该模型受到当时计算机性能的限制,仅模拟了包含两个对称滑移系的简化单晶的拉伸。得益于计算机技术和性能的迅猛发展,20世纪90年代前后,CPFEM得到了大量的应用。该方法基于单晶塑性变形和位错滑移的实验及理论,结合塑性变形的物理过程和多种晶体塑性公式,将有限元方法应用到介观晶粒尺度的连续体力学研究中。CPFEM结合Voronoi划分技术可以建立二维和三维的多晶模型,用于解决晶粒和亚晶的微观塑性变形问题。
目前所有CPFEM模型的本构关系描述都认为位错滑移是晶粒的塑性变形机制。它的优势在于可以根据用户自定义提供灵活的塑性流动和硬化公式,因而可以在复杂的边界条件下解决晶体的力学问题和织构的形成。国内外运用晶体塑性理论进行了多项研究:Li L等[6]模拟了多晶铝合金的纳米压痕,研究其晶体取向、晶粒形状以及晶界对于压痕的影响;聂君锋等[7]利用基于蠕变的晶体塑性有限元模型阐释了析出相尺寸对于蠕变性的影响;杨竞等[8]建立了FCC结构孪生诱导塑性钢耦合孪生的晶体塑性本构模型,对硬化参数进行了灵敏度分析;赵蒙等[9]以多晶位错滑移及塑性流动为基础,探究了TA15 钛合金在高温变形过程中介观层次上的形变不均匀性以及力学响应。晶体塑性有限元方法已被广泛运用于各种金属多晶材料的介观力学模拟之中。然而,国内外基于位错密度的晶体塑性有限元方法对于镍基In-conel 718合金的研究还较少开展。本实验研究的材料为镍基Inconel 718高温合金,为FCC结构,具有〈110〉{111}的12 个对称滑移系,其成分如表1所示。
本实验深入探讨基于位错密度的CPFEM以及其中的参数,并结合代表体积单元(Representative volume element,RVE)模型以及拉伸试验,对Inconel 718合金的晶体塑性参数进行了标定。为进一步使用该CPFEM研究Inconel 718合金介观晶粒尺度的力学性能提供了重要的参考和依据。
1 晶体塑性有限元方法
1.1 晶体塑性理论
晶体的变形可以分解为两部分,即晶格的伸展变形及转动F*和基于位错滑移运动引起的塑性变形Fp[10],总的变形梯度为:
晶体的弹性性能不受滑移的影响,所以弹性应力都是由F*决定的。Fp的变形率则由滑移速率决定,由式(2)给出:
式中:为滑移系α 的滑移速率,nα为该滑移系的滑移方向,mα为该滑移面的法线方向。基于晶格变形的滑移方向和滑移面法向分别为:
速度梯度L可以定义为:
对于面心立方(FCC)的晶体结构,塑性速度张量Lp可以用FCC所有12个滑移系的剪切速率之和来表示(α=1,2,…,12),如式(5)所示:
每个滑移系的剪切速率可以近似地用幂次定律表示:
式中:速率指数n表示应变速率敏感性。滑移系α的滑移面分解切应力τα表示名义应力在滑移面上的投影,也是塑性变形的驱动力。
滑移系α的强度τcα反映材料对塑性变形的抵抗能力,即在滑移系α达到参考速率所需的应力。基于Harder[11]的硬化方程,位错的切割力是塑性变形的最大阻碍。因此塑性剪切应变率和移动位错密度的平均效果相关。临界剪切应力τcα也就可以写成林位错ρF的关系式:
式中:A是描述不同滑移系之间相互阻碍作用的矩阵。滑移系α的位错的演变可以分解为位错增殖的一项和位错湮灭的一项。如式(8)所示:
式中:b是伯氏矢量,常数yc表示临界位错湮灭的长度,与动态回复相关。
1.2 基于ABAQUS开发的晶体塑性子程序
大型有限元通用软件ABAQUS可以使用子程序(Sub-routine)进行用户自定义的二次开发。其中,子程序接口规范定义了用户自定义材料力学性能(User-defined material mechanical behavior,UMAT)的标准Fortran接口,用于用户自定义材料子程序开发,CPFEM的本构可以通过UMAT输入到ABAQUS。在输入文件中,用关键词 “*USER MA-TERIAL”表示用户自定义的材料属性。运算时,ABAQUS在单元的每个积分点上调用UMAT子程序进行计算[12]。每个增量步开始时,主程序路径经由子程序接口进入UMAT,并将当前积分点的变量初始值传递给相应的变量。用户需要在程序中更新应力以及求解依赖状态变量,并计算一致切线模量来更新雅克比矩阵。在UMAT结束时,更新的变量值将通过接口返回主程序。晶体塑性有限元中,需要输入并更新的参数包括:当前滑移系强度τcα,滑移系的剪切应变γα,滑移系的滑移面分解切应力τα,当前滑移面法向m*α,当前滑移方向n*α,所有滑移系上累计剪切应变,以及每个滑移系上林位错的密度ρFα。通过增量步的变化,积分点上的应力应变以及其他状态参数以增量的形式被逐步求解。
其中,需要在ABAQUS输入文件中定义的输入量为:材料的弹性常量C11、C12、C44;晶体的滑移系包括滑移面和滑移方向,全局坐标和对应的材料取向,应变速率敏感性n,参考剪切应变速率,初始滑移系强度τ0,临界位错湮灭的长度yc,常数K,以及伯氏矢量b。
2 不同晶体塑性参数对材料力学性能的影响
晶体塑性参数由于涉及材料的微观性质,往往不易直接测量,因此需要研究各个材料属性对于材料整体力学性能的影响。其中,N.Zhou[13]利用相场和第一性原理模拟了In-conel 718的相变并计算了弹性常量C11、C12、C44,其结果和实验数据吻合较好。此外FCC晶体的伯氏矢量通过计算得到。其他参数仍无法简单确定,需要进一步研究它们对于材料力学性能的影响。
2.1 单轴拉伸的单晶模型
为了研究晶体塑性各个参数对于材料力学性能的影响,建立如图1所示的单轴拉伸模型,试样呈圆柱的棒状(直径2mm,长5mm),共由2412 个六面体单元组成(C3D8R)。限定模型的自由度使其不会发生刚性位移和转动。轴向拉伸方向为晶体的[110]方向。模拟时采取控制变量法,逐一分析晶体塑性参数对于材料力学性能的影响。
2.2 单轴拉伸模拟结果及讨论
单轴拉伸模型得到的应力-应变曲线如图2所示。与黄永刚[14]的模型不同,该模型的应力-应变曲线不只考虑初始硬化模量h0、滑移系屈服应力τ0以及滑移系突破应力(屈服第二阶段)对力学性能的影响,材料的硬化主要由位错密度来控制,因此应力-应变曲线不是简单的两阶段屈服的应力-应变形式[14],可以很好地反映出非线性硬化过程,计算结果和试验得到的应力-应变曲线更加接近。
2.2.1 初始滑移系强度τ0
初始滑移系强度τ0对晶体拉伸性能的影响仅仅体现在屈服强度的变化,如图2(a)所示,随着初始滑移系强度的增加,材料的屈服点提高。利用最小二乘法进行线性拟合,得到它们的关系式:σs=12.049τ0+15.312。
初始滑移系强度表示在初始位错密度下滑移系抵抗滑移变形的临界切应力。根据Peirce等[5]的硬化理论,当临界分切应力超过初始滑移系强度τ0时,滑移系开动,材料变形进入塑性阶段。因此,材料的屈服点由τ0决定,当材料变形到达屈服点时,开始发生塑性变形,滑移系开动,随着材料的位错密度增加,滑移系的强度增大。
2.2.2 应变速率敏感性n
如图2(b)所示,应变速率敏感性n对晶体的屈服点没有影响,而对晶体的硬化过程有着显著的影响,随着n的提高材料的硬化速率减缓,且极限强度降低。
根据Hintchinson[15]提出的受应变速率敏感材料的硬化方程,如式(9)所示:
式中:f(x)是一个描述应力对应变速率相关程度的普适方程,n是应变速率敏感性,当n→∞时,代表应力与应变速率无关。因此,改变应变敏感速率可以显著地改变材料的硬化模式,当n较小时,材料的硬化过程非常显著,材料对应变速率十分敏感;而当n较大时,材料没有发生明显的硬化,取而代之的是恒定的应力值下的塑性变形,且材料对应变速率不敏感。
2.2.3 参考剪切应变速率
如图2(c)所示,参考剪切应变速率同样对晶体的屈服点没有影响,而是影响晶体的塑性硬化过程,随着的提高,材料更快进入第二阶段屈服,且极限强度降低。
根据施密特定律[16],对于应变速率敏感的晶体,滑移系的滑移速率由相应的滑移面分解切应力τ 决定,如式(10)所示:
式中:是参考剪切应变速率,g是描述当前滑移系强度的变量。可以看到滑移系的滑移速率与参考剪切应变速率成正比。因此,改变可以控制材料滑移系的滑移速率,即材料的剪切应变速率,当较小时,材料更难进入第二阶段屈服,而当较大时,材料很快从屈服阶段进入第二阶段屈服。
2.2.4 临界位错湮灭的长度yc
如图2(d)所示,临界位错湮灭的长度yc只对材料的极限强度产生影响,即当材料的位错的增殖和湮灭达到动态平衡时,拉伸性能达到其极限强度。
由1.1节的晶体塑性理论可知,材料的强度取决于林位错密度,而滑移系中的位错密度变化可以分为位错的增殖和湮灭两个部分,如式(11)所示:
随着临界位错湮灭的长度yc的增加,材料中的位错湮灭作用变得显著,减慢了材料在塑性变形中的位错增殖速度,也使材料的极限强度降低。因此控制临界位错湮灭长度yc可以控制材料的极限强度,材料的极限强度随着yc的增大而降低,yc和极限强度的关系是非线性的,如式(12)所示:
2.2.5 常数K
如图2(e)所示,常数K对于力学性能的影响与yc类似。常数K控制着位错的增殖,而位错的增殖与材料的硬化以及极限强度相关。随着K增大,位错增殖效果减弱,同样减慢了材料在塑性变形中的位错的增殖速度,因此硬化的速率变慢,极限强度也随之下降。
综上所述,各个晶体塑性参数对于材料力学性能的影响可以用图3来表示。晶体塑性参数通过控制材料的屈服强度和硬化过程来实现对材料力学性能的控制。其中,初始滑移系强度τ0控制材料的屈服强度;应变敏感速率n和参考剪切应变速率通过剪切应变速率影响材料的硬化过程;临界位错湮灭的长度yc和常数K通过位错密度影响硬化过程。
3 晶体塑性参数的标定
由2.2小节的讨论可知晶体塑性的各个参数对于材料的力学性能都有显著的影响,同时这些参数都是和材料本身相关的,因此在使用晶体塑性方程前,需要对材料的各个参数进行标定。参数的标定可以通过比较晶体塑性代表体积单元(Crystal plasticity representative volume element,CPRVE)模型的单轴拉伸模拟结果和拉伸试验的实验结果来完成。拉伸试验采用标准拉伸试样,在Zwick Z020万能材料试验机上进行。
3.1 使用Voronoi图建立RVE模型
模拟完整晶粒的拉伸试样的拉伸过程对计算机的要求非常苛刻同时也耗时过长。因此,需要使用RVE模型以最少的晶粒数反映宏观试样的机械力学行为[17,18]。等轴晶粒可以通过三维Voronoi图划分得到[19],其中规定规整度α为最小晶核间距和等价规则截角八面体或十四面体间距的比值。晶核分布的规整度越高,晶粒的规整度α也越高[20]。十四面体划分是最规整的划分,这种划分最接近开尔文多面体,即用最小外表面积的多面体填充整个空间。在14面体划分中,晶核位置在体心立方(BCC)的晶格上。十四面体有14个面,其中8个是正六边形,还有6个是正方形。两个相邻的十四面体共用一个正六边形表面,它们的晶核距离dreg可以用式(13)表示:
式中:Dmean是Voronoi图中的平均晶粒体积。规整度α 则可以被定义为:
δ则是Voronoi图中最小的晶核间距。需要注意的是,α从0到1对应Voronoi图划分的不规则度减小,从完全无序的划分到完全规则的十四面体划分。同时,也可以通过控制最小的晶核间距δ来得到想要的三维Voronoi图的规整度。
因此生成三维Voronoi划分的算法是:考虑在一个三维直角坐标系中,体积为V的立方体区域内,通过依次随机给定晶核的X、Y、Z坐标来确定晶核的位置。这时加入一个限制条件,只有当该晶核到其他晶核的距离大于δ时才接受该随机坐标,否则再次随机坐标值,这样循环往复直至最后N个晶核位置全部产生。需要注意的是,由于控制了最短距离,所得到的规整度会比设定的稍大。但是,随着晶核数量N的增大,规整度会越发接近设定值。通过以上方法得到的平均晶粒度可以表示为:
3.2 RVE模型的敏感度
为了验证RVE模型的可靠性,分别改变RVE模型内晶粒的数量N和晶粒的规整度α,来检验RVE模型的敏感度。
3.2.1 晶粒数量N对于RVE模型的影响
在控制晶粒大小基本不变的情况下,分别建立包含100、250、500、750、1000个晶粒的RVE模型进行单轴拉伸模拟,拉伸方向为X方向,如图4所示,不同的颜色代表不同的晶粒,每个晶粒的取向都是随机得到的,且互不相同。模拟拉伸得到的应力应变曲线如图5所示。可以看到,当模型内晶粒较少时,晶粒数对于RVE模型的影响较大,随着晶粒数的增多,模型在相同应变下的应力值显著增大,然而当晶粒数大于750时,晶粒数对于RVE模型的影响几乎可以忽略不计。因此,750个晶粒是RVE模型晶粒数的临界值。
3.2.2 晶粒的规整度α对于RVE模型的影响
运用Voronoi划分在相同的体积内生成1000个晶粒,设置晶粒的规整度分别为0.2、0.4、0.6和0.8,如图6所示,同样进行单轴拉伸试验的有限元数值模拟,拉伸方向为X方向。得到的应力应变曲线如图7所示,结果表明最大应力不是单调变化的,且规整度对于RVE模型的计算结果影响非常小,当应变量为25% 时,最大应力和最小应力值相差仅1.6%。由于晶粒规整度可以反映等轴晶的晶粒形貌的各向异性,因此该计算结果表明RVE模型降低了单个晶粒的各向异性效应,能准确反映出等轴晶材料在宏观上各向同性的特点。
3.3 Inconel 718晶体塑性参数的标定
用上述RVE模型的拉伸模拟结果和标准拉伸试验比较的方法来标定Inconel 718合金的晶体塑性参数。RVE模型的大小为4.4mm×4.4mm×1.1mm,包含了1000个随机取向的等轴晶粒,模型由98568 个六面体单元(C3D8R)构成,模型的X、Y、Z方向均有一个面限制了相关自由度以避免刚性位移,拉伸沿X方向,如图8所示。
通过改变晶体塑性参数得到大量模拟结果和实验拉伸曲线的比较,最后可以标定出一组晶体塑性参数,见表2。该组晶体塑性参数所得的模拟结果和实验结果的对比如图9所示。在弹性阶段,模拟和实验结果吻合较好,这说明模拟所用的弹性常量C11、C12、C44是比较准确的。而在塑性阶段,模拟的曲线也基本上和实验结果相近,但是受到其本构方程的限制,无法使曲线完全重合。尽管如此,其最大误差仍能控制在小于5%以内。因此标定出的晶体塑性参数已经足以能够反映Inconel 718的力学性能。
使用标定过的晶体塑性参数可以进一步计算Inconel718合金中介观晶粒尺度的力学性能,包括晶粒取向的影响,晶粒形貌的影响,晶界的影响以及其中析出相的影响,也可以加载循环载荷等其他复杂的边界条件,研究微区在复杂载荷下的力学行为。
4 结论
(1)基于位错密度的晶体塑性有限元方法可以准确地反映材料的性能。其中各个参数对材料的力学性能均有显著的影响:初始滑移系强度τ0影响材料的屈服点;应变敏感速率n可以显著地改变材料的硬化方式;参考剪切应变速率可以控制材料的剪切应变速率;临界位错湮灭的长度yc影响材料的极限强度;常数K控制着位错的增殖。
(2)使用三维Voronoi划分可以得到多晶的RVE模型。晶粒数对于RVE模型的影响较大,随着晶粒数的增多,模型在相同应变下的应力值明显增大。而当晶粒数超过其临界值750个时,晶粒数对RVE模型几乎无影响。此外,规整度α对于RVE模型的结果影响不大。
(3)通过对比RVE模拟和拉伸试验的结果对Inconel718合金的晶体塑性参数进行标定。利用标定得到的参数进行RVE模拟,计算结果和试验结果误差小于5%。
(4)标定过的晶体塑性参数可以应用于进一步研究In conel 718合金介观晶粒尺度的力学性能,包括晶粒取向的影响,晶粒形貌的影响,晶界的影响以及其中析出相的影响。
摘要：运用ABAQUS有限元分析软件对基于位错密度的晶体塑性有限元方法(CPFEM)及其晶体塑性参数进行了深入的研究。结果表明,CPFEM晶体塑性本构可以准确地体现材料的力学性能。通过讨论不同晶体塑性参数,得到各个参数可以分别控制材料的屈服强度、硬化过程、剪切应变速率、极限强度等性能。此外,为了标定材料的晶体塑性参数引入多晶的代表体积单元(RVE)模型,并讨论了晶粒数以及晶粒规整度对于RVE模型的影响。结果表明,RVE模型的晶粒数达到临界值750个时能够体现等轴晶的宏观各向同性。结合晶体塑性RVE模拟和拉伸试验结果,对Inconel 718合金的晶体塑性参数进行标定,晶体塑性有限元的模拟结果和实验结果的误差小于5%。证明经过标定的晶体塑性参数可以准确反映Inconel 718的力学性能,也使得进一步研究该合金介观晶粒尺度的力学性能成为可能。
关键词：晶体塑性,有限单元法,代表体积单元,Inconel 718高温合金
数值有限元方法 篇2
基于有限元的一类stefan问题数值模型研究
该文给出基于有限元方法的一类一维stefan问题的数值求解过程及算法.模型的建立基于已知的相变界面和固定边界处测得的温度和热流.模型的精度通过与Neumann获得的解析解的`比较而得到验证.文中所讨论的模型可以用于反Stefan问题中自由边界的实时跟踪或者控制.最后,比较了已有的有限元模型,给出了仿真结果.
作 者：王新房 汪春华 吴光超 WANG Xin-fang WANG Chun-hua WU Guang-chao 作者单位：西安理工大学自动化与信息工程学院,陕西,西安,710048刊 名：自动化技术与应用英文刊名：TECHNIQUES OF AUTOMATION AND APPLICATIONS年，卷(期)：26(5)分类号：N945.12关键词：模型 有限元法 相变
数值有限元方法 篇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 一阶声波波动方程
假设体力为零, 三维一阶声波方程有具体形式如 (1) 式所示:
其中, Vx、Vy、Vz为速度分量, P为地压分量, 为介质密度, vÁ为体积模量, v为纵波速度。
2 高阶有限差分格式
在用交错网格有限差分数值模拟时, 将应力和速度分量分别设置在 和t时刻进行计算[5]采用的空间网格节点分布如图1及表1所示。
相对于空间高阶差分, 时间上的高阶差分会极大地增加计算量, 而时间频散在地震勘探所需的时间采样上并不明显[6], 因此基于运算效率和精度的综合考虑, 本文在时间上采用二阶差分格式, 空间上采用高阶差分格式。
对于微分算子 其2N阶有限差分格式如下[4]:
式中 为差分方程待定系数, 它可通过方程 (3) 进行求解。
综合以上分析可得三维一阶声波方程的2N阶交错网格差分格式如下:
3 吸收边界条件
有限差分数值模拟中常采用的完美匹配层 (PML) 边界吸收条件, 其基本思想是将波场分裂成垂直和平行于传播方向的两组, 在人工边界处快速衰平行界面法向传播的平面波, 对其它方向不衰减, 从而达到减小边界反射的目的, 该方法可很好地吸收不同方向不同频率的边界反射[7,8,9]。对于一阶声波方程, 由于速度分量的计算仅是声压分量与其对应的坐标轴的一阶空间导数, 因此无须对其进行分裂处理。只要对声压分量进行分裂, 将其分解为Px, Py, Pz三个部分, 即:
与其对应的PML边界条件方程可写为:
此外, 三个速度位移方程的PML形式:
(6) ~ (11) 式中d (x) 、d (y) 、d (z) 分别为x, y, z三个方向的衰减因子, Collino推导了其具体表达式如下[10]:
式中vmax为最大纵波速度, 为匹配层厚度, R为理想状态下介质的反射系数。综上分析可得三维一阶声波方程的PML边界条件的差分格式如下:
4 稳定性条件
通过傅里叶分析方法, 可得到三维交错网格高阶有限差分格式的稳定性条件如 (5) 式所示[11]。
式中, t为时间采样间隔, Vmax为计算区域内最大纵波速度, x、y、z分别为x、y、z三个方向的网格间距。
5 模型试算
为了验证算法的正确性, 本文采用均匀介质模型进行了数值模拟和分析。图2所示为一均匀介质模型, 模型尺度为1200m×1200m×1200m, 网格剖分大小为10m×10m×10m, 介质速度为2000m/s, 震源位于坐标 (600m, 600m, 600m) 处。采用主频为30Hz的雷克子波用为震源, 时间步长为1ms, 模拟记录长度为1s, 差分阶数为时间2阶, 空间10阶。
图3是未加吸收边界进t=400ms的波场快照, 可见未加边界处理时, 地震波在模型边界处理发生了反射, 干扰了地震记为, 图4是经PML边界处理后地波场快照, 可见经边界处理后, 在边界处, 地震波得到了有效地吸收, 未发生边界反射。
图4是PML边界处理前后半空间均匀介质的地震记录对比。通过对比可以看出边界处理前, 直达波在边界处理发生了多次反射, 干扰了地震记录, 而经过边界处理后, 地震记录中的边界反射得到了有效吸收, 且不受任何角度和频率的限制。
6 结论
本文从理论上推导了三维一阶声波方程高阶交错网格有限差分正演的差分格式, 并采用完全匹配边界条件对模型边界进行了处理。正演结果表明, 该方法具有空间频散小的特点, 可以很好地吸收衰减模拟边界产生的边界反射, 不受反射角度和频率范围的限制。
摘要：在地震波数值模拟过程中, 采用高阶交错网格有限差分法可以有效地压制频散, 提高模拟精度。本文推导了三维一阶声波方程的空间任意偶数阶有限差分格式, 并利用完全匹配边界条件对三维地质模型边界进行处理。数值模拟结果表明, 该方法可以有效地吸收人为边界反射, 不产生任何边界反射, 取得了良好的吸收效果。
关键词：一阶声波方程,交错网格,PML边界条件
参考文献
[1]杜增利, 李亚林, 尹成, 高宏亮.虚谱法一阶应力-速度方程地震模拟[J].石油地球物理勘探, 2009.10, 44 (5) :637.
[2]陈耿毅, 余钦范, 蔡希玲等.地震模拟技术新进展-第67届EAGE年会论文综述[J].勘探地球物理进展, 2005, 28 (6) :439~448.
[3]张永刚.地震波场数值模拟方法[J].石油物探, 2003.6, 42 (2) :143~148.
[4]和少伟.弹性介质地震波场的数值模拟[D].荆州:长江大学.2012.4:1~2.
[5]董良国, 马在田, 曹景忠, 王忠华等.一阶弹性波方程交错网格高阶差分解法[J].地球物理学报, 2000.5, 43 (3) :413.
[6]邓帅奇.全空间弹性波场数值模拟与逆时偏移成像方法研究[D].徐州:中国矿业大学.2012.5:23.
[7]牟永光, 裴正林.三维复杂介质地震数值模拟[M].北京:石油工业出版社, 2005:182.
[8]王永刚, 刑文军, 谢万学, 朱兆林.完全匹配层收边界条件的研究[J].中国石油大学学报 (自然科学版) , 2007.2, 31 (1) :19.
[9]Collino F and Tsogka C.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic hetero geneous media.Geophysics, 2001, 66 (1) :294~307.
[10]夏凡, 董良国, 马在田.三维弹性波数值模拟中的吸收边界条件[J].地球物理学报, 2004.1, 47 (1) :132.
地面雷达数值风洞技术实现方法 篇5
地面雷达数值风洞技术实现方法
基于Fluent等软件平台和standard k-e湍流模型,介绍了地面雷达数值风洞技术的实现方法,同时对模型建立、网格划分等关键技术进行了详细介绍,通过实例进行了计算并将计算结果与试验结果进行了对比分析和验证,结果表明,采用数值风洞的理论和方法,实现地面雷达风洞试验的数值化是可行的.,误差是可以接受的,最后对雷达系统数值风洞技术的未来发展作了展望.
作 者：梅启元 胡长明 秦国良 MEI Qi-yuan HU Chang-ming QIN Guo-liang 作者单位：南京电子技术研究所,江苏,南京,210013刊 名：电子机械工程英文刊名：ELECTRO-MECHANICAL ENGINEERING年，卷(期)：24(5)分类号：V211.74关键词：Fluent 湍流模型 地面雷达 风洞试验 数值风洞
多种微分方程数值计算方法分析 篇6
关键词：偏微分方程；数值计算；有限元方法；有限差分法
一、引言
在工程应用、科学研究中，所建立的数学模型很大程度上都是偏微分方程。但是，除了极少数特殊类型的偏微分方程能用解析的方法求得其精确解外，大多数情况下要得出解的解析表达式是非常困难的，因此，就需要应用数值方法来近似逼近精确解。
大规模的计算中，往往要用到计算机，然而，电子计算机只能储存有限个数据，做有限次运算，所以任何一种用计算机解题的方法，都必须要把连续问题离散化，最终化成有限形式的线性代数方程组。
二、有限元方法的基本理论
1.基本步骤及理论。有限元法是求解边值问题的数值方法。有限元法求解偏微分方程的基本思想就是传统的Ritz-Galerkin法，但是它运用样条函数提供了一种选取“局部基函数”或“分片多项式空间”的技巧，克服了Ritz-Galerkin法选取基函数的固有困难，它已成为求解偏微分方程，特别是线性椭圆型偏微分方程的一种有效的数值方法。
有限元方法的基本步骤可以归纳如下：
（1）把原问题转化为变分形式。
（2）选定单元形状，对区域进行剖分。
（3）构造节点基函数，形成有限元空间。
（4）以某种方法给出单元各状态变量的离散关系，形成单元刚度矩阵，并且组装成总刚矩阵，有限元法最终导致联立方程组，求解各节点处的函数值。
（5）收敛性及误差估计，对于计算所得的结果，将通过与设计准则提供的准许值比较来评价并确定是否需要重复计算。
2.数值算例。下面以一类二维带复合边界条件的偏微分方程 进行有限元解法的说明。
-αΔu+βu=s，x∈Ωu=d，x∈Dα■+?酌u=r，x∈?祝 （1）
其中，■为u沿单位外法向的方向导数，Ω为有界区域，边界D∪?祝，?鄣Ω，α，β，?酌，r，s，d为给定常数。
本文对区域采取三角剖分，在每个单元（e）上，在该单元上的有限元解可以表示为：u=uiφi+ujφj+ukφk。其中，φi为节点（x，y）的基函数，选取如下：
φi=ai+bix+ciy （2）
满足：φi（xi，yi）=1，φi（xj，yj）=0，φi（xk，yk）=0，将上述条件代入（2），可具体求得基函数的表达式。同理，其他节点处的基函数也是如此形式及求法。
下面将是在变分形式下，推导出方程（1）的有限元解法的代数方程组形式，即方程（1）的弱解。
设u是方程（1）的解，则u在任意单元（e）上也满足方程（1），记单元（e）为Ω，则变分形式如下：
■（-αΔu+βu）vdσ=■svdσ，?坌v∈H10（Ω'） （3）
运用Green公式，（3）可以化为：
■α?塄u·?塄v+βuvdσ=■svdσ （4）
令u=■uiφi，v=φj，代入上式，整理得：
■ui■[α（■■+■■）+βφiφj]dσ=
■sφjdσ （5）
其中，对j=1，2，3均成立。写成矩阵的形式即：
（αA+βB）U=b （6）
其中，U=（u1，u2，u3）T为待求矩阵，
Aij=■■■+■■dσ （7）
Bij=■φiφjdσ，bj=■sφjdσ （8）
根据上述，对每一个单元（e）都可以求出如式（7），（8）形式的矩阵，即为单刚矩阵，然后将各单元的单刚矩阵组装成总刚矩阵，并根据边界条件删除相关的节点，即可得到最终的线性方程组。
三、有限差分法基本理论
1.基本步骤及理论。有限差分法是微分方程和积分方程数值解的一种重要方法，基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组。
有限差分法的基本步骤可以归纳如下：
（1）对求解区域进行网格剖分。
（2）构造逼近微分方程定解问题的差分格式（差分格式有多种，可以根据方程的类型选择适当的差分格式）。
（3）差分方程的解法，差分解的存在唯一性，收敛性及稳定性的研究分析。
2.数值算例。下面以一类一维四阶半线性抛物型方程进行有限差分法的说明。
u1+Δ2u=|u|p，（x，t）∈（0，1）×（0，1]u（x，0）=u0（x），x∈（0，1）u=0，■=0，t∈（0，1]，x=0，x=1 （9）
其中，p＞1，Ω是一维有界区域，?鄣Ω为Ω的边界。
下面是对式（9）进行构造差分格式。对时间区间[0，1]，N等分，时间步长Δt=■；对空间区域[0，1]作M等分，空间步长h=■，且满足Δt=o（h2）；记xm=mh，tn=nΔt，网格Ih={xm|m=0，1，2，…，M}，φ为定义在Ih上的网格函数，定义如下差商：
δ-φm=■，m=1，2，…，M；
δ+φm=■，m=0，1，…，M-1；
δ+δ-φm=■，m=1，2，…，M-1；
δ+2φm=■，m=0，1，…，M-2；
δ+2δ-φm=■，m=1，2，…，M-2；
δ+2δ-2φm=■，m=2，3，…，M-2。
设Umn是精确解u在（xm，tn）处的离散解，则利用上述差商的定义，给出方程（9）的有限差分格式为：求{Umn|m=0，1…，M；n=0，1，…，N}满足
■+δ+2δ-2Umn=|Umn-1|p，m=2，3，…，M-2；n=1，2，…，NUm0=u0m，m=0，1，…，M （10）U0n=UMn=0，■=■=0，n=0，1，…，N 其中，u0m=u0（xm）。该格式的收敛性分析见文献。
四、结论
作为两种广泛应用的数值解法，有限元方法和有限差分法在其基本思想上有着共同的地方。从上述论述中可以看出，有限元法和有限差分法都是将连续问题离散化，其离散主要步骤都是：首先，对求解区域进行网格剖分，用有限个节点来代替连续的区域；其次，将微分算子离散化，从而把微分方程的定解问题化成线性代数方程组的问题来求解。显然，这两种方法的主要区别在于有限元方法是从定解问题的变分形式出发，用Ritz-Galerkin方法推导出相应的线性代数方程组，而且基函数要特定选取，而有限差分法是从定解问题的积分或者微分形式出发，用数值积分公式或者数值微商推导出相应的线性代数方程组。
参考文献：
[1]李荣华.偏微分方程数值解法[M].北京：高等教育出版社，2005.
[2]王烈衡.有限元方法的数学基础[M].北京：科学出版社，2004.
[3]张文生.科学计算中的偏微分方程有限差分法[M].北京：高等教育出版社，2006.
[4]李德元，陈光南.抛物型方程差分方法引论[M].北京：科学出版社，1998.
土中爆炸衬砌有限元数值模拟 篇7
1 有限元建模
对于爆炸形成的地下球形衬砌,考虑爆炸过程中的材料不均匀,将衬砌分成内层、中层和外层三种力学性能有差异的简单混凝土材料组成的结构体,并研究处于一定埋深的土体包容环境下的应力分布和变形情况,考虑到球形衬砌及工作环境具有在空间对称性的特点,将衬砌模拟简化为平面问题,并取半结构进行分析,半结构及材料分布如图1a)、单元离散如图1b)所示。在图1b)中,衬砌外土体在达到一定厚度后,衬砌受此厚度外土体影响不大,因此在模拟分析时外层土体厚度取衬砌半径的一定倍数便可,本文取衬砌外土体厚为衬砌外半径的15倍。
2 材料参数
根据对爆炸衬砌材料的力学性能试验[3],确定爆炸衬砌各层材料用于模拟分析的力学参数见表1。
3 破坏准则
考虑到混凝土拉压力学性能不同,并且破坏与围压有密切关系等特征,选用能反映这些特点的混凝土Hsieh-Ting-Chen四参数破坏准则[4]作为混凝土的破坏判据。
其中,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.
数值有限元方法 篇8
喷丸强化可有效抑制零部件裂纹的产生和发展,提高其抗疲劳和应力腐蚀能力,近年来在模具、汽车、航空等重要领域得到了广泛的应用。然而,喷丸试验成本昂贵,当前已开始使用计算机数值仿真技术进行研究,如建立准静态模型将有限元法应用于喷丸过程模拟[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.
数值有限元方法 篇9
手持式电动工具内部零件之间常采用过盈配合[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.
数值有限元方法 篇10
某水电站将修建在青海省一条河流上,坝址是选取目前的中坝址方案。要研究的工程部位属于该电站预选坝址方案中的右岸导流洞进口边坡。边坡整体坡向为SW214°35′2″,平均坡角43.3°,边坡高度为195 m。本区位于青藏高原北部,干旱少雨,植被较不发育,出露的岩层主要是二叠系下统灰黑色千枚状板岩。据最新观测资料该河现今水位为2 595.7 m,正常蓄水位是2 715 m,谷顶高程2 901 m左右,与现今河水位高差达到305.3 m。
20世纪80年代以来,高边坡稳定性问题已经成为我国最具现实意义的重大工程地质问题之一[1]。导流洞进口边坡这类高边坡的稳定性对水电站的施工安全和今后的运营都至关重要,所以分别从岩体质量和有限元数值计算上对该边坡的稳定性进行研究。
2 边坡岩体质量研究
合理的岩体质量分级,对于客观反映边坡岩体的固有属性、深入认识岩体力学特性和合理选取参数、制定岩体工程设计和施工方案、采取合理的工程处理措施是十分重要的[2]。对工程边坡进行岩体质量分级,也可以从宏观上判断其稳定性。目前在工程边坡领域,国内外常用的综合岩体质量分级方法有南非学者比尼奥斯基(Bieniawski)的RMR法[3];我国学者在“八五”期间提出了边坡岩体质量分级的CSMR(Chinese System for SMR)分类体系,并将其初步应用于我国的水电工程边坡中[4]。野外实地调查现有的勘探平硐,对右岸导流洞进口边坡岩体质量进行研究,具体评价结果如下表1所示。
3 边坡稳定性有限元分析
以往研究边坡稳定的极限平衡方法,没有考虑到边坡内部的应力与应变的关系。随着计算机技术的飞速发展,能很好模拟应力与应变本构关系的有限单元法,近年来成为分析边坡稳定性的新趋势。
按有限元计算得到的应力结果对区域内各单元进行强度判别,根据破损区的范围和分布评价边坡的稳定性[5]。根据应力分析结果进行边坡的整体稳定性分析,Donald提出了折减抗剪强度进行有限元分析,根据节点位移发展评价边坡稳定性的NDM方法[6]。在国内应用方面,桂重等人通过研究实际工程证明了有限元强度折减寻找边坡主要失稳模式方法的可行性[7];张燕等人采用有限元法分析引水发电洞围岩开挖后是否支护两种情况的围岩稳定性,对设计与施工提出了支护建议[8]。
3.1 计算条件
右岸导流洞工程边坡区域内构造比较简单,长大的断裂构造不发育,主要表现为小型的裂隙和挤压带。该边坡的卸荷也不太发育,仅在坡表有浅层的卸荷,深约5～6 m。
因此基于以上风化卸荷带的划分,该边坡岩体主要由3种材料构成,由坡表向坡内分别是卸荷岩体、弱风化岩体和微新岩体。这3种材料的物理力学参数结合试验指标,再参照类似工程经验,最终选取结果如表2。
3.2 计算模型的建立
通过有限元计算验证该边坡的稳定性,建立导流洞边坡岩体应力-应变有限元数值分析模型如图1。该模型的范围,底部和顶部边界别分取2 487m和2 740m的高程位置,两侧边界从坡脚处开始沿导流洞轴线方向向坡内取长度206 m。模型的两侧边界和底部边界均采用位移约束边界。
计算模型岩体由上述3种材料单元构成,材料的本构关系采用的是摩尔-库伦准则。该模型采用三角形参数应变单元,模型分析容差0.001,计算时将其视为二维平面应变问题,其规模为:单元数1 661,结点数886(开挖前);单元数1 464,结点数780(开挖后)。
3.3 计算工况的选取
采用有限元2D软件,模拟了边坡开挖过程和地震效应,进行边坡应力位移场仿真计算,共选取2种工况作为计算方案:①天然工况,仅以边坡岩体的自重应力场作为初始应力场;②天然+地震工况,据地震局颁布的《中国地震动峰值加速度区划图》(GB18304.6.1.2001),工程区50年超越概率10%的地震动峰值加速度为0.15 g,相应地震基本烈度为Ⅶ度。
3.4 计算结果与分析
3.4.1 天然工况的成果分析
(1)应力场特征。计算结果表明,天然状况下边坡最大主应力分布受重力场控制,呈现从上至下,由坡面向坡体内部逐渐增大的规律分布,在模型计算范围内最大可达5.60 MPa;开挖后,边坡内部的应力发生重新分布,于开挖坡脚处发生应力集中,最大主应力值一般为3.50～4.20 MPa,最大能达4.55 MPa(图2)。在马道转折处也出现了小范围的应力集中区域,其最大主应力值为1.05 MPa左右。
在自然边坡的坡表均未出现有拉应力;开挖后,在工程边坡的下部坡表和进水底板面都出现了拉应力,其中在进水底板面前缘的拉应力约为-0.10 MPa(图3)。
(2)开挖后位移变形特征。边坡开挖后产生位移变形的部位主要位于导流洞的进水底板前缘,最大水平位移为3.80 mm,最大垂直位移约为7.70 mm(图4);最大总位移约为8.40 mm(图5)。
(3)开挖后安全系数特征。利用有限单元法,通过强度折减,使模型达到不稳定状态时,有限元静力稳态计算将不收敛,边坡出现底部到顶部的塑性贯通区,此时的折减系数就是传统意义上的安全系数[9]。
首先选取初始折减系数,然后计算运行,模型开始自动不断地折减其上所有岩体的强度参数(黏聚力C和内摩擦角φ),使系统达到不稳定状态,有限元将不收敛。从边坡最大总位移与强度折减系数(Strength Reduction Factor,缩写为SRF)的关系曲线(图6)可以直观的看出,SRF增加到1.78以后,最大总位移量突然急剧的上升,有限元计算呈现收敛的特征,因此该边坡的稳定性安全系数即为1.78。图7表示整个边坡的剪应变特征,图中展示当SRF=1.78时,岩体抗剪强度不能承受其所受的剪应力,坡体内出现了一条贯通边坡上下的塑性区,此时边坡将处于极限稳定状态。
3.4.2 天然+地震工况的成果分析
地震只能是在某一特定时刻发生,并作为一种外动力作用施加于边坡,为了考虑地震对边坡应力场、位移场的影响,应在开挖末期施加地震力[10]。
(1)应力场特征。在地震作用下,坡体内的应力场没有显著的变化,只是坡脚处的最大主应力由4.50 MPa增大到4.80 MPa,说明应力集中更加明显;而马道转折处的最大主应力则大致减少了0.05～0.15 MPa(图8)。而最小主应变化也不太大。
(2)位移变形特征。施加地震力后,边坡坡体内的位移场也没有发生很大的变化,主要位移变形部位还是在导流洞的进水底板前缘,只是最大水平位移变为4.40 mm,但最大垂直位移还是7.70 mm(图9);最大总位移则增大到了8.80 mm(图10)。
(3)安全系数特征。从图11可以得出,SRF增加到1.69以后,最大总位移量突然急剧的上升,有限元计算呈现不收敛的特征,因此该边坡在地震时的稳定性安全系数即为1.69。图12表示了对应的整个边坡的剪应变特征。
4 结 论
通过以上的分析,得出如下结论:
(1)岩体质量分级的关键问题是,对岩石强度、RQD、结构面间距以及风化程度等核心分级指标取值方法的合理选取,这需要结合野外详细调查进行,使得分类结果与实际相符合。最终按照CSMR法,右岸导流洞边坡的岩体质量评价结果是:31～49 m的洞段是Ⅲa类岩体,其余洞段都以Ⅲb类为主。边坡岩体质量一般。
(2)有限元数值模拟结果表明,边坡开挖后应力集中最明显的部位出现在工程边坡的坡脚;而位移变形相对突出的部位为导流洞的进水底板前缘。但开挖引起的位移值不超过1 cm,因此边坡的变形特征是以卸荷回弹为主。
(3)地震作用下,边坡内的应力场和位移场变化多不大,只是起到对以上现象的轻度加强作用。地震对边坡水平位移的影响较之垂直位移要显著的多,这说明地震作用相当于对边坡施加了一个水平方向的荷载作用。
(4)该水电工程规模为一等大Ⅰ型工程,主要建筑物级别为一级,按《水利水电工程边坡设计规范》(DLT 5353/2006)的规定,该边坡属于Ⅰ级枢纽工程区边坡,其在持久状况下的设计安全系数为1.30～1.25。采用有限元强度折减法分析该边坡的稳定性,得到边坡开挖后的安全系数,在天然工况下为1.78,地震工况下为1.69,都满足规范对边坡安全设计的要求,可见边坡整体处于稳定状态。
综上所述,右岸导流洞进口边坡,除了坡表存在一些受裂隙切割的不稳块体外,未有整体失稳的破坏迹象。参照《水利水电工程地质勘察规范》(送审稿)的规定,该类边坡开挖洞室时,一般的支护方式是:喷混凝土、系统锚杆加钢筋网;采用TBM掘进时,需及时支护;跨度为20～25 m时,还需刚性支护。
参考文献
[1]黄润秋,邓荣贵等.高边坡物质全运动过程模拟[M].成都:成都科技大学出版社,1993.
[2]李天斌.岩质工程高边坡稳定性及其控制的系统研究[D].成都:成都理工大学,2002.
[3]Bieniawskizt Z.T.Engineering Rock Mass Classification[M].New York:Intercience Publication.Viley,1993.
[4]陈祖安.中国水利发电工程(工程地质卷)[M].北京:中国水利出版社,2000.
[5]钱家欢.土工原理与计算(上册)[M].北京:水利电力出版社,1982.
[6]Donald I B.Finite element assessment of slope Stability[R].De-part ment of Civil Engineering Monash University,Melhourne,Australia,1993,D1-D18.
[7]桂重,陈建康,吴震宇,等.基于可靠度的边坡稳定风险分析新方法[J].中国农村水利水电,2008,(9):114-116.
[8]张燕,段亚辉.江坪河水电站引水发电洞围岩稳定性分析[J].中国农村水利水电,2009,(5):78-81.
[9]赵尚毅,郑颖人,时卫民,等.用有限元强度折减法求边坡稳定安全系数[J].岩土工程报,2002,24(3):333-336.
数值有限元方法 篇11
摘要：本文从案例教学法的发展现状出发，结合数值计算方法课程的特点，简述了数值计算方法课程引入案例教学法的必要性，提出了案例教学法的研究内容和研究方法，最后总结了案例教学法在课程建设中的作用。
关键词：数值计算方法、案例教学法、案例资源库
中图分类号：O241-4
一、 引言
案例教学法是一种以案例为基础的教学法。教师在教学中扮演着设计者和激励者的角色，鼓励学生积极参与讨论。案例教学法起源于上世纪二十年代，由美国哈佛商学院所倡导，当时是采取一种很独特的案例形式的教学，这些案例都是来自于商业管理的真实情境或事件，通过此种方式，有助于培养和发展学生主动参与课堂讨论，实施之后，颇有成效。
二、《数值计算方法》课程运用案例教学法必要性
《数值计算方法》课程既有数学类课程中理论上的抽象性和严谨性，又有实用性和实验性的技术特征， 是一门理论性和实践性都很强的学科。目前利用传统的教学方法在数值计算方法课程的教学中反映出如下几点问题：
（1）学习兴趣淹没在冗长的公式推导和理论分析之中。
（2）学生个体差异被忽略在统一的授课内容和进度之中。
（3）课程内容难以体现贯通培养课程体系的桥梁性作用。
三、《数值计算方法》课程案例教学法的研究内容
对数值计算方法案例教学做系统的理论研究和应用研究，揭示数值计算方法案例教学中的规律性，教师要把繁杂、枯燥的理论知识与相应的实践环节有机地结合起来，既要让学生学习理论知识，又要激发学生的学习兴趣，提高认知水平，培养创新能力。在对《数值计算方法》课程进行案例教学法的研究中，主要从以下几个方面进行。
首先，对《数值计算方法》教学案例资源的建设，即对案例教学内容的设计与组建。主要是其中包括
① 案例的适用性的分析。分析原始素材作为案例是否适合教學选用，选择与知识点相匹配，适合教学选用的素材作为案例，并进行适用性的分析。比如，在介绍拟合问题的时候我们在选择教学案例的时候，就要考虑到利用拟合原理解决的实际问题很多，例如，热敏电阻电阻值的变化规律，材料学中混凝土泌水率的曲线拟合问题，轧钢板型问题……，但是要筛选出适合教学使用的，相对完整的案例作为教学素材，并在选用前进行适用性分析。
② 把原始素材制作成和课程匹配的案例。按照规范好的格式对已有原始素材进行整理、收集和分类汇总。在教学过程中，平时注意积累和收集与数值计算方法相关的工业生产、经济生活中案例素材，将这些案例素材整理，汇总并进行分类，比如，利用插值法解决的工业生产中的问题归为一类，诸如矿井突出后的瓦斯涌出量的计算、变形抗力模型等。
③ 案例的分层筛选。依据学生的学习要求，对案例进行分层筛选，比如，第一层案例适用于工科专业学生教学选用，以工科所学专业知识作为背景的案例，第二层适用于数学和信科专业学生教学选用，第三层适用于研究生教学选用，也就是说，选取特定的案例为特定的教学对象服务。
其次，对案例教学环节的设计
对课堂的案例教学环节进行设计，包括如引入、分析和总结。案例的引入，能激发学生产生对学习内容的兴趣，案例的分析，能提升学生对问题的分析能力，案例的总结，引领学生对未来的学习进行展望。通过案例教学能有效调动学生学习的积极性，启发学生对学习内容进行深入的探究。
最后，案例资源库的建设
建设与案例教学相关的资源库，将经过精心筛选和设计的案例放在网络平台上，供学生学习使用，真正实现在网络平台上的资源共享。
四、主要研究方法：
案例教学法代表着当代教育中比较新颖、颇有前景的一种教学方法，是指教师提炼和采用现实生活中已经发生的一些典型案例对原理、理论和道理进行解释，在教师的指导下运用多种方式启发学生独立思考，以加深学生对它们的理解和记忆。对将案例教学法应用于《数值计算方法》课程的主要研究方法为：
（1）调查分析法
采用调查分析的研究方法，在制定案例实施计划之前做好准备工作。运用座谈、问卷、网上留帖等手段，收集学生对案例背景材料的客观看法和疑问。通过对调查资料的整理分析，做到教师对案例材料的有的放矢，使学生得到正确引导。
（2）归纳总结法
对案例材料内容研究的过程和结果进行判断和评价。判断和评价其是否达到原实施计划要求，可以从授课现场气氛、学生课后反映、校内外数学建模比赛水平分析得出。进而归纳总结出案例内容的数学知识深度、采取的组织形式和下一步计划是否需要修正，以及修正的内容和方法。
（3）文献法
案例教学需要教师能够提炼和采用现实生活中已经发生的一些典型案例对计算方法中的原理、理论进行解释，因此教师可以通过阅读有关文献，掌握学科前沿，从工程实践中凝练科学问题，寻找新的思路。
（4）经验借鉴法
任课教师团队积极开展交流与学习，除了集体备课，对案例进行分析讨论和筛选，交流教学心得和体会外，还应向本校有教学经验的老教师探讨请教，并且还可以到校外参加教学研讨和交流，借鉴经验。
随着计算机软硬件技术的不断发展，科学计算与理论研究、科学实验一起构成了当今科学研究的三大方法，数值计算方法是科学计算的重要基础，也是理工科大学生和研究生必须掌握的核心课程。案例教学法强调自主探索和实践，从实际应用提炼科学问题，更加贴合数值计算方法的课程性质，充分激发学生的学习兴趣，提高了学生对数值计算方法思想精髓的理解水平和应用能力。案例教学法无疑也给教师提出了新的挑战，需要教师投入更多的精力，不再是孤立地在一门课程的理论教学上下苦功，更要融会贯通地开展全程化教学。在教学中必须认真思考并合理解决的问题需要在今后的教学实践中不断补充、丰富和完善。
参考文献
[1] 李庆扬，能超，易大义.数值分析（ 第 5 版）.清华大学出版社，2008年
[2] 刘春凤. 实用数值分析教程. 冶金工业出版社，2006年
[3] 杜廷松. 关于《数值分析》课程教学改革研究的综述和思考 [J]. 大学数学， 2007，23 .
数值有限元方法 篇12
关键词：隧道,有限元方法,数值模拟,强度安全系数
随着社会的发展,如何解决日益拥挤的交通问题成为城市建设中的重头戏。其中,公路隧道等地下工程成为解决这一难题的关键。现行隧道及地下工程结构设计,要求在经济合理的前提下满足安全性、适用性、耐久性,即满足结构的可靠性。因此以可靠度理论为基础的极限状态设计对隧道及地下工程更为必要。隧道衬砌的投资约为隧道全部投资的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,检查其是否满足规范要求
(1)式中:Nu为偏心受压构件的极限承载力;K为截面的抗压(或抗拉)强度安全系数;
当由抗压强度控制,即
Nu=φαRabh (2)
(2)式中: φ为构件纵向系数,隧道衬砌取1;Ra为混凝土极限抗压强度;b为截面宽度,通常取1 m;h为截面厚度。α为轴力的偏心影响系数,按以下经验公式确定:
当由抗拉强度控制,即
(4)式中:φ为构件纵向系数,隧道衬砌取1;R1为混凝土极限抗拉强度;b为截面宽度,通常取1 m;h为截面厚度。
采用以上公式,对隧道拱圈单元的安全系数进行了计算。根据以上数据列表统计了变厚度衬砌在拱脚、拱腰、拱顶单元的安全系数,并绘制了单元强度安全系数图。
在对不同厚度的衬砌安全系数的计算过程中出现的最小安全系数为4.751 , 发生在厚度为45 cm的处于拱腰的衬砌断面,即衬砌断面厚度为45 cm时的初期支护的受力是最不利的,但对照中国《公路隧道设计规范》中的规定,初期支护在该地质条件及施工条件下仍然是安全的。
根据上述内力结果计算衬砌结构安全系数。对衬砌的安全性能进行检验。根据《公路隧道设计规范》规定,混凝土偏心受压构件按破坏阶段进行强度验算。具体计算方法为根据材料的极限强度,计算出偏心受压构件的极限承载力,与实际内力相比较,得出截面的抗压(或抗拉)强度安全系数,检查是否满足规范要求[15]。最后绘出安全系数图(内层实线为初期支护中心线,外层虚线为安全系数)。衬砌支护的安全系数如图10所示。衬砌厚度变化下的结构安全系数如表3所示。
3 结 论
通过对二维有限元数值模拟分析,得出以下主要结论:
1)从衬砌结构的轴力、弯矩图可以看出,对于同一厚度的衬砌结构,边墙的轴力比较大,而边墙与仰拱连接处的弯矩和剪力都特别大,有应力集中现象。
2)在衬砌厚度变化的情况下,随着衬砌厚度增加,衬砌各项内力均呈上升趋势,从而造成衬砌强度安全系数增大,可见衬砌厚度是影响其强度安全系数的一个重要因素。
同时可以看出,用有限元数值方法对隧道开挖支护进行模拟能较真实地反映隧道结构的内力和变形情况,也可以用于对二次衬砌结构安全性能进行评价。
【数值有限元方法】推荐阅读：
有限元数值模拟11-12
有限元数值模拟仿真08-06
数值计算方法10-07
数值计算方法实验大纲12-10
数值计算方法结课心得07-06
一种新的模拟渗流运动的数值方法10-10
激光与靶相互作用体汽化模型的数值计算方法11-24
时域有限元方法07-02
有限元分析方法11-13
非线性有限元方法11-19