三维有限元法（精选十篇）
三维有限元法 篇1
国内外专家学者对盾构法施工的研究方法可归纳为:经验公式法、实测数据回归、室内模拟试验、数值模拟法等途径[1]。对隧道盾构法施工过程进行了三维数值模拟分析, 对盾构施工过程中隧道围岩的应力、位移和地表的沉降及管片结构受力情况进行了研究, 为以后的设计和施工提供相关依据。
1 工程概况
本工程为陕北某铁路专用线一长约200 m的隧道, 由于其进口位置地势起伏不大, 为使计算方便, 本文取进口位置长60 m进行分析。地层从上至下, 根据地层的物性参数不同将其分为3层, 第 (1) 层为黄土, 厚度约10 m, 第 (2) 层为强风化砂岩, 厚度约20 m, 第 (3) 层为弱风化砂岩, 厚度大于20 m。隧道埋深约15 m, 洞身位于第 (2) 层强风化砂岩中, 隧道衬砌内径为5.4 m, 衬砌厚度为0.6 m。各围岩的分布及主要物理、力学性质见表1。
2 数值模拟
本模型取宽为60 m, 长为60 m, 高为45 m进行模拟计算, 共计11 920个单元, 节点13 264个。单元采用ANSYS中的Solid45单元, 其模型及单元划分见图11。。
本次数值模拟的约束边界条件为:前后、左右施加水平方向的位移约束, 顶部为自由端, 下部为竖直方向的位移约束。
本模型采用全断面开挖, 每次开挖长度为6 m, 随后进行衬砌施工。在数值模拟分析过程中, 利用ANSYS软件中的生死单元, 模拟隧道在开挖过程中的真实受力情况。
3 计算结果及分析
在开挖前在重力作用下Y方向的位移见图2。图3~图5为施工0 m~6 m, 6 m~12 m, 12 m~60 m时各地层在Y方向的位移云图。从图中可以分析出:各地层竖向位移的最大值出现在隧道的拱顶和仰拱处, 并且随着隧道的开挖, 其竖向位移逐渐变大。当隧道贯通时, 其中拱顶产生最大下沉约0.343 m, 仰拱处产生最大上隆约0.343 m。从中可以看出盾构法比其他施工方法在隧道施工过程中的优势。
图6~图8为隧道施工0 m~6 m, 6 m~12 m, 12 m~60 m时地面沉降云图。从图中可以看出:隧道向前推进的过程中, 已开挖段的地表会产生细微的沉降, 其沉降量为0.002 m~0.025 m, 同时盾构机在向前推进时, 对前方土体产生一定的挤压作用, 使得前方一定范围内地表出现隆起的现象, 其隆起量可以忽略不计 (0.001 m~0.006 m) 。由此可见盾构法在浅埋隧道施工中对地面的影响非常小, 所以在许多地下工程中, 尤其是城市地铁施工, 盾构法的优势更加明显。
图9~图11为隧道施工0 m~6 m, 6 m~12 m, 12 m~60 m时衬砌的应力云图。从图中可以看出:衬砌应力最大值主要在左、右两侧, 说明衬砌主要受力是水平方向的力。洞口位置的应力值要比中间段大。
4 结语
通过对隧道盾构法施工的三维有限元分析, 得到以下结论:
1) 各地层竖向位移的最大值出现在隧道的拱顶和仰拱处, 并且随着隧道的开挖, 其竖向位移逐渐变大。
2) 盾构法在浅埋隧道施工过程中对地面的影响非常小, 所以在许多地下工程中, 尤其是城市地铁施工, 盾构法的优势更加明显。
3) 衬砌应力最大值主要在左、右两侧, 说明衬砌主要受力是水平方向的力。洞口位置的应力值要比中间段大。
摘要：根据有限元的基本原理, 对隧道盾构法施工过程进行了三维数值模拟分析, 研究了盾构施工推进过程中隧道围岩的应力、位移和地表的沉降及衬砌结构受力情况, 为以后的设计和施工提供相关依据。
关键词：ANSYS,有限元,数值模拟,隧道,盾构法
参考文献
三维有限元法 篇2
航空发动机转子叶片三维有限元振动特性分析
叶片是航空发动机的主要零件之一,结构及承载情况十分复杂.在实际使用中,由于叶片的振动破坏而造成发动机失效甚至飞机失事的例子时有发生.因此,对叶片的振动特性分析就显得尤为重要.根据叶片的实体结构建立了叶片的三维有限元模型,编制了叶片振动特性分析的有限元程序(Fortran语言,约8000余句).采用子空间迭代法计算叶片的频率.利用开发的`软件对影响压气机叶片振动特性的几个因素进行了分析.经算例考核证明,计算模型和程序是正确的、有效的,具有较高的计算精度.为转子叶片的疲劳损伤及可靠性分析,提供了一个有效的、灵活的研究手段,具有应用价值.
作 者：贺威 黄宝宗 HE Wei HUANG Bao-zong 作者单位：沈阳农业大学,高等职业技术学院,沈阳,110004刊 名：沈阳农业大学学报 ISTIC PKU英文刊名：JOURNAL OF SHENYANG AGRICULTURAL UNIVERSITY年，卷(期)：37(5)分类号：V232.4关键词：叶片 静频 动频 振动特性 有限元
大地电磁三维矢量有限元正演模拟 篇3
关键词：矢量有限元；大地电磁；正演模拟；张量阻抗
中图分类号：P631 文献标识码：A
文章编号：1674-2974（2016）10-0119-07
Abstract：Starting from the Maxwell equations， this article studied the boundary conditions of 3D MT. By using the weighted residual method， we derived the three-dimensional MT finite element equation. The three-dimensional vector finite element hexahedral meshing mode was introduced and the basis functions were selected. Then we derived the three-dimensional magnetotelluric vector finite element stiffness coefficient matrix and discrete format. A three-dimensional vector finite element magnetotelluric forward Matlab program was done. The apparent resistivity curve of the dimensional COMMEMI 3D-1 model matches the international standard test data， which proves the correctness of 3D magnetotelluric forward program. With the analysis of high and low resistivity anomalies， it shows that tensor impedance map can roughly determine the anomaly characteristics， which enriches the magnetotelluric response characteristics of expression.
Key words：vector finite element method；magnetotelluric； forward modeling； impedance tensor
大地电磁（MT）是以电离层激发的天然交变电磁场为场源，在地表观测相互正交的电场、磁场分量来获取地电构造信息的一种重要地球物理勘探方法[1].MT不需要庞大的发射源设备，只需采用比较轻便的接收设备，野外工作方便、成本低，被广泛应用于地壳和上地幔电性结构的研究，在石油天然气勘探、矿产资源勘探、工程与环境普查等领域，发挥着举足轻重的作用[2-9].可以预见，三维MT勘探技术是地球物理中深层领域的研究热点及今后MT的发展趋势，而三维MT正演是理解MT勘探物理现象并认识地质体电磁响应规律的有效手段，显然尤其重要.
尽管矢量FEM拥有诸多优点，但在地球物理的电磁法正演领域中，其应用并不多见，尚需要进一步完善.目前的研究主要包括：Yoshimura 等 [10]开展了矢量FEM的MT响应数值模拟，并将矢量FEM的计算结果与交错网格FDM的计算结果进行了对比；Mitsuhata等 [11]利用矢量FEM和节点FEM耦合的方法对三维MT数值模拟；Nam [12]采用不规则六面体矢量FEM直接计算电场，研究了起伏地形下MT的电阻率和相位的变化规律；刘长生等[13]将完全非结构化四面体单元引入到矢量有限元中，实现了三维大地电磁h-型自适应矢量有限元正演；王烨[14]开展了高频率大地电磁法矢量有限元正演，并采用改进的威尔金森方法求解大型病态方程组，提高了迭代速度；顾观文等[15]开展了矢量有限元法MT三维地形数值模拟，研究了地形起伏下三维阻抗张量的变化规律；杨军等[16]采用非结构四面体单元的三维矢量FEM实现了海洋可控源电磁数值模拟；苏晓波等[17]采用规则六面体单元的三维矢量FEM实现了大地电磁数值模拟，并对网格剖分的重要性进行了研究.
在前人基础上，作者推导了三维大地电磁矢量FEM正演的离散形式，应用矢量FEM算法计算了三维COMMEMI 3D-1国际模型[18]的MT视电阻率及模型张量阻抗，研究了高低阻异常体的电磁响应特性，有效地指导了MT的资料解释.
图5（a），（b）分别为f=0.1 Hz时XY模式与YX模式下矢量FEM正演视电阻率曲线.分析图5（a），（b）可知，两幅图中的矢量FEM曲线与COMMEMI所提供的数据都能够很好地吻合，说明无论是在低频还是高频部分，应用矢量FEM开展三维大地电磁正演，都具有较高的精度，同时也验证了矢量FEM算法及程序的正确性.
图6为应用矢量FEM正演计算COMMEMI3D-1模型得到的张量阻抗.由图可见，10 Hz与0.1 Hz两个频率下的张量阻抗形态基本一致，10 Hz的数值较0.1 Hz要大.对比图中4个不同的张量阻抗，可以发现，图6（a），（d），（e）和（h）中两个频率下的Zxx与Zyy分为四瓣，且阻抗值较小，四瓣的中心反映了异常体的边界，而图6（b），（c），（f）和（g）中的Zxy与Zyx阻抗值较大，反映了入射场的特性.根据张量阻抗理论可知，当构造为二维构造时，Zxx和Zyy为零，即当异常体走向方向越长，Zxy与Zyx越小，Zxy与Zyx差异也越大.由此，张量阻抗分解后，无需做反演即可以判断出异常体的简单特性.
为了进一步认识大地电磁的响应特性，对比高、低阻异常体张量阻抗的不同，在图3中COMMEMI3D-1测试模型的基础上，仅将低阻异常体改为1 000 Ω·m高阻异常体.其他参数均与国际模型相同.应用三维矢量FEM开展三维高阻异常体模型的张量阻抗研究.
图7为应用三维矢量FEM正演的10 Hz大地电磁张量阻抗图.分析图7（a）与图7（d）可知，高阻异常体张量阻抗中的Zxx与Zyy同样分为四瓣，且阻抗值较小，其四瓣的中心反映了异常体的边界.由于异常体x方向与y方向的比值为1∶2，图7（b）中的Zxy与图7（c）中的Zyx差异较大.对比高低阻异常10 Hz时的阻抗相位Zxx，虽然两者都为四瓣，但是阻抗值正负值的分布正好相反，低阻异常体四瓣的中心向外辐射，幅值变小趋于0；而高阻异常体四瓣的中心向外辐射，幅值变小趋于0之后会发生反转之后再次趋于0.Zyy具有相同的规律.对比Zxy和Zyx，高阻异常体中心仅出现一个闭合异常形态，而低阻异常体则形态更为复杂.
4 结 论
1） 介绍了三维矢量有限元区域剖分方式，对矢量FEM插值基函数以及单元插值方式进行了阐述，应用Galerkin算法，推导了三维矢量FEM大地电磁方程离散格式，编制了矢量FEM三维MT的Matlab模拟程序.
2） 设置三维COMMEMI 3D-1模型进行矢量FEM的计算，模拟结果与COMMEMI提供的数据拟合效果很好，验证了矢量有限元程序的正确性.通过对比高低阻异常体的张量阻抗，分析了不同异常下张量阻抗的特点，进一步认识了MT的响应特性.
参考文献
[1] 柳建新，童孝忠，郭荣文，等. 大地电磁测深勘探：资料处理反演与解释[M]. 北京：科学出版社，2012：1-12.
LIU jian-xin， TONG Xiao-zhong， GUO Rong-wen， et al. Magnetotelluric sounding exploration： data processing inversion and interpretation [M].Beijing： Science Press， 2012：1-12.（In Chinese）
[2] 底青云，王若. 可控源音频大地电磁数据正反演及方法应用[M]. 北京：科学出版社，2008：1-8.
DI Qing-yun， WANG Ruo. Controlled source audio magnetotelluric data inversion and methods application [M].Beijing： Science Press， 2008： 1-8.（In Chinese）
[3] 谭捍东，余钦范，JOHN B，等. 大地电磁法三维交错采样有限差分数值模拟[J]. 地球物理学报，2003，46（5）：705-711.
TAN Han-dong， YU Qin-fan， JOHN B， et al. Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method[J]. Chinese Journal of Geophysics， 2003， 46（5）：705-711.（In Chinese）
[4] 陈辉，邓居智，谭捍东，等. 大地电磁三维交错网格有限差分数值模拟中的散度校正方法研究[J]. 地球物理学报，2011，54（6）：1649-1659.
CHEN Hui， DENG Ju-zhi， TAN Han-dong， et，al. Study on divergence correction method in three-dimensionalmagnetotelluric modeling with staggered-grid finite difference method [J]. Chinese Journal of Geophysics， 2011， 54（6）： 1649-1659. （In Chinese）
[5] 李焱，胡祥云，杨文采，等. 大地电磁三维交错网格有限差分数值模拟的并行计算研究[J]. 地球物理学报，2012，55（12）：4036-4043.
LI Yan， HU Xiang-yun， YANG Wen-cai， et al. A study on parallel computation for 3D magnetotelluric modeling using the staggered-grid finite difference method [J]. Chinese Journal of Geophysics， 2012， 55（12）： 4036-4043. （In Chinese）
[6] 徐凌华，童孝忠，柳建新，等. 基于有限单元法的二维/三维大地电磁正演模拟策略[J]. 物探化探计算技术，2009，31（5）：421-425.
XU Ling-hua， TONG Xiao-zhong， LIU Jian-xin， et al. Solution strategies for 2D and 3D magnetotelluric forward modeling based on the finite element method [J]. Computing Techniques for Geophysical and Geochemical Exploration， 2009， 31 （5）： 421-425. （In Chinese）
[7] MOGI T. Three-dimensional modeling of magnetotelluric data using finite- element method [J]. Journal of Applied Geophysics， 1996， 35（2）： 185-189.
[8] 梁生贤，张胜业，吾守艾力，等. 复杂三维介质的大地电磁正演模拟[J]. 地球物理学进展，2012，27（5）：1981-1988.
LIANG Sheng-xian， ZHANG Sheng-ye， WU Shouaili， et al. Magnetotelluric forward modeling in complex three- dimensional media [J]. Progress in Geophys， 2012， 27（5）： 1981-1988. （In Chinese）
[9] 金建铭.电磁场有限元方法[M]. 西安：西安电子科技大学出版社， 1998：164-177.
JIN Jian-ming. Finite element method of electromagnetic field [M]. Xian： Xian Electronic University Press， 1998：164-177.（In Chinese）
[10]YOSHIMURA R，OSHIMAN N. Edge-based finite element approach to the simulation of geoelectromagnetic induction in a 3-D sphere [J]. Geophysical Research Letters， 2002， 29（3）： 1039-1042.
[11]MITSUHATA Y， UCHIDA T. 3D magnetotelluric modeling using the T-Ω finite-element method [J].Geophysics， 2004， 69（1）：108-119.
[12]NAM M J， KIM H J， SONG Y， et al. 3D magnetotelluric modelling including surface topography [J]. Geophysical Prospecting， 2007， 55（2）： 277-287.
[13]刘长生，汤井田，任政勇，等. 基于非结构化网格的三维大地电磁自适应矢量有限元模拟[J]. 中南大学学报：自然科学版，2010，41（5）：1855-1860.
LIU Chang-sheng， TANG Jing-tian， REN Zheng-yong， et al. Three-dimension magnetotellurics modeling by adaptive edgefinite-element using unstructured meshes [J]. Journal of Central South University： Science and Technology，2010， 41（5）： 1855-1860. （In Chinese）
[14]王烨. 基于矢量有限元的高频大地电磁法三维数值模拟[D]. 长沙：中南大学地球科学与信息物理学院，2008：35-62.
WANG Ye. A study of 3D high frequency magnetotelluric modeling by edge-based finite element method [D]. Changsha： Central South University. School of Geosciences and Info-Physics， 2008：35-62. （In Chinese）
[15]顾观文，吴文鹂，李桐林. 大地电磁场三维地形影响的矢量有限元数值模拟[J]. 吉林大学学报：地球科学版，2014，44（5）：1678-1686.
GU Guan-wen， WU Wen-li， LI Tong-lin. Modeling for the effect of magnetotelluric 3D topography based on the vector finite-element method [J]. Journal of Jilin University： Earth Science Edition， 2014，44（5）：1678-1686. （In Chinese）
[16]杨军，刘颖，吴小平. 海洋可控源电磁三维非结构矢量有限元数值模拟[J]. 地球物理学报，2015，58（8）：2827-2838.
YANG Jun， LIU Ying， WU Xiao-ping. 3D simulation of marine CSEM using vector finite element method on unstructured grids [J]. Chinese Journal of Geophysics， 2015， 58（8）： 2827-2838. （In Chinese）
[17]苏晓波，李桐林，朱成，等.大地电磁三维矢量有限元正演研究[J].地球物理学进展，2015，30（4）：1772-1778.
SU Xiao-bo， LI Tong-lin， ZHU Cheng， et al. Study of three-dimensional MT forward modeling using vector finite element method [J].Progress in Geophysics， 2015，30（4）：1772-1778. （In Chinese）
三维有限元法 篇4
随着光纤通信和光信号处理技术的发展,人们越来越关注光集成技术。目前已有许多数值分析方法分析集成光波导,包括线性法、模匹配法、耦合模理论、光束传输法(BPM)、时域有限差分法和有效折射率法等。由于BPM[1]简单方便,准确性高,因此BPM是目前光波导器件研究与设计领域最流行的方法,也是最适合于开发通用光波导器件CAD软件的方法。从模拟空间的角度, BPM可分为二维BPM和三维BPM。
对于一般的用SiO2制作的阵列波导光栅(AWG),我们可以采用有效折射率方法将三维情况转化为二维情况,再进行计算;而对于InP制作的AWG,用有效折射率方法进行计算就很不准确,此时最好用三维BPM来对光波导器件进行分析。
本文将介绍三维有限差分光束传输法(FD- BPM),并采用三维FD-BPM对光场在阵列波导中的传输过程进行计算,给出模拟结果。
1 三维FD-BPM
三维FD-BPM[2,3]是在亥姆赫兹方程的基础上进行推导。三维亥姆赫兹方程如式(1)所示:
将式(1)写成差分形式,得到式(2)和式(3),式中,i表示x方向;j表示y方向,此两式仍为三对角方程组。
合并相应项后,最后求解方程为式(4)和式(5)。
式中,
由于是三维,x、y方向都为未知数,要用方向交替隐含追赶法对其求解。这种方法的实质是把时间步长Δt 分成两个半步, 前半步在x方向用隐格式, y方向用显格式;后半步在y方向用隐格式, x方向上用显格式;如此循环,逐步求解。
2 三维FD-BPM在AWG模拟计算 中的应用
2.1 AWG器件的基本结构
AWG器件主要是由输入/输出波导、AWG及连接它们的两个对称平板波导组成,输入/输出波导连接平板波导的一端,均匀地排列在一个半径为R的罗兰圆上,每条阵列波导正对中心输入/输出波导出射端位置O,均匀地排列在以O为圆心、半径R1=2R的光栅圆周上,这种设计采用的即是罗兰圆结构设计。阵列波导是AWG器件的分光部分,相邻波导的光程差是中心载波波长的整数倍,阵列波导数为2M+1。AWG平面结构图如图1所示。
2.2 三维FD-BPM在AWG模拟计算中的应用
一般AWG模拟计算[4,5],对光场经过阵列波导部分采用的处理方法是直接在输入平板的输出光场上加上位相差,这样就没有考虑相邻阵列间的光场串扰作用。为了更准确地计算传输光场,本文对阵列波导部分的模拟计算采用的处理方法是先使用三维FD-BPM计算输入平板的输出光场在经过一段直波导阵列后的光场分布,然后再加上位相差。设计一个1×8 AWG,图2为用三维方法计算的1×8 AWG中输入平板波导的输出光场经过一段直波导阵列传输后的三维图形,可以看出光场被很好地限制在直阵列中传输,达到了稳定的传输状态。
图3是采用三维FD-BPM时1×8 AWG的光谱响应图。图4是没有计算光场在直阵列中传输时1×8 AWG的光谱响应图。比较图4与图3可以发现,图4所示光谱变窄了,光功率变小了,而且相邻信道光谱的交点降低了,不同信道的光谱相差很大,这是因为没有考虑阵列波导间的串扰。由此可见,加入对直阵列波导的光场计算提高了计算的精度,使计算结果更接近于实际情况。
3 对三维FD-BPM步长的讨论
由于三维FD-BPM中把时间步长Δt分成两个半步, 因此在计算同样的传播长度时,比二维FD-BPM要多花一倍的时间。但是通过多次改变步长计算发现,三维FD-BPM对纵向步长是不敏感的,纵向步长取1 μm和取5 μm时,阵列波导输出端各阵列中心的场强和相位没有很大区别。因此在用三维FD-BPM进行计算时,我们可以采用加大步长的方法来减少计算时间,使得在计算精度得到提高的情况下,计算时间没有受到很大影响。
4 结束语
在AWG的光场模拟计算中,本文采用三维FD-BPM实现光场在阵列波导中的传输模拟,这种方法加入了对直阵列波导的光场计算,提高了计算的精度,使计算结果更接近于实际情况。本文还对三维FD-BPM的步长进行了计算和讨论,计算结果表明,三维FD-BPM对纵向步长不敏感。因此在采用三维FD-BPM进行计算时,可采用加大步长的方法来减少计算时间,使得在计算精度提高的情况下,计算时间没有受到很大的影响。
参考文献
[1]刘辛,鲁平,刘德明.一种修正的有限差分光束传播方法[J].中国激光,2005,32(4):511-513.
[2]Fujisawa Takeshi,Koshiba Masanori.Full-vector fi-nite-element beam propagation method for three-di-mensional nonlinear optical waveguides[J].Journal ofLightwave Technology,2002,20(10):1 876-1 884.
[3]王子谦,祝宁华.一种三维有限差分光束传输法[J].中国激光,1994,21(11):909-913.
[4]刘辛,鲁平,刘德明.修正的光束传输法在阵列波导光栅中的应用[J].光子学报,2005,34(7):1 015-1 018.
三维有限元法 篇5
在分析济宁市水文地质条件的基础上,对引起地面沉降诸因素进行了分析,并建立了地面沉降量与地下水位变幅之间的`相关关系.集中过量开采地下水是引起济宁市地面沉降的主要原因.在此基础上,建立了准三维地下水流模型和一维地面沉降模型.通过水力联系建立地下水与地面沉降耦合数值模型,运用有限元法对地下水渗流场和地面沉降量进行模拟,并对和地面沉降进行了预测.
作 者：魏加华 崔亚莉 邵景力 陈占成 文唐章 作者单位：魏加华,崔亚莉,邵景力(中国地质大学水资源与环境工程系,北京 100083)
陈占成,文唐章(鲁南地质工程勘察院,山东六允州 27)
三维有限元法 篇6
【关键词】 下颌第一磨牙；牙体缺损修复；三维有限元模型
【中图分类号】R783.3 【文献标志码】 A 【文章编号】1007-8517（2016）04-0111-02
本次研究主要采用充填全冠修复和桩核全冠修复两种修复方法对牙体缺损进行修复，并采用金合金桩核全冠修复、镍铬合金全冠修复、复合树脂充填以及银汞合金这四种材料，观察使用不同材料和修复方法时，牙本质Mohr的应力分布情况，比较牙本质最大的Mohr值。针对下颌第一磨牙牙体缺损修复的三维有限元力学进行深入的分析。
1 资料与方法
1.1 建立三维有限元模型 对全口牙裂教具中右下颌第一磨牙进行三维激光扫描，获得数据，建立三维有限元基础模型。根据相关文献以及解剖学理论提供的根管壁和髓室的厚度、牙体和牙周组织之间的关系以及髓腔结构等，借助CAD软件，参考金属冠修复标准，建立三维几何模型。构建右下颌第一磨牙远中颌面的三维集合模型，借助有关软件划分网格，建立四节点四面体单元的三维有限元基础模型。以此为基础，在保证远中缺损龈壁和近中牙体壁不变的前提下，在髓室顶水平以上颊、舌向扩大洞型，建立舌侧壁和颌面颊厚厚度为1毫米的最大缺损模型，简称模型1。在髓室底平面向颊、舌向扩大洞型，当舌侧壁颌面端厚1mm的缺损时，颊侧壁为2.7mm，简称模型2。模型1和模型2都采用复合树脂和银汞合金充填，进行镍铬合金全冠修复。
1.2 组织和材料的力学参数 泊松比是指材料在单向受拉或受压时，横向正应变与轴向正应变的绝对值的比值，也叫横向变形系数，它是反映材料横向变形的弹性常数，而弹性系数是相互联系的两个指标增长速度的比率，它是衡量一个变量的增长幅度对另一个变量增长幅度的依存关系。组织和材料的力学参数如表1所示。
结果计算标准。根据Mohr强度理论公式进行计算，输出牙本质Mohr的应力云图[2]，公式为：σrm=σ1σtσcσ3≤[σ3]，[σt]指单向拉伸许用应力，取值40兆帕；[σc]指单向压缩许用应力，取值257兆帕；σ1是指最大主应力；σ3是指最小主应力。
1.3 观察指标 观察比较两种模型充填和桩核全冠修复在四种载荷下的牙本质Mohr应力分布规律和牙本质最大Mohr值。如果牙本质拉伸极限强度（40MPa）小于牙本质最大Mohr值，可以判断为牙体折裂[3]。
2 结果
2.1 牙本质最大Mohr值 在水平载荷、斜向载荷、垂直载荷三种不同的载荷下，模型1应优先选择银汞合金充填全冠修复和复合树脂充填全冠修复两种修复方法，如果这两种方法无法实施，可以根据实际情况选择金合金桩全冠修复或者镍铬合金桩核全冠修复方法。在垂直载荷和斜向载荷两种不同的载荷下，模型2应优先选择金合金桩核全冠修复和镍铬合金桩核全冠修复，其次选择复合树脂或者银汞合金充填全冠修复。
2.2 牙本质Mohr应力云图结果比较 充填全冠修复和桩核全冠修复比较，由充填全冠修复将应力集中区从牙颈部和冠边缘逐渐转移到根部。两种充填全冠修复方式相比，复合树脂充填全冠修复应力分布小于银汞合金充填全冠修复。最大载荷下，模型1中，4种修复颌面观Mohr应力云图如图1所示。模型1和模型2相比，模型1（max1）在平髓室顶牙体组织处容易产生应力集中区，模型2（max2）在平髓室底牙体组织处容易产生应力集中区。在图1中，暖色代表数值较大，冷色代表数值较小。
3 讨论
近年来，随着医学技术的不断成熟，磨牙根管治疗手术的水平也有了很大的提高，大大提高了下颌龋坏第一磨牙的保存率。临床修复时，如果牙体缺损并没有到髓室壁，应该尽量选择全冠修复，因为髓室壁牙体组织的保留有助于减少修复后牙体组织的应力[4]。如果髓室壁牙体组织受到部分缺损，应该选择桩核全冠修复方法，并且尽量使用金合金材料，比镍铬合金材料的效果更好。如果受到外部环境因素的影响，不具有进行桩核全冠修复的条件，应尽量采用复合树脂作为充填材料，相对于银汞合金充填而言，复合树脂充填效果更好。如果水平载荷非常大，或者患者有磨牙的习惯，此时，无论采用何种修复方法，在进行修复之前，都要减轻颌力，才可以保证修复不受影响[5]。同时，如果髓室壁牙体组织遭到一定程度的破坏，尽量不要使用以镍铬合金为材料的桩核全冠修复方法。在本次研究中，斜向载荷、垂直载荷、最大载荷下，两类修复牙本质最大Mohr值都低于牙本质拉伸极限强度，说明在咀嚼能力正常的情况下，这两种修复方法都可以使用。两种充填全冠修复在牙本质最大Mohr值方面没有较大的差异，相对而言，在复合树脂充填全冠修复方法下，应力集中面积要小于银汞合金充填全冠修复。如果发现髓室壁牙体组织遭到非常严重的破坏，应力主要在髓室底集中，应尽量使用复合树脂充填全冠修复方法。金合金桩全冠修复和镍铬合金桩核全冠修复两种修复方法相比，镍铬合金桩核全冠修复的牙本质最大Mohr值更大，所以应尽量采用金合金桩全冠修复方法。
参考文献
[1]田力丽，梁伟，李凌旻，等. 下颌第一磨牙牙体缺损修复后静力与冲击动力的三维有限元分析[J]. 华西口腔医学杂志，2007，06：595-598-602.
[2]田力丽，梁伟，李凌旻，等. 纤维桩与金属桩核修复磨牙牙体缺损的三维有限元应力分析[J]. 口腔颌面修复学杂志，2008，04：280-282-290.
[3]刘涛，耿海霞，刘建彰，等. 下颌第一磨牙任意两壁缺损不同充填材料全瓷冠修复的三维有限元应力分析[J]. 临床口腔医学杂志，2015，02：76-79.
[4]连克乾，王安训，许鸿生，等. 三维有限元法分析带桩嵌体修复下颌第一磨牙残冠的等效应力分布[J]. 中国临床康复，2005，30：112-114.
[5]徐红梅，贾静，朱晓英，等. 根管治疗后下颌磨牙采用各种充填材料修复的三维有限元分析[J]. 口腔颌面修复学杂志，2014，02：75-79.
三维有限元法 篇7
1 材料和方法
1.1 材料与设备
1.1.1 建立下颌骨颏部正中骨折模型
利用Ansys软件建立颏部正中骨折三维有限元模型并将骨折线的宽度设定为1 mm[2],将骨折线上的单元属性设定为肉芽属性,从而模拟骨折初期情况,并在下颌骨骨折模型上建立坚强内固定系统(图 1~5)。
1.1.2 坚强内固定系统
4 孔小型接骨板(中国西安中邦公司),长、宽、厚分别为28.0、 4.0、 1.0 mm。固定螺钉长、直径分别为7.0、 2.0 mm。
1.2 接骨板放置的位置
1.2.1 颏部二维双板
①下缘接骨板位于下颌骨下缘上方约3 mm处,分别在其上方3 mm(距下颌骨下缘约10 mm)处、 5 mm(距下颌骨下缘约12 mm)处、7 mm(距下颌骨下缘约14 mm)处安放接骨板, 2 板平行,模拟3 种固定方式。
A:下颌骨颏部正中骨折; B:常规二维双板固定; C:三维双板固定; D:常规二维三板固定; E:三维三板固定
1.2.2 颏部三维双板
①下颌骨下缘下方安放接骨板,分别在其上方距下颌骨下缘约10 mm 处、12 mm 处、14 mm 处安放接骨板模拟3种固定方式。
1.2.3 颏部二维与三维三板
①下缘接骨板位于下颌骨下缘上方约3 mm处,在其上方4 mm(距下颌骨下缘约11 mm)和7 mm(距下颌骨下缘约14 mm)处安放2 块接骨板, 3 板平行。②下颌骨下缘下方安放接骨板,在其上方距下颌骨下缘约11 mm和14 mm处安放2 块接骨板,此2 板平行。
1.3 试验条件假设
1.3.1 材料力学参数
根据参考文献设定材料属性,其中皮质骨、松质骨设为各向异性,牙齿设为各向同性。具体材料力学参数如下。
皮质骨[3]:
E1=E2=13 GPa V12=V21=0.22
E3=19 GPa V31=V32=0.42
G12=5.3 GPa V23=V13=0.29
G23=G31=5.9 GPa
松质骨[3]:
E1=E2=273 MPa V12=V21=0.19
E3=823 MPa V31=V32=0.335
G12=115 MPa V23=V13=0.105
G23=G31=5.9 GPa
牙齿[4]:E=20 290 MPa V=0.30
肉芽组织[2](骨折断层愈合初期属性):
肉芽组织的强度极限[5] σb=0.047 5 MPa
钛板及螺钉[6]:E=103 400 MPa V=0.35
钛的屈服极限[6]: σs=200~260 MPa
(注:E 杨氏弹性模量,G 剪切弹性模量,V 泊松比,ób 强度极限,ós屈服极限)。
1.3.2 模型的边界约束条件和载荷条件
将颞下颌关节处(髁突顶点)以铰链支撑处理;对下颌牙功能尖的节点进行约束,用以模拟各种咬合情况,对于不同咬合状态下施加相应的咀嚼肌力,并按坐标轴方向予以分解,在各种工况下加载。正常人日常咀嚼食物所需 力为29.4~294 N,约为最大力的一半,故本研究中正中咬合时肌力值取其1/2。其他咬合状态,肌力取可以相应产生100 N的平均咬合力的数值。模拟4 种咬合状态:正中咬合、前牙咬合、单侧磨牙咬合、单侧前磨牙咬合。翼外肌对分析闭颌运动时下颌骨的应力分布影响较小,可略去以简化计算[7]。
1.4 计算方法及计算内容
在下颌骨骨折坚强内固定模型中,比较颏部骨折时,不同固定方式、不同咬合状态下及不同咬合力下骨折段位移、骨断层的受力情况及接骨板的受力情况。
2 结 果
2.1 颏部正中骨折时骨折段的位移
从表 1~2得知:①颏部正中骨折时以上2 种内固定方法,骨折段位移均未超过骨折愈合所需的条件(小于0.15 mm)。②在各种咬合工况下, 2 板间距离越大,骨折段位移越小。③在各种咬合工况下,不同的固定方法,骨折段在正中咬合时位移最大。④将下缘钛板固定于下颌骨下缘下方骨面时在上板位置相同时骨折段位移小于平行双板。
2.2 颏部正中骨折时骨折断层应力分布
从表 3~4得知:① 颏部正中骨折时以上2 种内固定方法在各种咬合工况下,骨断层最大Von- Mises应力都未超过肉芽组织的强度极限(σb=0.047 5 MPa)。② 在各种咬合工况下,骨断层最大Von- Mises应力从大到小排列顺序均为: 2 板间距离越大, 骨断层最大Von- Mises应力反而越小。③ 在各种咬合工况下,常 规 双 板固定时 骨折 段在 正中咬 合时Von- Mises应力最大,而三维双板固定时骨折段在前牙咬合时Von- Mises应力最大。
2.3 颏部正中骨折时板钉应力
从表 5~6得知:颏部正中骨折时以上2 种内固定方法在各种咬合工况下,接骨板上应力都远小于接骨板的屈服极限(钛的屈服极限σs=900~1 000 MPa),应力集中区域为接骨板中部。
2.4 颏部正中骨折三板固定时骨折段位移
从表 7得 知: ① 2 种 固 定 方式在各 种 咬 合工 况
下,骨折段位移均未超过了骨折愈合所需的条件(小于0.15 mm)。②在各种咬合工况下,平行三板固定后的骨折段位移大于将下缘钛板固定于下颌骨下缘下方骨面的三维固定。③在各种咬合工况下,骨折段在正中咬合时位移最大。④增加板钉数目可明显减小骨折段位移且将下缘钛板固定于下颌骨下缘的三维固定优于平行固定。
2.5 颏部正中联合骨折三板固定时骨折断层应力
从表 8得知:颏部正中骨折时固定方式在各种咬合工况下,骨断层最大Von- Mises应力都未超过肉芽组织的强度极限(σb=0.047 5 MPa)。
2.6 颏部正中骨折三板固定时板钉应力
从表9 得知: 2 种固定方式在各种咬合工况下,接骨板上应力都远小于接骨板的屈服极限(钛的屈服极限σs=900~1 000 MPa)。
2.7 不同固定方式在不同力作用下骨折段位移情况
我们选取二维双板中固定效果最好的(双板相距7 mm,距下缘14 mm)和三维双板中固 定 效 果 最 好的
(上板距下缘14 mm),并与二维三板和三维三板一起进行比较在颏部正中骨折正中咬合时不同力作用下骨断层位移,见表 10。 从表 10得知:①二维平行固定方法无论双板和三板当正中咬合力小于190 N时骨折段位移均小于0.15 mm,未超过骨折愈合所需的条件。②三维固定方法无论双板和三板当正中咬合力小于210 N时骨折段位移均小于0.15 mm,未超过骨折愈合所需的条件。
3 讨 论
小型钛板沿张力带坚强内固定技术是目前治疗颌面部创伤与修复重建领域的主要方法,但在临床工作中发现小型钛板存在强度不够、稳定性不足[1,8]及抗扭力强度差等缺点。
对于下颌骨颏部正中骨折来说传统板钉只是在下颌骨侧面固定,这种固定位置可以达到二维空间的有效固定,防止骨折段受力时产生的弯矩作用和剪切力作用,但不能防止骨折段受力时产生的扭矩,而在下颌骨颏部正中骨折时骨折段有明显的扭力[9]。
本研究在上缘板钉位置相对固定的情况下对比常规固定和将板钉置于下颌骨下缘下方的三维固定,在不同骨折部位、不同咬合情况下的骨折段位移情况、骨折段应力和板钉应力。通过分析所得数据从理论上探讨这种固定方法的可行性和效果。
研究结果显示颏部常规平行双板固定方式和将下缘钛板固定于下颌骨下缘下方的三维固定方式在各种咬合工况下,骨折段位移、应力均未超过了骨折愈合所需的条件,在板钉应力方面在所有工况下,接骨板上应力都远小于接骨板的屈服极限,因此板钉在各种工况下均不会发生折断。在各种咬合工况下,不同的固定方法,骨折段均在正中咬合时位移最大,并且双板相距距离越大骨折段位移越小,因此当下颌骨颏部正中骨折时无论采取哪种固定方式术后都应避免正中咬合,并且在保证不损伤牙根及其神经的前提下,无论采取那种固定方式固定装置都应尽可能地拉开距离,这与其它学者的研究结果相符合[10]。通过比较可知将下缘钛板固定于下颌骨下缘下方时在上板位置相同、在各种咬合工况下骨折段位移均小于平行双板,因此下颌骨颏部正中骨折时,采用将下缘钛板固定于下颌骨下缘下方的三维固定方法是一种可靠的内固定方法,优于传统固定方法。
治疗下颌骨骨折可采用一点或多点固定,通常是2 点固定。我们在颏部骨折模型上试图增加1 点形成3 点固定,并比较其与2 点固定在各种咬合工况下,骨折段位移、应力的区别。我们发现平行三板固定在各种咬合工况下骨折段位移和应力等方面都优于平行双板,但次于三维双板固定中的上板距下缘12、 14 mm和三维三板固定。并且在各种咬合工况下三维三板和三维双板中上板距下缘14 mm时在骨折段位移和应力等方面差距并不明显。
下颌骨颏部正中骨折时通过改变正中咬合力的大小计算出不同固定方法和不同板钉数目的极限咬合力。二维平行固定方法无论双板和三板当正中咬合力大于190 N时骨折段位移均大于0.15 mm,超过了骨折愈合所需的条件,而三维固定方法无论双板和三板当正中咬合力小于210 N时骨折段位移均小于0.15 mm,未超过了骨折愈合所需的条件。因此,从不同载荷方面也能证明将下缘钛板固定于下颌骨下缘下方的三维固定方法是一种可靠的内固定方法,优于传统固定方法。
有研究表明骨段不稳固是导致术后感染的重要原因[11]。本研究表明将下缘钛板固定于下颌骨下缘下方骨面的三维固定骨折段位移小,而位移越小固定越稳固、越有利于骨折的愈合减少术后并发症的发生。
因此我们认为在下颌骨颏部正中骨折将下缘钛板固定于下颌骨下缘下方骨面的三维固定方式是安全可靠的,并且在各种工况下骨折段位移、应力等方面均优于常规平行板钉固定方式,可有效的防止因扭矩作用使骨折段产生的移位,使固定更为稳定,更有利于骨折的一期愈合。
虽然下颌骨三维固定方法在临床实践中[12]已证实其可行性,但缺乏生物力学的研究。本研究利用有限元法从生物力学角度对将下缘钛板固定于下颌骨下缘下方的三维固定方式进行定量分析,为临床治疗方案提供理论依据。
参考文献
[1]Ellis E,Walker L.Treatment of mandibular angle fracturesusing two nonconpression plate[J].J Oral Maxillofac Surg,19945,2:1032-1036.
[2]卢军,斯方杰,李明,等.下颌骨坚强内固定应力遮挡作用的有限元力学分析[J].口腔颌面外科杂志,1995,5(1):23-26.
[3]Hart RT,Hennebel VV,Thongpreda N,et al.Modeling thebiomechanics of the mandible:A three-dimensional finiteelement study[J].J Biomech1,9922,5(3):261-266.
[4]Batenberg RH,Raghoebar GM.Madibular overdentures sup-ported by two or four endosteal implant.A prospectivec,om-parative study[J].Int J Oral Maxillo fac Surg,1998,27(6):435-441.
[5]Hahn M,Nassutt R,Delling G,et al.The influence of ma-terial and design features on the mechanical properties oftranspedicular spinal fixation implants[J].J Biomed MaterRes,2002,63(3):354-362.
[6]Wagner A,Krach W,Schicho K,et al.A 3-dimensional fi-nite-element analysis investigating the biomechanical beha-vior of the mandible and plate osteosynthesis in cases of frac-tures of the condylar process[J].Oral Surg Oral Med OralPathol Oral Radiol Endod,2002,94(6):678-683.
[7]周学军,赵志河,赵美英,等.包括下颌骨的颞下颌关节三维有限元模型的建立[J].实用口腔医学杂志,2000,16(1):17-19.
[8]高学琴,李新明.坚强内固定联合颌间固定治疗下颌骨骨折96例分析[J].中国误诊学杂志,2007,7(15):3586-3587.
[9]于子莹,刘春丽.下颌骨骨折三维力学计算模型研究[J].中国实医药,2008,3(23):172.
[10]Edward Ellis.下颌骨骨折的固定处理[J].中华口腔医学杂志,2007,42(5):259.
[11]Mosbah MR,Oloyede D,Koppel DA,et al.Miniplate re-moval in trauma and orthognathic surgery a retrospectivestudy[J].Int J Oral Maxillofac Sury,2003,32(2):148-151.
桑郎拱坝三维有限元应力分析 篇8
关键词：应力,有限元,拱坝,等效应力法,拱梁分载法
1 概 述
桑郎水库工程是集发电、灌溉、供水为一体的综合利用工程, 灌溉面积1 266.7 hm2, 电站装机容量1.26万kW, 水厂供水水源2万t/d。坝址位于贵州省望谟县桑郎镇上游4km的“V”型狭谷河段, 控制集水面积650 km2, 多年平均流量10.7 m3/s, 水库正常蓄水位505 m, 最高校核洪水位505.36 m, 总库容1 480万m3。坝型为碾压混凝土拱坝。碾压混凝土拱坝为变圆心双曲拱坝结构, 坝体剖面顶宽5m, 坝顶高程507.0 m, 底宽17.5 m, 坝高90 m, 属高坝;其厚高比为0.194, 属薄拱坝。坝顶弧长123.115 m, 坝顶中心角92.815°。采用坝顶溢流, 溢流段布置在坝中部, 堰顶高程500 m, 设3扇5 m (高) ×8 m (宽) 弧形工作闸门, 溢流净宽24 m, 最大泄洪流量596 m3/s, 单宽流量24.8 m3/s。
桑郎拱坝属于薄拱坝, 按照规范可以采用拱梁分载法。但是拱梁分载法无法考虑坝身孔口对坝体应力分布的影响, 对于地基及两岸坝基的复杂性对坝体应力的影响亦不能准确计算, 所以对桑郎拱坝应进行三维有限元应力计算。利用有限元对坝体应力进行计算时, 在坝身溢流段及坝体角缘位置不可避免的会出现一定程度的应力集中, 对出现的应力集中应采用“有限元等效应力”法进行修正处理。
本文运用三维有限元法建立拱坝和地基的整体模型, 模拟拱坝逐级加载的施工过程, 计算分析拱坝在各种荷载组合的不同工况下的应力。并采用“有限元等效应力”法进行处理, 消除线弹性有限元计算出的拱坝角缘附近的应力集中, 使其应力状况满足大坝有限元应力控制标准。最后将拱坝角缘处用“有限元等效应力”求得的应力值和拱梁分载法的计算结果进行比较, 验证了有限元结果的准确性。
2 有限元等效应力法
等效内力为结点径向对应的单位单元上的内力, 单位单元定义为拱圈中心线宽度为1 的梁截面, 单位高度为拱径向拱截面及上下游面共同组成的体单元。作用在该单元上的内力通过局部坐标下的有限元应力沿厚度积分而得。有限元等效应力具体计算步骤为:
(1) 设拱坝的整体坐标系为 (x′, y′, z′) , 计算坝体应力的梁拱交点的局部坐标系为 (x, y, z) 。将有限元法计算的整体坐标系中的应力{σ′}=[σ′x, σ′y, σ′z, τ′xy, τ′yz, τ′zx]T, 经坐标变换, 得到局部坐标系的应力{σ}=[σx, σy, σz, τxy, τyz, τzx]T;
(2) 沿单位高度拱的径向截面和在中心线 (或坝轴线) 取单位宽度的水平梁截面, 对{σ}的有关分量进行积分, 得到拱圈和梁的内力 (包括梁的竖向力、切向剪力、径向剪力、弯矩、扭矩和拱的轴向力、径向剪力、弯矩) ;
(3) 按材料力学方法计算坝体应力。
3 有限元计算
3.1 荷载组合及材料参数
本文计算了拱坝在3种荷载组合下的应力状况, 按照《混凝土拱坝设计规范》 (SL282-2003) 的规定, 混凝土拱坝3种设计荷载组合见表 1。
坝体及地基岩体材料均采用线弹性材料, 材料参数见表 2。
3.2 有限元计算模型
参照同类型工程的经验, 有限元模型范围取为:大坝上游、左右岸拱端及底部地基取1倍坝高, 下游后取1.5倍坝高。在整个计算域内, 四周施加法向位移约束, 底部施加全部位移约束。
按坝体的实际体型进行建模, 并大致考虑了溢流堰的开口形式和闸墩的布置, 但不考虑具体的溢流堰形状, 采用8节点6面体等参单元对坝体及基础进行有限元离散。三维有限元整体网格及坝体细部模型如图 1。单元总数共计34 720个, 其中坝体单元7 360个, 节点总数共计42 041个。
3.3 计算结果
分别给出拱坝在3种荷载组合下的应力分布, 找出最大主应力及其出现的位置, 然后采用有限元等效应力法对其进行处理, 比较处理前后最大主应力的值, 并与拱梁分载法所得结果进行比较, 验证应力状况是否满足大坝有限元应力控制标准。应力等值线图的单位为Pa, 压应力为负, 拉应力为正。
(1) 基本组合下主应力分布。
由图 2、图 3第一主应力等值线图和第三主应力等值线图可知, 上游坝面受拉区主要分布在左右两坝肩区域, 上游坝面第三主应力完全表现为受压;下游坝面第一主应力大部分为拉应力, 受压区主要分布在左右两坝肩区域, 下游坝面第三主应力完全表现为受压。
坝体高拉应力分布于两个区域:一个是462～480 m高程上游坝面左右两坝肩处, 其分布还是较大的, 但绝大部分区域的拉应力仍然是低于0.7 MPa的, 仅在近上游两拱端的微小区域内出现了高于1 MPa的拉应力;另一个是下游坝面的中心附近, 但高于1 MPa的拉应力的范围也不大。高受压区域主要出现在下游坝面左右两坝肩及坝趾处, 但真正达到5.0 MPa的区域也仅限于下游两拱端的微小区域及坝趾区域内。相对整个坝体而言, 无论是高拉应力区域还是高压应力区域, 其区域都是较小的。
下游面溢流堰顶区域由于有尖角而发生拉应力集中, 最大拉应力达3.46 MPa;由于该处为钢筋混凝土结构, 拉应力由钢筋承担, 对大坝整体应力并无影响。
(2) 特殊组合一下主应力分布。
由图 4、图 5第一主应力等值线图和第三主应力等值线图可知, 上游坝面第一主应力基本表现为受压, 受拉区主要分布在左右两坝肩区域, 上游坝面第三主应力完全表现为受压;下游坝面第一主应力大部分为拉应力, 受压区主要分布在左右两坝肩区域, 下游坝面第三主应力完全表现为受压。
拉应力分布在462～480 m高程上游坝面左右两坝肩处, 绝大部分区域的拉应力仍然是低于0.5 MPa的, 仅在近上游两拱端的微小区域内出现了高于0.7 MPa的拉应力, 整体水平较低。高受压区域主要出现在下游坝面左右两坝肩及坝趾处, 但真正达到5.0 MPa的区域也仅限于下游两拱端的微小区域及坝趾区域内。相对整个坝体而言, 无论是高拉应力区域还是高压应力区域, 其区域都是较小的。
(3) 特殊组合二下主应力分布。
由图 6、图 7第一主应力等值线图和第三主应力等值线图可知, 上游坝面第一主应力基本表现为受压, 受拉区主要分布在左右两坝肩区域, 上游坝面第三主应力完全表现为受压;下游坝面第一主应力大部分为拉应力, 受压区主要分布在左右两坝肩区域, 下游坝面第三主应力完全表现为受压。
高拉应力分布在两个区域:一个是462～480 m高程上游坝面左右两坝肩处, 其分布还是较大, 但绝大部分区域的拉应力仍然是低于0.7 MPa, 仅在近上游两拱端的微小区域内出现了高于1 MPa的拉应力;另一个是下游坝面的中心附近, 但高于1 MPa的拉应力的范围也不大。高受压区域主要出现在下游坝面左右两坝肩及坝趾处, 但真正达到5.0 MPa的区域也仅限于下游两拱端的微小区域及坝趾区域内。相对整个坝体而言, 无论是高拉应力区域还是高压应力区域, 其区域都较小。
下游面溢流堰顶区域由于有尖角而发生拉应力集中, 最大拉应力达3.40 MPa;由于该处为钢筋混凝土结构, 拉应力由钢筋承担, 对大坝整体应力并无影响。
(4) 拱梁分载法和有限元等效应力计算成果比较。
为进一步验证有限元计算的成果, 现将拱坝角缘处用“有限元等效应力”求得的应力值和拱梁分载法的计算结果进行比较。拱梁分载法的计算结果来自《望谟桑朗水库工程初步设计报告》。表3为3种荷载组合下有限元计算所得的最大主拉应力和最大主压应力等效前后的应力值。取基本组合和特殊组合一下拱梁分载法与有限元等效应力法所得坝体角缘的主应力进行比较, 如表4～7为3种工况下拱梁分载法与有限元等效应力法所得坝体角缘的主应力比较表。
显然, 通过有限元等效应力的计算, 其拉应力及压应力的集中现象有了较大的消除, 无论是基本组合还是特殊组合下的最大主拉应力及最大主压应力均是满足大坝有限元应力控制标准的。
由以上比较成果分析可知, “有限元等效应力”求得的应力值和拱梁分载法计算的应力沿坝高的分布规律大致相近, “有限元等效应力”求得的主拉应力值一般稍大一些, 主压应力值基本相当, 符合一般规律。
4 结 语
本文完成了一拱坝在3种荷载组合下的应力状况计算分析, 主要分析结论可归纳如下。
(1) 有限元等效应力法是拱坝设计规范新规定的坝体强度校核方法, 此方法的应用将有助于发挥有限元计算能反映坝体开孔、体形局部变化和复杂地基对坝体应力影响的优势。
(2) 经过上述分析, 在有限元计算中, 其应力集中现象是难以避免的, 基本组合下的最大有限元主拉应力及最大有限元主压应力分别达到了1.67 MPa和6.0 MPa, 特殊组合一下的最大有限元主拉应力及最大有限元主压应力也分别达到了0.75 MPa和7.8 MPa, 特殊组合二下的最大有限元主拉应力及最大有限元主压应力也分别达到了2.0 MPa和6.1 MPa。基本组合、特殊组合一及特殊组合二下的最大主拉应力所对应的等效应力分别为1.27 MPa、0.53 MPa及1.44 MPa, 是满足有限元应力控制标准的, 而三种荷载组合下的有限元最大主压应力所对应的等效应力分别为3.42 MPa 、4.59 MPa及3.57 MPa, 也是满足有限元应力控制标准的。
参考文献
[1]SL282-2003, 混凝土拱坝设计规范[S].
[3]肖伟荣, 苏志敏, 唐涛.有限元等效应立法在拱坝设计中的应用[J].云南水力发电, 2005, 21 (1) , 36-39.
[4]李守义, 周伟, 苏礼邦, 等.基于ANSYS的拱坝等效应力研究[J].水利发电学报, 2007, 26 (5) , 38-41.
三维有限元法 篇9
管片是盾构隧道管片环的主体结构。从施工到正常使用的过程中,管片一般承受的荷载包括主要荷载(竖向与水平土压力、水压力、自重、超载重和周围土层的弹性抗力)、次要荷载(内部荷载、施工期荷载和地震效应)和特殊荷载(相邻隧道的影响、地基沉降的影响和其他荷载)。其中施工荷载主要是盾构千斤顶顶力、注浆压力和管片举重臂操作的荷载[1]。
以广州地铁为例,目前已完工的盾构隧道管片含钢量为(128～165) kg/m3,采用不同的含钢量,将会使工程投资有较大的区别[2]。因此,研究管片在施工阶段及正常使用阶段的力学特性,即应力、裂缝分布及薄弱部分,为配筋构造方式的优化提供依据,从节约工程投资及优化管片设计来说均有明显意义。
足尺试验是研究管片力学特性的最好方法,但要在试验中要模拟管片所受的复杂的荷载模式、量值及其边界条件存在一定困难。随着有限元理论与计算机技术的不断发展,采用有限元方法分析隧道衬砌受力已经成为一种适应性很强的方法。利用有限元法分析管片,能对管片施加相应于不同阶段的形式及大小的荷载,利用边界条件的变化模拟施工中出现的结构面不平整等情况。现利用三维有限元法,以广州地铁使用的管片为研究对象,使用通用有限元软件ADINA,分析了单块管片在各种荷载下的力学特性,以期为管片设计及配筋形式的优化提供一定依据。
1有限元分析模型的建立
1.1管片尺寸
广州地铁采用的是五加一拼装方式,即三块标准块管片(A1～A3),两块邻接块管片(B、C)和一块封顶块(也称为楔形块)管片(K),管片外直径6 m,内直径5.4 m,管片宽度1.5 m,厚度0.3 m。各环之间采用错缝形式拼装,相互错开±18°,标准衬砌环管片具体布置及相应编号见图1,错缝拼装形式见图2,标准衬砌环各管片的角度尺寸见图3。管片间及管片环间采用弯曲螺栓连接,组成错缝拼装的盾构隧道。
1.2混凝土本构模型及主要参数
ADINA软件的混凝土模型有三个基本特性[3]:①利用非线性应力-应变关系反映材料在上升压应力下的软化特性;②用破坏包络线定义材料在拉应力下开裂及压应力下压碎的特性;③能模拟材料在开裂及压碎后的特性。在分析结果中,程序将利用三种不同的标识表示三种不同的状态:张开裂缝、闭合裂缝和压碎。
管片采用C50混凝土,ADINA程序需要输入的参数及其量值如下:泊松比0.20,弹性模量3.45×107 kPa,极限抗拉强度6 072 kPa,极限抗压强度60 720 kPa,极限抗压强度对应的应变-0.002 028,残余抗压强度-30 360 kPa,残余抗压强度对应的应变-0.005 075。
1.3荷载、边界条件及计算模型
在施工阶段,管片环缝面承受的千斤顶顶力通常为1 600 kN/m[4],标准块管片弧长为3.4 m,侧面有8块传力垫(图4),则每块传力垫承受顶力为680 kN。计算中将千斤顶顶力设定为一般情况下千斤顶顶力的2倍,即每块传力垫1 360 kN,以考虑由于盾构机转弯或偶然千斤顶顶力不均匀的情况。
采用接触算法,在管片环缝面与传力垫之间建立接触面,模拟管片传力垫之间的相互接触。管片纵缝面既不是铰支端,也不是固定端。为准确地模拟相邻管片纵缝面的复杂边界条件,在计算模型中,在管片纵缝面设置了支承台用以模拟相邻的管片[5],支承台面与管片纵缝面之间采用接触面建立相互作用关系,使管片纵缝面能在不同荷载作用下能相互挤压或相互分离(图5)。
在正常使用阶段,每块管片均受到邻接两块管片及周围土层的作用,包括土压力、水压力和周围土层的弹性抗力。周围土层的弹性抗力是管片环受荷变形后,衬砌结构向周围土层相对移动而产生,实际上反映为管片受到的被动土压力,形式上也表现为管片外弧面变化的压力荷载。各种外荷载的分布形式为抛物线[6]。每块管片与其他管片的相互作用也使用支承台进行模拟(图6)。
混凝土管片、支座和传力垫按照文献[7]的图纸建模,模型各部分均采用三维实体单元模拟。为剔除其他因素的影响,充分研究管片本体的力学特性,同时参考文献[5,7],把单块管片在施工阶段及正常使用阶段受到的各种荷载概化为如表1所列5种荷载,共10种计算模型。
2计算结果及分析
2.1计算结果
计算结果见表2,各计算模型开裂形式见图7。在图7中,所有裂缝均为张开裂缝,圆环所在处均表示管片张开裂缝的位置
模型1～模型10的有限元网格如图6所示。
2.2计算结果的讨论
通过对单块管片三维有限元分析,可知管片有如下力学特性。
(1)在千斤顶顶力作用下,当管片环缝面不平整导致传力垫与管片间有间隙时,管片出现沿盾构推进方向的裂缝。参考图7,并结合文献[4]的计算情况和文献[8]所列的管片开裂类型及实际情况,可知提高管片制作和安装精度,对缓解施工过程中管片的应力集中,减少管片在施工状态时所需配筋有重要的意义。
(2)所有计算模型中,只有模型6管片外弧面出现裂缝。该模型最大压力与最小压力之比达到6:1,在单块管片的实际使用过程中难以产生如此悬殊的不均匀压力,故不妨认为管片外弧面在正常使用情况下不会出现裂缝。模型8的荷载为沿管片宽方向变化的“山”形荷载,在管片内弧面形成了大量沿管片长方向的连通裂缝。
(3)除(2)中提到的模型6和模型8外,其余计算模型,即模型4、模型5、模型7、模型9和模型10,其裂缝均只发生在纵向、环缝的手孔和螺栓孔附近,即不受钢筋保护的素混凝土保护层部分。只需加强薄弱部位的局部强度,就能很好抑制正常使用状态下裂缝在手孔和螺栓孔的出现与扩展。
(4)模型4、模型5、模型7、模型9和模型10中,极限荷载的平均值或峰值从(909.4～3 750.0) kN/m2。文献[6]反演计算管片法向压力值为(117.0～260.0) kN/m2。文献[9]对广州地铁某段盾构区间隧道进行实测,管片外侧土压力从(18.8～210.8) kN/m2。无论是反演计算结果还是实测数据,均远小于计算中所施加的极限荷载。因此,正常使用过程中,管片的素混凝土保护层部分一般不会出现如图7所示明显应力集中与裂缝开展。由此可以推断,管片配筋在正常使用过程中作用不明显,配筋的作用更多是用以承受千斤顶顶力带来的劈裂拉力[4]。
3结论
研究了施工状态和正常使用状态下单块管片的应力、裂缝分布、极限荷载与薄弱之处,从而明确,管片的制作和安装精度极大的影响了管片在施工状态时应力分布,环缝面的不平整更容易造成管片的应力集中而产生开裂,裂缝主要沿盾构推进方向发展,螺栓孔也存在开裂情况。正常使用状态下管片的薄弱之处为螺栓孔,裂缝主要在螺栓孔及其附近区域出现。由于有限元分析中所施加荷载的峰值远大于反演计算及实测的土压力结果,因而可以认为在正常使用状态下,管片内钢筋主要起架立作用,管片在施工阶段所需配筋量在一般情况下成为控制配筋量。只要能提高管片的制作和安装精度,用构造措施加强螺栓孔附近混凝土的抗裂性能,管片在使用过程中出现裂缝及破损的几率就会大大降低,也能节省盾构区间工程投资及隧道使用过程中的维修费用。
参考文献
[1] Working Group No.2 International Tunneling Association.Guidelinesfor the design of shield tunnel lining.Tunneling and UndergroundSpace Technology.2000;15(3):303—331
[2]李志南,陈丽娜.地铁盾构隧道管片配筋型式探讨.地下铁道新技术文集2003—中国土木工程学会隧道及地下工程学会地下铁道专业委员会第十五届学术交流会论文选集,成都,2003
[3] ADINA R&D Inc.Theory and modeling guide volume I:ADINASolids&structures.Watertown,2005
[4]李志南.关于地铁盾构隧道管片制作误差的讨论.中国土木工程学会快速轨道交通委员会学术交流会地下铁道专业委员会第十四届学术交流会论文集,北京,2001
[5]刘烈钢.盾构隧道管片抗弯吊装孔抗拔试验方法.广东水利水电.2005;6(增):39—40
[6]朱合华,崔茂玉,杨金松.盾构衬砌管片的设计模型与荷载分布的研究.岩土工程学报,2000;22(2):190—194
[7]广州市地下铁道设计研究院.广州市轨道交通三号线工程施工图设计(天河客运站～华师站区间),广州,2002
[8]竺维彬,鞠世健.盾构隧道管片开裂的原因及相应对策.现代隧道技术,2003;40(1):21—25
三维有限元法 篇10
关键词：圆筒式,水电站厂房,三维孔洞结构,有限元分析,大化水电站
0引言
大化水电站工程是红水河规划开发的第六个梯级水电站。扩建工程位于原枢纽左岸, 将原来的土坝部分拆除以布置厂房, 设计装设一台单机容量为110 MW的混流式水轮发电机组, 主要建筑物有河床式厂房、引水渠、尾水渠、刺墙坝、接头土坝等。
圆筒式主厂房左右侧接土石坝, 承受较大的土压力。主厂房长83.22 m, 高78.50 m, 上部33.98 m为完整的圆筒结构。圆筒内径33 m, 高程筒壁厚2.0～3.0 m。
大化扩建工程厂房结构上主要有3个特点:厂房基础存在t1、t2、t3等多条分布范围较大的破碎带, 且部分区域弱风化、强风化地质分界线较低;厂房与土坝直接连接, 厂房作为挡水结构, 承受较大的填土荷载和水荷载, 水土合力最高达1.05 MPa;主厂房采用圆筒式, 其受力特点与一般的矩形厂房有很大区别。圆筒式厂房在国内仅有广西乐滩水电站、洛东水电站扩建项目等工程使用过。为研究厂房结构在各种工况下的安全性和构件的应力及其配筋, 有必要对大化厂房结构进行三维有限元分析。
1计算条件、模型及工况
1.1计算软件
本项目采用的计算软件是大型通用计算软件ANSYS。ANSYS软件是融结构、热、流体、电磁、声学于一体的大型通用有限元分析软件。该软件具有强大的前处理及后处理功能, 它的图形界面和交互式操作大大简化了计算模型的创建过程, 同时在计算之前, 可通过图形显示来验证模型的几何形状、材料及边界条件;在后处理中, 其计算结果可以采用多种方式输出, 比如计算结果排序和检索、彩色云图、等值线、动画显示等等。
1.2计算假定、模型及边界条件
计算假定:副厂房、安装间分别置于主厂房的上部, 对主厂房而言为有利荷载, 且不承受水压力及填土等荷载, 可以不予考虑。同时, 不考虑1 m以下的孔洞。
计算模型:选取主厂房为分析对象, 基础分别向上游、下游和底部适当延伸。有限元模型共有实体单元69 733个, 结点107 230个。图1和图2分别给出了内外圆筒结构及其荷载示意图、蜗壳及其附近混凝土模型。
模型坐标系:X轴为水流向, Y轴为铅垂向, Z轴为水平横河向, 如图1所示。
计算边界条件:基岩的上、下游端面和左、右两侧按法向连杆模拟, 底部三向约束, 其他表面按力边界或自由边界考虑。
1.3材料参数
厂房的一期混凝土和二、三期混凝土分别采用C20、C25混凝土, 材料参数根据规范[1]采用;岩石主要力学指标如表1所示。另外填土的天然容重为19.3 kN/m3, 饱和容重为20.5 kN/m3, 凝聚力25 kPa, 内摩擦角为15°;两侧填土在主厂房上游端处高程为174.50 m, 在主厂房下游端高程为138.50 m。
1.4计算工况、荷载及其组合
共计算了5种工况, 计算工况、特征水位、荷载及其组合如表2所示, 其中渡汛工况厂房蜗壳二、三期混凝土未浇筑。
2计算成果与分析
通过有限元分析, 得到该水电站厂房坝段在各种工况下的位移和应力分布规律。由于厂房结构复杂, 取厂房关键部位控制工况结果进行分析如下。
2.1厂房总位移分析
图3为校核工况的厂房总位移图, 总位移大小为18.24～50.84 mm。就厂房整体结构而言, 厂房各个方向位移分布较连续均匀。就构件而言, 上游墩墙特别是下部承受较大的水压力及填土荷载, 变形较明显;圆筒上下游承受非对称荷载, 变形较明显;下游右侧墙承受水土荷载, 且墙体较薄, 变形较明显;基础软弱且不均匀, 位移较不均匀。
2.2圆筒应力分析
图4为渡汛工况圆筒环向拉应力等值线图。圆筒的环向拉应力主要分布于圆筒上部的左右两侧, 圆筒的环向主拉应力在0.8 MPa以内, 范围较大。圆筒变形和应力主要由于圆筒直接或间接承受上下游和左右侧不对称荷载所致。从计算结果来看, 圆筒结构的优点是其受力条件好, 可以大大减少工程量。缺点是不能分缝, 布置比较困难, 施工比较复杂, 且直径太大时其受力条件的优势不能显现。
2.3底板应力分析
厂房底板出现一定的应力集中, 但区域较小, 主要分布于底板t1挤压破碎带、t3挤压破碎带附近, 主要为顺水流向应力。相对其他工况而言, 校核工况底板受地基软弱和扬压力双重影响最大, 其拉应力分布较广, 最大拉应力值达3.32 MPa, 除在拉应力分布区域配筋外, 建议在t1挤压破碎带附近补强。
2.4墩墙应力分析
墩墙第一拉应力主要分布于上游部分的内外侧, 最大拉应力为1.33 MPa。主要原因在于墩墙外侧水土荷载很大, 竣工工况墩墙内外荷载差比较大, 且二、三期混凝土未浇筑, 墩墙的拉应力分布较广, 且应力值较大。
2.5厂房蜗壳分析
图7、图8为蜗壳周边混凝土环向应力等值线图。竣工工况蜗壳及其附近混凝土内外承受较大的荷载差, 环向拉应力主要分布于内环的左右两侧和外环的下游侧, 且分布较广较深, 最大环向拉应力1.32 MPa。渡汛工况蜗壳二、三期混凝土未浇筑, 厂房下部蜗壳附近的筒体较单薄, 承受较大的水土荷载, 最大环向拉应力为2.35 MPa, 分布较深;蜗壳浇筑后, 其应力水平明显降低。
2.6配筋分析
对出现拉应力的区域进行配筋分析和配筋, 并对出现较大拉应力的区域如底板、蜗壳等进行了抗裂分析。结果表明, 配筋满足抗裂要求。鉴于篇幅要求在此不详列。
3结语
(1) 用三维有限元方法对厂房整体结构进行分析, 充分反映了厂房空间特性对结构应力和变形的影响, 对工程的设计和优化有一定的指导意义。
(2) 圆筒结构承受上下游和左右侧非对称荷载, 主要在顶部左右侧出现拉应力, 且拉应力不大。说明圆筒结构受力条件好, 可以大大减少工程量;但是不能分缝, 布置比较困难, 施工比较复杂, 且直径不宜太大。
(3) 计算分析表明:在厂房圆筒、上下游墩墙、底板、圆筒与墩墙交界处等出现了一定的拉应力, 数值一般不大, 分布范围有限;渡汛工况蜗壳附近应力状态不利, 建议从结构、施工等方面予以注意;整个结构的压应力均在混凝土的抗压强度以内。
(4) 厂房结构各部位最大拉应力并不由一个工况控制, 因此应综合考虑各工况的应力成果进行配筋, 尤其是对拉应力较大的部位配以足够的钢筋, 必要时辅以结构和施工措施, 以确保结构安全。
参考文献
[1]SL/T 191-96, 水工混凝土结构设计规范[S].
[2]SL 266-2001, 水电站厂房设计规范[S].
[3]顾鹏飞.水电站厂房设计[M].北京:水利电力出版社, 1985.
【三维有限元法】相关文章：
三维模型有限元05-25
三维时域有限差分法06-19
有限元法论文06-18
塑性有限元法05-10
有限元计算法06-20
随机有限元法06-21
有限元法课程的教与学06-05
西安盾构法施工引起的地面沉降有限元分析09-11
基于三维学习法的中学语文论文05-01