应力模拟

应力模拟（精选九篇）

应力模拟 篇1

关键词：曲线论,摩擦损失,沿程荷载分布力

在工程实践中,处理具有空间曲线索形预应力筋的预应力损失时,按二维函数考虑,忽略空间效应;或者将曲线分别投影到平面和立面上成二维函数并计算曲率转角,然后再累加各面上计算得到的预应力转角摩擦损失。这在一定条件下满足工程的精度要求,但对于布置较复杂的预应力不能准确地反映其对结构的实际作用,此外计算结果也存在较大的误差。

目前,有部分商业软件已经采用空间曲线函数来计算空间预应力筋,如ANSYS土木工程模块预应力筋线型编辑器中应用的贝塞尔曲线(Bezier curve)和文献[1]中的双圆曲线(bicircular curve)。文中根据曲线论和预应力损失相关公式推导出空间预应力筋的摩擦损失,然后根据预应力对结构作用的特点,求解出预应力沿程荷载分布力,为编写三维预应力程序提供理论依据。

1 预应力筋摩擦损失计算

基本公式推导。设某一空间预应力筋的曲线方程为:

P(t)=(x(t),y(t),z(t)) (1)

取预应力筋一微段作为分析对象,由沿其切线方向的力的平衡条件可知:在预应力筋微段两端作用的集中力分别为N(t)和N(t+dt),方向是沿该端点切线方向(见图1),大小等于预应力筋的截面积乘以该点初应力。

N(t+dt)=N(t)-dN(t) (2)

由曲线论知[4],空间曲线在某点的曲率为$\phi (t)\text{d}l=|{\overrightarrow{{\rm P}}}^{\prime }\times {\overrightarrow{{\rm P}}}^{″}|/|{\overrightarrow{{\rm P}}}^{\prime }|{}^3$,该微段长度$\text{d}l(t)=|{\overrightarrow{{\rm P}}}^{\prime }|\text{d}t$,弯曲角$\text{d}\theta (t)=\phi (t)\text{d}l=(|{\overrightarrow{{\rm P}}}^{\prime }\times {\overrightarrow{{\rm P}}}^{″}|/|{\overrightarrow{{\rm P}}}^{\prime }|{}^2)\text{d}t$,则摩擦损失值为:

dN(t)=N(t)(μdθ(t)+kdl(t)) (3)

把dl(t),dθ(t)公式代入式(3),得:

$\text{d}{\rm N}(t)={\rm N}(t)(\mu |{\overrightarrow{{\rm P}}}^{\prime }\times {\overrightarrow{{\rm P}}}^{″}|/|{\overrightarrow{{\rm P}}}^{\prime }|{}^2+k|{\overrightarrow{{\rm P}}}^{\prime }|)\text{d}t$ (4)

其中,μ为预应力与管道的摩擦系数;k为管道每米局部偏差对摩擦的影响系数。实际计算摩擦损失[5]σ1时只需要计算从张拉端至计算截面曲线孔道部分夹角之和$\theta ={\displaystyle \int |}{\overrightarrow{{\rm P}}}^{\prime }\times {\overrightarrow{{\rm P}}}^{″}|/|{\overrightarrow{{\rm P}}}^{\prime }|{}^2\text{d}t$,从张拉端至计算截面曲线孔道长度$L={\displaystyle \int |}{\overrightarrow{{\rm P}}}^{\prime }|\text{d}t$,则可得:

$\sigma {}_1=\sigma {}_{\text{c}\text{o}\text{n}}[1-e{}^{-(\mu {\displaystyle \int |}{\overrightarrow{{\rm P}}}^{\prime }\times {\overrightarrow{{\rm P}}}^{″}|/|{\overrightarrow{{\rm P}}}^{\prime }|{}^2\text{d}t+k{\displaystyle \int |}{\overrightarrow{{\rm P}}}^{\prime }|\text{d}t)}]$ (5)

其中,σcon为锚下控制应力。

2 预应力沿程荷载分布力的确定

假想将预应力筋和结构“脱离”,用作用力和反作用力来模拟预应力筋与结构间的相互作用[6]。将预应力筋看作一个脱离体,作用在该脱离体上的力系必定是一个静力平衡力系。按力的平衡可求得作用于脱离体上的力,再将该力反加到结构上,这样就模拟了预应力筋对结构的作用。

基本公式推导。沿预应力筋延伸方向分别取微段,考虑预应力筋内部存在初应力,即预拉力,由力的平衡条件可分析出预应力筋沿程力的分布情况。取预应力筋弯曲部分的一微段作为分析对象(见图1),由力的平衡条件可求得该处结构对预应力筋的作用力$\overrightarrow{F(t)}$为法向力和切向力的矢量和,具体计算公式如下:

$\overrightarrow{F(t)}=\overrightarrow{{\rm N}(t+\text{d}t)}+\overrightarrow{{\rm N}(t)}$ (6)

考虑$\overrightarrow{v(t)}$为曲线在t点的单位切向量,方向为曲线的行进方向。

$\overrightarrow{F(t)}={\rm N}(t+\text{d}t)\overrightarrow{v(t+\text{d}t)}-{\rm N}(t)\overrightarrow{v(t)}$ (7)

$\overrightarrow{v(t+\text{d}t)}$泰勒展开式为:

$\overrightarrow{v(t+\text{d}t)}=\overrightarrow{v(t)}+\overrightarrow{v^{\prime }(t)}\text{d}t+\overrightarrow{v^{″}(t)}\text{d}t{}^2/2$ (8)

将式(4),(5),(6),(8)代入式(7)并忽略高阶小量得:

$\overrightarrow{F(t)}={\rm N}(t)\overrightarrow{v^{\prime }(t)}\text{d}t-{\rm N}(t)(\mu |{\overrightarrow{{\rm P}}}^{\prime }\times {\overrightarrow{{\rm P}}}^{″}|/|{\overrightarrow{{\rm P}}}^{\prime }|{}^2+k|{\overrightarrow{{\rm P}}}^{\prime }|)\overrightarrow{v(t)}\text{d}t$ (9)

其中有:

$\overrightarrow{v(t)}={\overrightarrow{{\rm P}}}^{\prime }/|{\overrightarrow{{\rm P}}}^{\prime }|$ (10)

$\overrightarrow{v^{\prime }(t)}=\frac{{\overrightarrow{{\rm P}}}^{″}|{\overrightarrow{{\rm P}}}^{\prime }|-{\overrightarrow{{\rm P}}}^{\prime }|{\overrightarrow{{\rm P}}}^{″}|}{|{\overrightarrow{{\rm P}}}^{\prime }|{}^2}=\frac{{\overrightarrow{{\rm P}}}^{″}|{\overrightarrow{{\rm P}}}^{\prime }|{}^2-(|{\overrightarrow{{\rm P}}}^{\prime }\cdot {\overrightarrow{{\rm P}}}^{″}|){\overrightarrow{{\rm P}}}^{\prime }}{|{\overrightarrow{{\rm P}}}^{\prime }|{}^3}$ (11)

最终将式(10),(11)代入式(9)得到沿程荷载分布力集度计算公式为:

$\overrightarrow{Q(t)}={\rm N}(t)\frac{{\overrightarrow{{\rm P}}}^{″}|{\overrightarrow{{\rm P}}}^{\prime }|{}^2-({\overrightarrow{{\rm P}}}^{\prime }\cdot {\overrightarrow{{\rm P}}}^{″}){\overrightarrow{{\rm P}}}^{\prime }-\mu |{\overrightarrow{{\rm P}}}^{\prime }\times {\overrightarrow{{\rm P}}}^{″}|{\overrightarrow{{\rm P}}}^{\prime }-k|{\overrightarrow{{\rm P}}}^{\prime }|{}^3{\overrightarrow{{\rm P}}}^{\prime }}{|{\overrightarrow{{\rm P}}}^{\prime }|{}^3}$ (12)

编写有限元程序时常需要求出某段预应力筋t∈[a,b]对结构的作用力,即沿程荷载为:

3算例

对于平面问题进行分析。不考虑摩擦损失,控制张拉力为F。这段曲线的函数表达式如下:

通过式(12)可以计算得沿程荷载分布力集度为:

Q(x)即为文献[7]中内荷载集度的表达式。因此,可以说文献[7]所推导的直接内载法为文中所推导的沿程荷载分布力集度公式的特例。

4结语

预应力对结构的作用是由预应力筋的径向力、预应力筋端部的锚固力和预应力束对结构的摩擦力组成的。

文中所推导的沿程荷载分布力公式能全面地反映预应力沿程荷载分布情况,满足自平衡条件并且计算简单易于理解,可以为空间有限元程序的编写以及预应力结构理论分析和设计提供理论依据。

此外,通过算例的比较,说明了考虑第一类损失得到的沿程荷载分布力计算公式的通用性和工程实用性。本公式适用于弯桥等预应力布置复杂的结构。在准确提供预应力筋曲线方程的基础上,能很好地反映预应力筋对结构的作用力。

参考文献

[1]Park S K,Kim K S.A new design method for prestressed concretecontinuous flat slabs[J].Struct Method for Tall Spec Build,2004,13(4):265-276.

[2]Donald Hearn,M Pauline Baker.Computer Graphics,C version[M].2nd ed.Beijing:Tsinghua University Press,1998.334-345.

[3]施法中.计算机辅助几何设计与非均匀有理B样条[M].北京:高等教育出版社,2001.273-280.

[4]苏步青.实用微分几何引论[M].北京:科学出版社,1986.19-28.

[5]李国平.预应力混凝土结构设计原理[M].北京:人民交通出版社,2000.31-33.

[6]叶见曙,袁国干.结构设计原理[M].北京:人民交通出版社,1999.223-242.

应力模拟 篇2

等应变率加载下冻土粘弹塑性应力应变本构模拟

借助模型分析,通过等应变率下对不同含水量、不同温度的兰州粉结冻土试样在单轴加载条件进行应力测试,得到应力应变曲线.根据粘弹塑性体模型,结合试验曲线特征,选择应力应变的方程形式.以其中一个组典型的.试验结果为例进行拟合数据处理,获得拟合的应力应变曲线,确定方程中的参数.再将方法用于不同情况下的实验曲线处理,得到10%和20%两种含水量:-10℃和-20℃两种试验温度下冻土试样的测试的应力.应变三阶段全过程本构方程.结果表明,由此得到的方程能很好地描述冻土的材料特性.

作 者：李健康 LI Jian-kang  作者单位：沙洲职业工学院,江苏张家港,215600 刊 名：沙洲职业工学院学报 英文刊名：JOURNAL OF SHAZHOU POLYTECHINICAL INSTITUTE OF TECHNOLOGY 年，卷(期)： 12(3) 分类号：P642.14 关键词：人工冻结土   粘性材料   应力应变曲线   本构方程  

应力模拟 篇3

关键词：激光；热应力；数值模拟；ANSYS

中图分类号：TG335文献标识码：A文章编号：1007-9599 (2010) 13-0000-02

The Numerical Simulation

of Metal Plate Thermal Stress under the Action of Laser

Mao Jingfeng,Jiang Houman

(National Defense University,Institute of Optical Science and Engineering,Changsha410073,China)

Abstract:This paper finite element software for ANASYS make brief introduction,on this basis,focusing on the use of the software on the metal plate under the laser temperature field and thermal stress field of process simulation,and gives the simulation results.Finally,the use of finite element software on the metal plate to simulate the thermal stress the significance and purpose are explained.

Keywords:Laser;Thermal stress;Numerical simulation;ANSYS

激光作用于金属平板时，平板表面由于吸收了大量的激光能量，将会引起温度升高、熔融、气化、喷溅等现象，激光与金属平板作用的具体过程与激光参数（能量、波长及脉宽等）、材料特征和环境条件等密切相关。

求解激光作用下结构体热力耦合问题，主要有解析方法和数值方法两大类。解析方法首先建立数学模型，然后从数学模型出发，利用各种数学工具，采用各种技巧，得到某种条件下数学问题的解析形式的解答。随着工程技术的发展，出现了大量复杂的数学物理问题，比如大型结构的强度分析和震动分析等，求解这些复杂的问题，传统的解析方法显得无能为力，需要借助于数值方法，有限元法正是在这种情形下应运而生的。

一、有限元简介

有限元方法是目前在结构分析中应用最广泛的一种方法，另外它还被广泛应用于各种场问题如热传导、电磁场、流体流动等等，目前，有限元方法已成为研究热力耦合问题最重要的数值方法。ANSYS是一个通用的有限元分析软件，可以模拟热与结构应力之间关系的耦合问题的分析。采用ANSYS对激光加热时金属平板的热应力进行数值模拟，能综合考虑实际条件，较好地模拟激光与物质的相互作用。

热应力是由于相互接触的不同结构体或同一结构体的不同部分之间的热膨胀系数不匹配，在加热或冷却时彼此的膨胀或收缩程度不一致，从而导致热应力的产生[2]。

热应力问题实际上是热和应力两个物理场之间的相互作用，属于耦合场分析问题。运用ANSYS有限元软件对结构体进行热力分析可采用直接法和间接法两种不同途径。直接法是直接采用具有温度和位移自由度的耦合单元，同时得到热分析和结构应力分析的结果；间接法是指先进行热分析，然后将求得的节点温度作为体载荷施加到结构应力分析中。

二、数值模拟

（一）模拟条件

设定金属平板几何尺寸（长×宽×厚）为：120mm×120mm×10mm。平板材料采用常用金属 钢作为研究对象，该材料密度为 ，导热率为30W/m℃，比热容为544J/kg℃，线膨胀系数为13.6×10-6/℃，弹性模量为191Gpa，剪切模量为73.3Gpa，泊松比为0.3，熔点为1524℃，设定初始温度为25℃。假设激光作用在金属平板的中心，光斑直径为30mm，功率密度为1500W/cm2，加热时间为5s。

（二）模拟过程

本文采用间接法对激光辐照下金属平板的热应力进行数值模拟，其基本步骤如下：

1.实体建模，首先根据分析类型和特点选择无中间节点的3-D热分析单元Brick 8node 70作为分析的单元类型，然后根据前面给出的平板材料参数对材料的性能参数进行定义，并根据平板的几何尺寸通过体素在图形窗口中建立平板的实体模型。为便于激光辐照的加载，需要在平板模型的中心位置创建一个与激光光斑直径相同的圆柱体，然后将圆柱体与平板模型进行布尔运算。创建好的实体模型如图1所示。

图1平板的实体模型 图2施加热流密度后的有限元模型

2.网格划分，手动指定实体模型上各条线段被划分的数量，然后程序自动按要求对实体模型进行网格划分。为了提高模拟精度，同时节省计算机资源，对实体模型的不同区域指定了不同的网格密度，模型的中心附近区域在激光辐照过程中温度变化梯度较大，网格划分的更稠密。完成网格划分后的实体模型即已转化为有限元模型。

3.加载求解，首先对平板的初始温度进行设定，然后对激光能量进行加载。将激光能量以热流密度的方式施加在有限元模型的中心部位，热流密度的数值根据激光能量的大小和平板材料对激光的吸收率确定。施加热流密度后的有限元模型如图2所示。

4.溫度场求解，首先对求解过程进行控制设置，主要包括求解类型设置、求解选项设置以及求解时间子步的设置等，完成设置后计算机将自动完成求解。完成求解后可通过通用后处理器对模型的温度场分布进行查看，如图3所示，此外还可通过时间历程后处理器对模型上任意指定点的温度-时间曲线进行显示，图4中显示了平板中心点温度随时间的变化关系曲线。

图3温度场分布等值线图 图4中心点温度随时间变化关系曲线

5.热应力求解，完成温度场的求解后重新进入前处理器，进行单元类型的转化，将热单元转换为相应的结构单元；然后对结构分析中的弹性模量、泊松比、线膨胀系数等材料属性进行设置；接着读入热分析结果将其作为载荷，并指定计算热应力时的参考温度为25℃；最后对热应力进行求解。平板的等效应力等值分布图如图5所示。

（三）模拟结果

根据2.1节设定的条件进行模拟，由模拟结果可知金属平板受激光作用后，平板中心点的温度最高，温度可达到341.6℃；由于激光作用时间较短，金属平板吸收的激光能量还来不及传导，激光光斑附近区域温度梯度变化较大，远离激光光斑的区域温度基本没有升高。金属平板内部最大应力值为789MPa，小于 钢的抗拉强度。

三、结论

运用该软件能够根据金属平板的几何尺寸、材料性能以及激光的功率密度、光斑大小、激光作用时间，材料对激光的吸收率等相关参数对受激光作用下金属平板内部的温度场和应力场分布进行有效数值模拟。通过模拟结果可以预先估计平板内部最高温度点及最大受力点的大小和位置，判断金属平板内部热应力是否超出材料应力极限。

此外，ANASYS有限元软件作为一种有效工具，也能够根据实际情况便捷、准确地对各种结构体受激光作用时的温度场、应力场进行数值模拟，可广泛应用于激光打孔、激光切割、结构强度分析等方面。

参考文献：

[1]陆健,倪晓武,贺安之.激光与材料相互作用物理学.北京:机械工业出版社,1996

[2]张朝晖.ANASYS热分析教程与实例解析.北京:中国铁道出版社,2007

地应力影响因素三维数值模拟分析 篇4

地壳岩体中的应力状态与人类工程活动的关系极大,它不仅是决定区域稳定性和岩体稳定性的重要因素,而且往往对各类建筑物设计和施工造成直接的影响。因此,地应力场研究也就成为大型工程建设中的一项重要课题。

在大型水电站工程中,针对坝址区岩体初始地应力场的反演已做了大量研究,方法已经成熟,而针对坝区河谷初始地应力场反演影响因素分析的并不多,本文以某水电工程为实例,结合三维数值模拟,分析河谷在下切过程中受断层、错动带以及侧向剥蚀等因素的影响地应力变化研究。

1 工程实例

某水电站位于金沙江上,河谷呈不对称的“V”字形,两岸为单斜山,左岸谷肩以上为斜坡地形,斜坡沿岩流层层面发育,倾向SE,倾角15°左右,谷肩以下为临江陡壁地形,陡壁间有2级～3级缓坡台地相连;右岸谷肩以上为缓坡地形,岩层倾向SE,倾角10°～15°,近谷肩附近覆盖层厚度5 m～10 m。谷肩以下为陡壁地形,岸坡高陡,间有狭窄的缓坡台阶。此段河流总体流向为由南向北流。

通过对实测地应力资料的分析整理,研究区地应力场实测应力场以NE向为主。在临空面附近的斜坡岩体中,实测最大主应力量值在5 MPa以内,方位角与河流向近于平行,倾角基本上与坡面近于平行,最小主应力方位角与河流向则近于垂直,倾角与坡面也近于垂直;最大、最小主应力从边坡表面向坡内随深度增大而增大,在断层穿越带应力产生分异,整体上以压应力为主;由于研究区范围较大,为了分析更准确,计算结果更接近实际值[1]。

2 三维有限元模型

现今河谷应力场主要是在区域应力场背景下,由于河谷下切、岩体卸荷导致一定范围内谷底、谷坡应力场调整所致。在建模过程中对模型概化作如下考虑[2]:

1)岩流层。

由于模型选取的范围大,为了便于模型的建立与计算,模型中没有具体区别不同类型的玄武岩和角砾熔岩,而将其作为同一种材料。

2)结构面。

为了突出重点,模型中仅考虑了错动较为强烈的F17,F14,F16,F19,F20,C2,C3。

3)地形。

考虑到研究区河谷应力场主要受河谷形状的影响,故建模时只考虑河谷的整体形态及右岸的深沟,而忽略了两岸边坡的微地形。

4)阶地。

研究区金沙江河谷在坝址区段经过了宽谷期和峡谷期两个过程,一共经过了五次快速下切,才形成目前河谷形态。在进行数值模拟计算时,宽谷期考虑为一次形成,而峡谷期则为四次下切形成。

5)计算模型。

模型X轴正方向对应地理坐标的正东向(垂直金沙江),Y轴正方向对应地理坐标的正北向(金沙江流向),Z轴正方向对应地理坐标竖直向上。模型中全部采用八节点四面体单元,总单元数为502 765个,节点数为695 726个(见图1,图2)。

前人研究表明,该水电站区域构造应力场的量级为8 MPa～12 MPa,方向为NE向。而我们对研究区地应力实测资料的统计分析的结果表明,厂房区地应力实测结果,深部最大主应力方向为N30°～80°E,量级在19 MPa～22 MPa之间;最小主应力方向为N64°～73°W,量级为6 MPa～9 MPa。因此,构造应力应沿N30°～80°E方向施加。在FLAC3D计算过程中,将该构造力分解为X方向的力、Y方向的力以及τxy,τxz,τyz施加到各个单元的节点上。

3 影响因素应力规律分析

计算后,从所筛选的两岸岸坡及地下厂房应力进行拟合,量值上,此次拟合计算70%模拟应力值相对误差小于20%,其中30%的模拟值相对误差达到40%左右。方向上,方位角走向基本一致。中间主应力及最小主应力倾角均拟合较好,低高程拟合效果较高。因此,计算结果具有一定的真实性,可以为下面影响因素分析提供数据依据。

3.1 断层及错动带的影响分析

模型中考虑了断层F17,F14,F16,F19,F20和层间错动带C2,C3。河谷应力场反演分析提取出来的数据表明,断层及错动带附近地应力都发生突变。

根据模拟结果提取出来的数据可知,研究区断层及错动带对地应力的影响可归纳为如下几点:

1)C2对地应力有放大作用,影响范围约为100 m。如图3～图5所示,在C2通过的高程,影响范围内岩体的应力增大,增大量值约1 MPa～3 MPa,影响范围外的左右两侧应力大小相等;

2)C3与C2不同,对附近岩体应力的影响则是先放大然后缩小。由图6～图8可以看出,C3上盘应力值明显大于下盘应力值,且上下盘应力量值差约为4 MPa。在C3的上盘应力值明显增大,但当过渡到下盘时应力值又突然减小;

3)F17,F19,F16附近岩体应力降低,如图4,图7所示。F14,F20对岩体应力影响不大,主要表现在F14,F20通过的地方岩体应力略有下降但不是很明显,如图3所示。

3.2 侧向剥蚀的影响分析

侧向剥蚀主要影响左岸应力分布,而对右岸应力分布的影响不是很明显。在河谷应力反演过程中,根据追踪线上的应力可以看出,第一阶段、第二阶段的宽谷河床岩体应力变化规律基本上相同,在这两个阶段研究区岩体主要受垂向剥蚀的影响,从第三阶段峡谷初期开始到第五阶段风化卸荷期,这三个阶段的应力变化规律基本上相同,研究区岩体在这三个阶段主要受侧向剥蚀的影响,岩体应力的这两次较大变化,主要体现在受层间错动带C3影响的地方,如图9～图12均受C3影响,在第一、二阶段受C3影响范围内的岩体在C3上盘岩体应力值明显减小,C3下盘岩体应力值明显增大,第三、四、五阶段C3附近岩体应力变化刚好与前两个阶段相反,即在C3的上盘则表现为岩体应力值明显增大,而C3下盘则表现为岩体应力值明显减小。

4 结语

1)在断层、断层与断层及断层与层间错动带交错部位岩体应力明显减小,在层间错动带C2上岩体应力值明显增大,在层间错动带C3上盘岩体应力值增大,C3下盘岩体应力值减小,且在断层及层间错动带发育的部位最大主应力的方向表现为平行结构面。侧向剥蚀主要影响左岸应力分布,而对右岸应力分布的影响不是很明显。

2)影响地应力的因素有很多,除上述两种外,风化剥蚀、地形地貌对地应力的影响也很大。

参考文献

[1]李攀峰.金沙江溪洛渡水电站坝区地应力场及地下洞室群围岩稳定性数值模拟[D].成都:成都理工学院硕士学位论文,2001.

[2]柴贺军,刘浩吾,王明华.大型电站坝区应力场三维弹塑性有限元模拟与拟合[J].岩石力学与工程学报,2002,21(9):1314-1318.

[3]张建国,张强勇.大岗山水电站坝区初始地应力场反演分析[J].岩土力学,2009,30(10):3071-3078.

冷弯薄壁C形钢的残余应力模拟 篇5

残余应力是使轴心受力构件提前进入非弹性阶段的主要因素,对稳定性影响比较大,是一个不容忽视的问题。因此,搞清楚残余应力对冷弯薄壁型钢轴心受压构件整体稳定的影响,是非常必要的[1,2]。

由于测量方法的限制,冷弯残余应力的研究曾一度停滞,直到1988年Weng和Pekoz[1]两位教授在康奈尔大学应用电化学腐蚀切条法成功测出冷加工板材的纵向残余应变,才使得冷弯残余应力的研究重新回到人们的视野。通过研究,人们发现冷弯构件的残余应力与热轧型钢的残余应力有着本质区别。但由于残余应力等初始缺陷对构件极限承载力影响的复杂性,我国《冷弯薄壁型钢结构技术规范》(GB50018—2002)[3]仍未直接考虑残余应力对构件整体稳定极限承载力的影响。杆件进入非弹性阶段以后,主要运用切线模量代替弹性模量处理(Et=τE),折减系数τ的取值是根据构件试验结果分析确定的,并没有从残余应力角度去分析和处理。本文为解决这个问题,采用了编写初始应力文件法给有限元施加残余应力。即在外部软件中根据残余应力的分布规律编写初始残余应力文件,通过ANSYS将该文件读入到有限元构件进行计算。模拟云图与计算假定吻合良好,从而证明了运用初始应力文件法模拟残余应力是可行的。

1 实测冷弯残余应力

1988年, Weng等[1]在康乃尔大学运用电化学腐蚀切条法,量测得到冷加工卷边薄壁槽形钢内外残余纵向应变的实测值,如图1所示。图1中εzrm的最大值在弯曲区内外壁处,外壁拉应变比0.5εs大;非变形区的εzrm值最小约为0.25εs;中间弯曲区的值较弯曲卷边区大,这里的εs是材料屈服时的应变值。1989年,赵峰[4]运用电化学腐蚀切条法得到卷边角钢的内外壁上冷弯残余应变的实测结果,如图2所示。弯曲区的εzrm值最大,中心弯曲区的内外壁为0.4εs,卷边弯曲区的内外壁为0.34εs。非变形区的εzrm值很小,为0.05εs。

从残余应力实测图形可以得出冷弯构件残余应力的一些基本规律[1,4]:构件外壁纵向为残余拉应力,构件内壁纵向为残余压应力;截面同一位置处,内外表面残余应力绝对值大体相等;构件转角处的纵向残余应力大于平板处纵向残余应力,纵向残余应力峰值出现在构件转角处;沿横截面周边其纵向拉、压残余应力分布是不均匀的。

2 有限元模拟分析

2.1 冷弯C形钢残余应力的计算模型

根据冷弯残余应力分布规律,本文采用以下冷弯残余应力计算模型(图3、图4):①截面外侧为残余拉应力,内侧为残余压应力;②构件转角处残余应力应力取峰值0.4fy,同一截面处内外壁残余应力拉压绝对值相等;③平板段的纵向残余应力均匀分布为0.2fy;④残余应力沿板件厚度方向呈线形变化,中性层处应力为0。

2.2 残余应力文件格式

在进行结果分析时,可以把残余应力看成一种初始应力(预应力)。ANSYS允许使用输入文件来把初始预应力指定为一种荷载,预应力荷载可用于静态分析或完全的瞬态分析,这种荷载只能在分析的第一荷载步中施加,分析可以是线性或非线性。ANSYS中只有以下单元支持初应力文件输入:PLANE2、PLANE42、SOLID45、PLANE82、SOLID92、SOLID95、LINK180、SHELL181、PLANE182、PLANE183、SOLID185、SOLID186、SOLID187、BEAM188和BEAM189。在初应力文件编制中,初应力可以指定在单元的中心或者单元的积分点处。初应力文件是一个ASCⅡ文件,可以用外部编辑器进行编辑,文件的扩展名为.ist。文件中的注释性语句的第一个字符必须是“!”号。每个单元记录的第一行使用字符串“EIS”标明,随后是单元编号和一个可选的当地位置符。这些项目必须以“,”隔开。每个单元记录的第一行后有好几行,将指定单元中每个应力点的单元应力记录。本文对应有五个积分点,所以有五行。ANSYS要求在每个应力记录中有6个应力张量分量(Sx、Sy、Sz、Sxy、Syz、Sxz)。如果一个单元坐标系(ESYS)是为一个单元定义的,初始应力的指定必须在这个坐标系中。具体格式见文献[6]。由于冷弯残余应力大小沿杆件长度方向不变,本文应用C++语言编写了残余应力转化程序,解决了因单元量巨大所导致的初应力文件编写困难的问题。在编写残余应力文件时,只需按初应力文件格式写出构件底(或顶)端截面每个单元应力记录,经程序转化就可直接形成构件所有单元应力记录,可以节约大量的时间。

2.3 有限元模型

编辑完残余应力文件后,用ISFILE命令读取初应力文件。图5为采用壳SHELL181单元模拟的冷弯薄壁C形钢残余应力云图,尺寸为h×b×a×l=160 mm×60 mm×20 mm×2.5 mm,材质为Q235钢,fy=235 MPa,Es=206 000$\frac{{\rm N}}{\text{m}\text{m}{}^2}$,泊松比为0.3。将每个单元的残余应力分别指定在单元积分点处。对于冷弯残余应力沿厚度方向变化,可采用非标准积分点位置加以模拟。沿单元厚度方向设5个积分点,厚度坐标从-0.5到0.5,积分点位置依次为0.5,0.3,0.0,-0.3,-0.5,见图6。冷弯残余应力沿厚度方向的变化,可通过把残余应力指定在这5个积分点上来描述。ANSYS计算所得的残余应力模拟云图如图5所示。

从残余应力云图得知,最大残余应力出现在构件转角部位且沿构件纵向分布;除去两端约束对构件残余应力略有影响外,中间平板处残余应力沿纵向大体呈均匀分布。模拟云图和计算模型吻合较好,能够比较真实地反映构件中的残余应力状态。对于各种不同的截面构件、不同残余应力分布方式,都可采用该方法模拟其残余应力。

2.4 屈曲分析中的应用

运用ANSYS软件对冷弯薄壁C形钢轴心受压构件进行屈曲分析步骤如下:运用SHELL181单元建立模型,划分网格,施加支座约束;稳定分析时,首先用ISFILE命令读入残余应力文件,进行特征值屈曲分析,求出极限荷载上限;接着进行非线性屈曲分析,根据已求得的.ist文件给有限元构件施加$\frac{l}{1000}$初始几何缺陷。使用弧长法求得构件的荷载—位移曲线,根据曲线判断其极限荷载值,计算出轴心受压整体稳定系数φ。运用该方法可以求得直接考虑冷弯残余应力影响的轴心受压构件的弹塑性极限荷载。

3 结语

本文通过外部软件编写初始应力文件、再通过ANSYS将其读入有限元构件,可以很好模拟冷弯薄壁C形钢的残余应力。从而可减少运用直接法测量残余应力[7]所耗费的人力财力。运用计算机模拟残余应力,可以为冷弯薄壁C形钢的其它研究提供便捷。

参考文献

[1] Weng C C,Pekoz T.Residual stress in cold—formed steel members.Journal of Structural Engineering,1990;116(4):1161—1165

[2] Quach W M,Teng J G,Chung K F.Finite element predictions of re-sidual stresses in press-braked thin-walled steel sections.EngineeringStructures,2006;(28):1609—1619

[3]中华人民共和国国家标准.钢结构设计规范(GB50017—2003),北京:中国计划出版社,2003:10

[4]赵峰.考虑残余应力影响冷弯薄壁型钢轴心压杆整体稳定的计算方法.工业建筑,1994;24(5):50—55

[5]刘华琛,赵峰.模拟冷弯残余应力及其在稳定分析中的应用.国外建材科技,2005;26(2):52—53

[6]武思宇,罗伟.ANSYS工程计算应用教程.北京:中国铁道出版社,2004

重轨水平模拟矫直后残余应力的分析 篇6

在我国铁路主要干线上, 60 kg/m重轨的铺设已明显上升, 在铁路运输中, 它的优势越来越明显。这样, 对重轨质量就提出了更高、更严的要求, 国家冶金局和铁道部专门出台了《铁路客运专线建设计划》和《时速200 km客运专线60 kg/m重轨暂行技术条件》, 首次对重轨的残余应力提出了要求。

重轨的残余应力是由重轨矫直机的结构参数和矫直工艺参数决定的, 因此对重轨矫直进行全面、综合的深入研究对提高重轨的质量、提高成材率 (国产重轨成材率只有85%) 、延长使用寿命等具有重要意义。

2 辊式矫直机的矫直原理

重轨的矫直原理是:重轨经过辊式矫直机时, 通过交错排列的矫直辊, 多次反复弯曲而多次发生反复的弹塑性变形, 从而使重轨的原始弯曲逐渐减小, 最终获得合格的重轨。

在矫直过程中, 重轨表层的变形如公式 (1) 所示:

undefined

在矫直过程中, 压弯量增大时残留量的差值减小, 当压弯次数增加时残留量的差值也减小, 递减量合适时残留量才能趋近于零。重轨在各矫直辊下的残余曲率如图1所示。

3 显示动力学有限元方程

重轨在矫直过程中存在材料非线性、几何非线性以及边界接触非线性, 既要计算弹性变形又要计算塑性变形, 目前在数值模拟中应用得最广泛的是有限变形弹塑性法。重轨矫直过程中的变形是大位移弹塑性变形, 由于变形的非线性, 一般要采取增量的方法进行分析, 把力和位移或应力与应变之间的非线性关系在每一个增量中处理为线性关系, 来求整个过程的近似解。重轨矫直过程中重轨与矫直辊不断接触、相对滑动、脱离, 为了真实反映这种变化过程, 必须使计算增量步减小到一定程度, 这样就涉及到运动方程对时间和空间的积分问题。根据时间积分方法的不同, 采用动力显式积分算法与经典的静力隐式平衡迭代求解方法相比, 动力显式积分算法具有计算内存少、速度快、时间步长小, 并无需定义接触单元, 使得三维接触处理简单, 同时无需求解刚度矩阵, 计算效率高, 收敛性好的优点。本论文采用了动力显式积分算法。

动力学的基本方程:

平衡方程:undefined域内) (2)

几何方程:undefined域内) (3)

物理方程:σij=Dijklεkl (V域内) (4)

边界条件:undefined (在Su边界上) (5)

undefined (在SF边界上) (6)

式中 ρ——质量密度;

μ——阻尼系数;

fi——单位体积力;

undefinedi、undefinedi——i方向的加速度和速度;

undefinedi——面力载荷;

nj——现时构形边界SF的外法线方向余弦;

ρundefinedi和μundefinedi——分别代表惯性力和阻尼力。

4 重轨矫直过程的数值模拟

在进行重轨模拟矫直的研究过程中, 考虑到有限元网格划分的特点及计算机的计算能力, 把模型进行简化, 轨底两端和轨头下颚两侧面的小圆弧简化为折线, 大圆弧简化为直线段。

建立水平辊辊式矫直机的模型为:上下各5个矫直辊, 即为10辊水平辊辊式矫直机。第10号矫直辊作为辅助辊, 没有压下量。为了减少单元数目, 节省计算时间, 把矫直辊简化成环形, 重轨矫直的有限元模型如图2所示。根据实际情况, 矫直辊定义成刚性体, 重轨定义成变形体。矫直辊定义成刚性接触体可以大大减少显式分析的时间, 这是由于刚性体内所有节点的自由度都耦合到刚性体的质量中心上去。

根据矫直机及被矫重轨的实际情况, 确定了矫直机的基本参数、矫直辊压下量参数。见表1:

5 模拟矫直后重轨纵向残余应力的大小及分布

重轨矫直过程的开始阶段, 进入矫直盲区, 没有进入正常的矫直阶段, 随着重轨的咬入长度增加, 开始进入正常的矫直阶段。在模拟矫直时, 下排辊无压下量, 上排辊有压下量, 且压下量越来越小。重轨在第三辊处弯曲程度最大, 其弹塑性变形量也最大。随着重轨向前运动, 上排辊的压下量逐渐减小, 重轨所受的弯曲程度也逐渐减小, 弹塑性变形量也逐渐减小。当重轨矫到第七辊时, 重轨开始进入脱离矫直阶段, 即结束矫直过程。在整个矫直过程中, 矫直轨的压下量依次减小, 因此, 对矫直过程中重轨内部的应力变化的影响逐渐减小。分别截取距轨端0.9 m、1 m、1.1 m、1.3 m、1.5 m和1.6 m处六个截面, 分析其纵向残余应力的大小及分布情况, 纵向残余应力云图如图4所示。

由图4、图5可知, 重轨纵向残余应力在整个截面上的分布为:轨头表面与轨底表面中央为正值, 轨底两边缘为负值, 轨腰处为正负相间。这种分布情况与文献[6]所描述的情形基本上是一致的, 只是残余应力值的大小略有不同。最大拉应力出现在重轨的轨头表面与轨底表面中央部位, 其值为157 MPa, 最大压应力出现在重轨的轨腰部位, 其值为139 MPa。EN规定重轨纵向残余应力不得超过250 MPa, 我国铁道部颁发的《时速200 km客运专线60 kg/m重轨暂行技术条件》中也规定, 轨底最大纵向残余拉应力应小于等于250 MPa。模拟矫直后重轨残余应力值是符合标准的。

6 实验分析

实验是用应变片组桥法, 如图6所示。

测试样轨选取某厂复合矫直后的重轨来进行试验。测定部位选在样轨轨底中心线附近, 沿纵向贴一电阻应变片, 并将温度补偿片贴在同材质的补偿板上, 测定部位详见图7。在《时速200 km客运专线60 kg/m重轨暂行技术条件》中规定了残余应力的评定方法, 首先, 把电阻应变片贴在轨底表面, 然后将贴有应变片的部分与重轨逐渐切割隔离, 用释放的应变值来评定残余应力。贴应变片的轨底表面处理和应变片使用方法均按要求进行。

残余应力的计算公式如式 (7) :

σ=-εE (7)

式 (7) 中 σ——残余应力值, MPa;

ε——释放应变值, με即10-6;

E——弹性模量, 2.07e5MPa。

测试结果:U71Mn重轨轨底中心纵向残余拉应力为265 MPa、246 MPa、246 MPa、241 MPa、202 MPa, 均值为240 MPa;U75V重轨轨底中心纵向残余拉应力为236 MPa、237 MPa、242 MPa、242 MPa、239 MPa, 均值为239 MPa。

7 误差的分析

从图3中分析可以得出水平模拟矫直后重轨轨底表面处的纵向残余应力为157 MPa, 实际测量结果的平均残余应力值分别为240 MPa、239 MPa, 之间存在80 MPa左右的误差。出现的误差有以下几方面原因:

(1) 在建立重轨水平矫直的模型时, 是把模型进行了简化, 重轨的圆弧全部用小线段组成的多条边代替, 这与实际就有了差异。

(2) 在数值模拟过程中, 把矫直辊定义成多边形刚性体, 矫直辊与重轨的接触是不连续的, 而实际的矫直过程是连续的接触过程, 这种差异也影响了残余应力的分布和大小情况。

(3) 实验所采用的重轨是经过复合矫直的重轨, 其残余应力与论文中单一水平辊矫直后得到的重轨残余应力之间必然存在差别, 这是误差产生的主要原因。

(4) 在测量时, 由于测量仪器本身有系统误差, 这也对测试结果有一定的影响。

8 结论

铁路运输在我国经济发展中起着举足轻重的作用, 它的发展直接关系到我国国民经济的发展。重轨的残余应力是重轨成材率的重要指标之一, 研究重轨的残余应力有助于提高重轨成才率和重轨质量。通过对重轨矫直模型的建立、研究和分析, 得出其残余应力分布规律。

摘要：铁路运输在我国经济发展中起着举足轻重的作用, 它的发展直接关系到我国国民经济的发展。重轨的残余应力是重轨成材率的重要指标之一, 研究重轨的残余应力有助于提高重轨成才率和重轨质量。通过对重轨矫直模型的建立、研究和分析, 得出其残余应力分布规律。

关键词：重轨,矫直,残余应力,有限元模型

参考文献

(1) 郭华, 陈亚平, 等.矫直引起的钢轨头部变形分析 (J) .钢铁钒钛, 2000, 21 (4) :28-33.

(2) 林丽华, 陈立功, 等.残余应力测量技术及其发展动向 (J) .机械, 1998, 25 (5) :46-49.

(3) Finstermann G, Fischer FD, Shan G, Schleinzer G.Residual stresses in rails due to roll straightening.Steel Research1998;69:272-8.

(4) 杨海波, 汪家才, 王卫平.重轨矫直过程应力应变模型的确定与分析 (J) .北京科技大学学报, 1997.

(5) 邹家祥.轧钢机械 (M) .北京:冶金工业出版社, 2000.

应力模拟 篇7

人类工程活动作为危岩崩塌灾害产生的主要诱发因素之一, 应力场的变化直接反应了危岩崩塌灾害的形成过程, 对危岩的形成机理研究具有直接的指导意义。陈洪凯和唐红梅建立了危岩稳定性的计算方法[1], 并且运用断裂力学方法研究了危岩主控结构面在临界条件下的疲劳断裂寿命[2];姜克春利用数值模拟的方法, 从应力场变化角度, 对危岩链式规律发育的力学演绎过程进行了研究[3];郑安兴等人运用扩展有限元方法, 以重庆万州太白岩危岩为例, 分析了危岩主控结构面的断裂扩展行为[4];张瑞刚针对压剪破坏型危岩, 采用突变理论和功能理论对危岩突发性破坏机理进行了研究[5];李明等人基于平面应变假定, 采用有限元数值模拟软件对重庆洪崖洞危岩的发育过程进行了模拟计算[6];王根龙等人针对塑流-拉裂式崩塌灾害, 利用离散元模拟得到崩塌的变形破坏地质过程[7]。

本文以重庆市合川区三汇镇磨子岩W2危岩为例, 利用ansys有限元数值模拟软件, 分析其在采矿活动下的应力场变化。W2危岩作为磨子岩的代表性危岩体, 位于危岩带的中部, 体积约27.1万m3, 形状类似于板状, 其破坏方式为滑移-剪切式破坏。根据现场调查及资料分析, 危岩坡面近于直立, 岩层产状125°∠20°, 为反向坡。危岩体及基岩的组成岩性主要为三叠系下统飞仙关组 (T1f) 及二叠系上统长兴组 (P2c) , 磨子岩W2危岩正面图及地质剖面示意图如图1、图2所示。

1 磨子岩初始应力场有限元分析

针对采矿活动对山体及危岩应力场的影响, 以下以磨子岩W2危岩体作为物理原型作二维数值模拟。模型模拟材料为灰岩, 弹性模量为3×104MPa, 泊松比为0.25, 在暴雨条件下, 岩体的饱和容重为26.77kn/m3, 施加位移边界约束条件, 采用平面四节点四边形单元PLAN42将整个模型划分为2001个单元, 划分单元后的有限元模型如图3所示。

1.1 初始状态下陡崖的应力场分布

初始应力状态下, 受构造运动及河流下蚀切割作用, 陡崖顶部裂隙发育, 危岩体顶部与母岩脱离, 裂隙张开宽度达120-450cm, 陡崖内部岩体主应力矢量图如图4所示。由图可以看出:

(1) 数值模拟分析结果符合岩体内应力场的分布规律, 即越靠近临空面, 大主应力越接近平行于临空面, 小主应力则与临空面正交, 于表面处降为0, 远离临空面, 地应力逐渐恢复原始状态。

(2) 分析初始状态下的主应力分布规律可以看出, 陡崖内部岩体主应力从危岩体位置向陡崖内部呈发散型分布。危岩体内主应力较小, 大主应力从上到下逐渐减小, 从内至外逐渐减小;而小主应力则呈现出从上至下逐渐减小, 从内至外逐渐增大的规律。

(3) 在陡崖顶部, 岩体受到拉应力, 应力平行于临空面, 故易产生张拉裂缝。在陡崖侧面, 岩体受到挤压, 以压应力为主, 从上到下, 压应力逐渐增大, 拉应力基本为0。

图5、图6为初始状态下磨子岩W2危岩体陡崖大、小主应力等值线图, 分析陡崖内部主应力的分布情况可以看到, 岩体内大主应力全部为拉应力, 危岩体附近及陡崖底部大主应力较小, 陡崖顶部拉应力较大, 最大值为6.50MPa, 从陡崖内部向外呈现逐渐增大的趋势;岩体内小主应力均为压应力, 从山体表面到内部, 岩体受到的压应力逐渐增大, 临空面附近及陡崖底脚压应力较小, 山体内部岩体最大压应力为6.65MPa。对比大、小主应力可以看出, 山体内部岩体最大、最小主应力大小基本相等。

1.2 初始状态下危岩体应力场分布情况

分析初始应力状态下W2危岩体内部的应力分布情况, 其大小主应力等值线图如图7、图8所示, 由图可以看出, 危岩体后部及底脚处大主应力最大, 拉应力最大, 且岩体临空面附近较内部大, 变化范围在0-0.018MPa;危岩体内小主应力全为压应力, 危岩体后缘裂缝尖端及底脚压应力较大, 底部压应力为-2.95 MPa。危岩体内部岩体压应力上部较下部小, 变化范围在-38-0MPa之间。

2 采矿活动下磨子岩应力场演化规律分析

采矿活动下, 岩体结构体发生改变, 应力重分布, 巷道危岩释放潜能, 危岩向巷道空间移动。将采矿活动对磨子岩危岩稳定性的影响当做平面问题来考虑, 即不考虑煤柱对采空区的支撑作用, 分析采矿活动下危岩应力场的变化。

2.1 采矿活动下陡崖应力场分布情况

分析某种开挖条件下陡崖内部岩体主应力矢量图, 如图9所示, 采矿活动使山体内岩体应力方向发生偏转。巷道的采空形成新的临空面, 临空面附近危岩小主应力与临空面平行, 巷道开挖工作面应力集中, 容易发生巷道的垮塌。陡崖内部, 靠近采空区附近围岩以压应力为主, 陡崖后部远离采空区的区域应力方向无明显变化。且根据箭头长短可以看出, 危岩体内应力较小。

分析采矿活动下山体的应力等值线图10、图11, 对比初始条件下的应力场, 分析采矿活动对陡崖的影响。由应力等值线图可以看出:

(1) 采空区上覆围岩以及开挖工作面应力明显增大, 形成应力集中带。其中采空区上覆围岩拉应力最大, 最大拉应力达到8.71MPa;危岩体所受拉应力较小。

(2) 采矿活动对山体应力场的影响主要体现在小主应力方面。采空区上覆围岩主要受压应力的影响, 在采矿活动的影响下, 小主应力变化明显。模型所示开采程度下, 最大压应力达-34.8MPa, 较初始状态下的-6.65MPa, 压应力明显增加, 其中, 开挖工作面所受压应力最大。

2.2 采矿活动下危岩体在应力场分布

如图12、图13危岩体应力等值线图所示, 分析W2危岩体内岩体在采矿作用下应力的变化, 可以得出以下结论:

(1) 采矿活动下, 危岩体内大主应力变化不明显, 较陡崖内部岩体所受拉应力较小。该开挖条件下, 危岩体内大主应力在0-0.097MPa之间, 较初始条件下0-0.018MPa有所增大。由图还可以看出, 岩体后部裂缝尖端拉应力集中, 有扩展趋势, 由于灰岩抗拉强度较小, 当尖端所受拉应力超过其抗拉强度, 随着开挖的进行, 裂缝逐渐加宽加深并贯通成主控结构面, 最终发生岩体崩塌破坏。

(2) 危岩体内小主应力变化明显, 在-3.87-0MPa之间, 均为压应力, 较初始状态的-38-0MPa, 所受压应力明显减小。其危岩体小主应力等值线图如图13所示。由图可以看出, 危岩体后部裂缝尖端应力集中, 危岩体内压应力大致呈由临空面向岩体内逐渐递增的趋势。

根据以上应力分析可以看出, 危岩体后部裂缝尖端应力集中明显, 以下将采矿活动分为8个开挖工况, 得出随着工作面的推进, 危岩体后部裂缝的扩展曲线如图14所示。由图可以看出, 裂缝尖端在初始应力条件下受压, 但所受压力不大, 接近于0;随着开挖工作面的推进, 裂缝尖端所受压应力先增大后减小, 最终转变为拉应力, 当开挖至一定程度时, 拉应力的增长随开挖程度接近线性分布, 压应力最大值达到0.920MPa。

3 结论

选取磨子岩W2危岩体作为地质模型, 利用数值模拟软件, 通过对比分析采矿活动前后陡崖及危岩体内主应力矢量图及等值线图, 得出采矿活动对陡崖及危岩应力场的影响, 主要体现在以下几个方面:

(1) 采矿活动形成新的临空面, 临空面附近应力方向发生偏转, 开挖工作面应力集中, 容易发生巷道的垮塌。

(2) 采矿活动对陡崖及危岩体应力场的影响主要体现在小主应力方面, 压应力明显增加。初始状态下, 最大压应力为-6.65MPa, 而模型所示开采条件下, 最大压应力达-34.8MPa。

(3) 危岩体体内大主应力变化不明显, 采矿活动引起危岩体内小主应力明显增大, 压应力明显减小, 由初始状态的-38-0MPa增大至-3.87-0MPa, 并且呈由临空面向内递增的趋势。

(4) 危岩体后部裂缝尖端应力集中, 随着开挖的推进, 裂缝尖端应力呈现先减后增的趋势, 当开挖达到一定深度时, 与开挖深度基本呈线性关系。

参考文献

[1]陈洪凯, 鲜学福, 唐红梅, 等.危岩稳定性分析方法[J].应用力学学报, 2009, 26 (2) :278-282.

[2]陈洪凯, 唐红梅.危岩主控结构面疲劳断裂寿命计算方法[J].应用数学和力学, 2007, 28 (5) :575-580.

[3]姜克春.缓倾角层状岩质边坡破坏规律的力学机理研究[D].重庆交通大学, 2008.

[4]郑安兴, 罗先启, 沈辉.危岩主控结构面变形破坏的扩展有限元法模拟分析[J].岩土力学, 2013, 34 (8) :2371-2377.

[5]张瑞刚.危岩突变性破坏机制研究[D].重庆交通大学, 2012.

[6]李明, 陈洪凯, 叶四桥.重庆市洪崖洞危岩发育机理[J].中国地质灾害与防治学报, 2008, 19 (2) :1-5.

体外预应力效应的有限元模拟方法 篇8

应用有限元方法对体外预应力混凝土桥梁结构进行受力分析和数值模拟, 具有重要的工程应用和学术研究价值, 但目前对有限元解法的研究尚少。Wu等提出了基于平面应力有限元方法的预应力混凝土结构的分离式模型, 能较好的模拟结构应力场及混凝土和钢筋的接触面, 但单元数量较多, 使用不方便。直接采用杆单元模拟外预应力钢索更为简单和实用。文中给出了体外索单元的切线刚度矩阵并将其转换成两端带刚臂体外索单元的刚度矩阵, 在不增加单元和节点的条件下, 直接采用有限元方法较精确地模拟了体外预应力效应及体外索和混凝土梁体之间的传力机理。

1体外预应力混凝土桥梁结构的有限元分析方法

体外预应力混凝土桥梁从最初的简单构件到最终的承重结构, 经历了一个结构体系转换与受力变化的成型过程。体外预应力筋仅在转向构造处与混凝土梁体接触, 组成了一个内部超静定结构体系, 因此结构计算可以采用有限单元分析方法[3,4]。在有限元计算过程中, 必须能反映体外预应力混凝土桥梁结构的受力变形特点。体外预应力钢筋在桥梁结构中只承受拉力, 与一般的杆单元不同。在外荷载作用下, 除转向座外, 体外筋的变形与梁的挠度不一致, 这与有粘结体内预应力钢筋的受力也有很大差别。考虑体外索的这些特点, 体外预应力混凝土桥梁的计算模型见图1。

体外预应力钢筋在转向构造处与混凝土梁上节点的连接以刚臂模拟, 在其余位置与混凝土梁体之间不存在变形协调关系, 这样体外索单元就成为两端带刚臂的只受拉杆单元。张拉体外预应力钢筋后, 体外索对梁施加了力的作用, 当体外索与梁接触处的约束状况确定后, 各段预应力筋的张拉力都随之确定。体外预应力筋的受力作用便自动形成, 同时桥梁结构体系上增加数个两端带刚臂的体外索单元。对体外预应力效应的模拟转换为在结构有限元求解方程中正确地添加体外索单元。将体外预应力钢筋作为构件反映在结构分析模型中, 除能较真实的体现体外预应力筋的作用外, 还能解决体外预应力筋回缩和分批张拉引起的预应力损失、混凝土收缩徐变引起的预应力损失的计算问题。因此, 体外预应力效应有限元模拟的重点问题是正确建立体外索单元的刚度矩阵并在结构分析中将其转换为两端带刚臂单元的刚度矩阵。

2体外预应力索单元的切线刚度矩阵

对体外预应力混凝土桥梁结构, 应考虑受力过程中的几何非线性问题。结构几何非线性采用Updated Lagrangian[6]描述, 即参考构形是基于结构的前一相邻构形。设定t时刻单元应变为0, t+Δt时刻的单元应变为单元的应变增量, Piola-Kirchhoff应力增量即为Cauchy应力增量。于是, 根据虚功原理, 有增量形式的体外预应力索单元的切线刚度矩阵为:

[K]=[K]L+[K]G (1)

[K]L=∫V[B]T[D][B]dV (2)

$[{\rm K}]{}_G={\displaystyle {\int }_V\left[\frac{\partial {\rm N}}{\partial x}\right]}{}^{\text{Τ}}[\sigma ]\left[\frac{\partial {\rm N}}{\partial x}\right]\text{d}V$ (3)

式中:[K]L——常规的弹性刚度矩阵;

[K]G——几何刚度矩阵或初应力刚度矩阵, 其含义是由于结构几何变形引起总体刚度矩阵的修正;

[B], [D], N——应变矩阵、线性本构矩阵和位移形函数。

根据式 (2) 和式 (3) , 不难得出弹性刚度矩阵和几何刚度矩阵的表达式为:

$[{\rm K}]{}_L=\frac{E{}_{pt}A{}_p}{l{}_p}\left[\begin{array}{cccccc}1& 0& 0& -1& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ -1& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]$

(4)

$[{\rm K}]{}_G=\frac{{\rm N}{}_p}{l{}_p}\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& -1& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& -1& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]$

(5)

式中:Ept——体外索单元的切线弹性模量;

Ap——体外索单元的横截面积;

Np, lp——体外索单元的轴向拉力和长度。

矩阵[K]L和[K]G相加, 即得到体外预应力钢筋单元的切线刚度矩阵。

3体外预应力索单元切线刚度矩阵的转换

体外索单元的两端以刚臂和梁上节点联系, 如果专门推演带刚臂杆单元的刚度矩阵, 会使分析求解的计算量增加。本文根据体外预应力钢索和混凝土梁体之间的传力特点, 直接对式 (1) 中的刚度矩阵做一些转换, 得出两端带刚臂的体外索单元刚度矩阵, 可以直接用于桥梁结构的有限元分析计算。这种做法不增加节点和单元的数目, 且便于有限元程序的编制。

节点i和j为原混凝土梁上的节点, 节点a和b为体外预应力钢索上的节点。节点a与b之间为一个体外索单元, 两端的偏心距分别为ea和eb。i点与a点之间以及j点与b点之间均以刚臂连接。坐标系$\overline{X}$$\overline{Y}$为混凝土梁上的局部坐标系;坐标系X′Y′为体外索单元上的局部坐标系;而XY为结构的整体坐标系 (见图2) 。

在坐标系X′Y′中, 得到体外索单元的切线刚度矩阵为[K]′, 记[T]为坐标系$\overline{X}$$\overline{Y}$到X′Y′的转换矩阵, 则在$\overline{X}$$\overline{Y}$坐标系中, 体外预应力筋单元的刚度矩阵表达式为:

$[\overline{{\rm K}}]=[{\rm T}]{}^{\text{Τ}}[{{\rm K}}^{\prime }][{\rm T}]$ (6)

由于刚臂本身没有变形, 故可按不同坐标系表达刚臂连接两点之间的位移关系。基于混凝土梁上的局部坐标系, i点与a点之间的位移关系描述为:

$\{\overline{\delta }{}_a\}=\left\{\begin{array}{l}\overline{u}{}_a\\ \overline{v}{}_a\\ \overline{\theta }{}_a\end{array}\right\}=\left[\begin{array}{ccc}1& 0& e{}_a\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left\{\begin{array}{l}\overline{u}{}_i\\ \overline{v}{}_i\\ \overline{\theta }{}_i\end{array}\right\}=\left[\begin{array}{ccc}1& 0& e{}_a\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\{\overline{\delta }{}_i\}=A{}_a\{\overline{\delta }{}_i\}$

(7)

同理, j点与b点之间的位移关系为:

$\{\overline{\delta }{}_b\}=\left\{\begin{array}{l}\overline{u}{}_b\\ \overline{v}{}_b\\ \overline{\theta }{}_b\end{array}\right\}=\left[\begin{array}{ccc}1& 0& e{}_b\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left\{\begin{array}{l}\overline{u}{}_j\\ \overline{v}{}_j\\ \overline{\theta }{}_j\end{array}\right\}=\left[\begin{array}{ccc}1& 0& e{}_b\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\{\overline{\delta }{}_j\}=A{}_b\{\overline{\delta }{}_j\}$

(8)

因此, 体外索单元两端节点的位移和梁上两节点的位移间的关系由下式确定:

$\left\{\begin{array}{l}\overline{\delta }{}_a\\ \overline{\delta }{}_b\end{array}\right\}=\left[\begin{array}{cc}A{}_a& 0\\ 0& A{}_b\end{array}\right]\left\{\begin{array}{l}\overline{\delta }{}_i\\ \overline{\delta }{}_j\end{array}\right\}=[A]\left\{\begin{array}{l}\overline{\delta }{}_i\\ \overline{\delta }{}_j\end{array}\right\}$

(9)

若以$\{\overline{S}{}_p\}$表示体外索节点a与b上的杆端力, $\{\overline{S}{}_p\}$等效移至梁上节点i和j, 得到的杆端力为$\{\overline{S}{}_L\}$, 由虚功原理可得它们之间的关系为:

$\{\overline{S}{}_L\}=[A]{}^{\text{Τ}}\{\overline{S}{}_p\}$ (10)

在坐标系$\overline{X}$$\overline{Y}$中, 体外预应力筋单元的杆端力和节点位移之间的关系为:

$\{\overline{S}{}_p\}=[\overline{{\rm K}}]\left\{\begin{array}{l}\overline{\delta }{}_a\\ \overline{\delta }{}_b\end{array}\right\}$

(11)

将式 (9) 和式 (10) 代入式 (11) , 得到:

$\{\overline{S}{}_L\}=[A]{}^{\text{Τ}}[\overline{{\rm K}}][A]\left\{\begin{array}{l}\overline{\delta }{}_i\\ \overline{\delta }{}_j\end{array}\right\}=[\overline{{\rm K}}]{}_E\left\{\begin{array}{l}\overline{\delta }{}_i\\ \overline{\delta }{}_j\end{array}\right\}$

(12)

其中, $[\overline{{\rm K}}]{}_E$即为两端带有刚臂的体外预应力索单元在局部坐标系$\overline{X}$$\overline{Y}$中的单元刚度矩阵。

将$[\overline{{\rm K}}]{}_E$从局部坐标系$\overline{X}$$\overline{Y}$转换到整体坐标系XY, 需乘以转换矩阵[L]。

$[{\rm K}]{}_E=[L]{}^{\text{Τ}}[\overline{{\rm K}}]{}_E[L]=[L]{}^{\text{Τ}}[A]{}^{\text{Τ}}[\overline{{\rm K}}][A][L]$ (13)

式 (13) 表明, 两端带刚臂的体外索单元刚度矩阵可以由一般体外索单元刚度矩阵乘以转换矩阵得到。通过矩阵[K]E, 使体外预应力筋单元直接以构件的形式参与桥梁结构的有限元分析计算。

4结语

根据体外预应力钢筋和混凝土梁体之间的连接形式和传力特点, 得出两端带刚臂体外索单元的刚度矩阵可由一般体外索单元刚度矩阵转换而来。使得体外预应力钢筋作为构件反映在结构分析计算中, 方便地实现了体外预应力效应的有限元模拟。

摘要：对体外预应力钢索的传力特点进行了研究, 推导了体外索单元刚度矩阵在桥梁有限元分析中的转换公式, 从而得出两端带刚臂体外索单元的刚度矩阵, 实现了体外预应力效应的有限元模拟。

关键词：体外预应力,刚度矩阵,刚臂,有限元

参考文献

[1]Ehsani M, Saadatmanesh H.Behavior of externally prestressedconcrete griders[J].Struct Des, Anal Test, 1989 (8) :218-222.

[2]Fenves G.Nonlinear analysis of external prestressed bridegs[J].ASCE, 1986 (13) :192-201.

[3]李国平.预应力混凝土结构设计原理[M].北京:人民交通出版社, 2000.

[4]Lou T-J.Finite element modeling of concrete beams prestres-sed with external tendons[J].Engineering Structures, 2006 (28) :1919-1926.

[5]林冰花, 谭德盼.体外预应力技术的发展[J].山西建筑, 2007, 33 (8) :149-150.

应力模拟 篇9

1 试验系统

研究所采用的WYQ1000-I型地下工程综合模拟试验系统,是由中国矿业大学建工学院设计研制的。该试验系统具有4个通道(1个备用通道),采用液压加载、稳压、6面油缸组合加载的方法,可以完成对三维岩体的载荷控制,有效克服以往平面应变加载模型中存在的问题。

2 模型开挖与测试

2.1 模型初始应力条件加载

试验以泥岩围岩的全煤巷道为基础,以水平产状的岩层节理为围岩的裂隙分布形式,巷道的埋深为变量,在800 m埋深巷道的基础上,每次增加200 m埋深的荷载增量,依次模拟埋深800 m、1 000 m、1 200 m、1 400 m、1 600 m的巷道;将计算出的顶部荷载和侧压分成均等6级,逐步加载。为保证试件受力均匀,每次加载后稳定0.5 h,直至试验初始设计值达到0.476 MPa,相当于800 m埋深时,垂直压力应达到20 MPa,最大水平主应力为22 MPa,对应初始设计值应为0.524 Pa(水平主应力与垂直主应力的相关系数在深部地应力为1.1左右),最小水平主应力为8.6 Pa(0.204 MPa)。加到此预定荷载后稳定各加载板油压,保持荷载不变。与此同时,在每次施加各级荷载过程中注意观察模型变形和油缸活塞位移情况,并记录下每个加载面的加载压力。

根据试验系统及加压系统的技术指标,可知各荷载板设计载荷所对应的油压读数,如表1所示。

上下板为0.476/12.35×23.8=0.92 MPa; 前后板为0.204/12.35×26.7=0.44 MPa; 左右板为0.524/12.35×23.8=1.01 MPa。

2.2 开挖与加载

在模型上下、左右、前后利用液压缸分别控制加围压至设计荷载后,稳定一段时间,拆除模型的后荷载板中间的挡板,在其中心位置机械钻孔开洞以模拟巷道开挖。

具体开挖步骤为:模拟开采从模型前荷载板开始,向后荷载板推进,每4 h开挖30 cm(相当于实际每天推进9 m),直至整个巷道模型开挖完成。开挖后,根据试验安排,分别模拟无支护与气囊支护条件下的巷道受力变形破坏情况。

1)无支护状态,巷道开挖完成后,不进行任何支护阻力的设计,在初始地应力下保持其稳定24 h(实际相当于5 d);然后模拟不同埋深巷道载荷,使垂直方向的荷载每隔1.5 h增加0.25倍,分别模拟荷载集中系数K=1.25～2,与此同时做好巷道各阶段的变形观测以及变形破坏情况的记录。

2)气囊支护状态。整个过程与无支护状态相似,不同处有两点:一是巷道开挖完成后立即装上预先加工好的气囊,施加p=10 KPa的支护阻力(相当于实际0.4 MPa);二是与巷道垂直的水平方向荷载每隔1.5 h增加0.25倍,分别模拟侧压系数K=1.25～2.50情况下巷道围岩的变形破坏发展情况。

2.3 试验数据分析

2.3.1 应变分析

1)无支护条件下。 图1为模型铺设完成直至模型最终破坏的全过程应力变化曲线,图1中应变砖的埋设位置分别距开挖自由面3 cm。

从图1可以看出,无支护巷道在加载过程以及巷道变形破坏期间,最大最小主应力的变化趋势基本一致:①加载的初期,随着模型的压实其变化逐渐趋于一致。②巷道开挖初始阶段,由于巷道开挖打破了原有三向应力平衡状态,卸荷后随即发生内应力释放和回弹,引发巷道各部位应力调整和变形,其中底板应力下降约0.3 MPa,帮部约0.4 MPa,顶板约0.6 MPa;由此可见顶板围岩体卸荷所引起的应力差值最大,应力调整阶段所引起的顶板变形破坏会更为强烈。③应力调整阶段,底板与两帮的应力集中程度较低,由于得不到有效约束,应力持续减小,继续加载,底板由2 MPa下降为1 MPa,帮部由3 MPa降至1.5 MPa,下降幅度均达50%;而顶板的应力集中程度最高,巷道顶板首先产生破坏,随后拱顶上方出现较大的变形乃至破坏,这一过程中,顶板的应力变化则有两峰值点、一低谷段,经历了“集中—卸荷—再集中—再卸荷”的过程,说明顶板在冒落过程中会自动寻找新的压力平衡拱各图中17～18对应的时间段,应力会由0.8 MPa上升至0.9 MPa,但这种平衡较为短暂,随即伴随着外围承载结构的破坏和压力的传递,最后破坏区不断向深部扩展,巷道失稳破坏。

2)有支护条件下。图2为模型铺设完成直至模型最终破坏的全过程应力变化曲线,图2中各部位应变砖的埋设位置分别距开挖自由面3 cm。

从图2中可以看出,巷道在加载过程以及巷道开挖后的变形破坏期间,最大最小主应力的变化趋势是基本一致的:①加载的初期,与无支护巷道变化规律相同,随着模型的压实变化趋于一致。②巷道开挖阶段,巷道各部位围岩的应力都存在与之相对应的卸荷区(各图中13～18所对应的时间段);其中底板下降约为0.2 MPa,帮部约为0.6 MPa,顶板约为0.3 MPa,与无支护状态下的模拟试验相比可知,最大应力差值发生在最大主应力方向,该方向应是巷道维护的重点。③应力调整阶段,当支护施加后,围岩体应力在重新调整中,巷道周边岩体应力升高,应力环境产生了较为明显的转变,其中底和两帮表现最为明显,底板应力由1.7 MPa提升到1.9 MPa,帮部应力由2.1 MPa提升到2.5 MPa,顶板也由卸载后的0.8 MPa提升到1.0 MPa;0.01 MPa的支护力提升了0.4 MPa的围岩应力,说明巷道围岩开挖卸荷后,支护对岩体能提供一定的支护抗力,使破碎区围岩保持三向受压,寻求到一个新的平衡状态(18～20所对应的时间段)而不至于持续向深部发展;当碎胀作用力超过支护阻力时(20～22所对应的时间段),支护体系让压新的破坏区域随之产生,但浅部围岩在支护力作用下没有发生失稳冒落,而向深部围岩发展的破坏区域很快在新的区域寻求到平衡点,因此可以看出支护对围岩二次稳定结构的形成起着至关重要的作用。

2.3.2 位移及破坏特征分析

模型在不同支护条件下的图像采集效果如图3所示。

3 结 论

通过对有支护和无支护两种情况下巷道围岩应力及位移的分析可知,施加一定阻力的柔性支护能使围岩体处于良好的工作状态,使两种承载结构都能充分发挥作用,支护作用机理主要有以下4点:

1)延缓和减少围岩中裂纹的产生、扩展与贯通,防止已破碎围岩塌落。

2)加快围岩承载结构的形成和稳定。由于支护力的存在使围岩更大程度的表现出三向受压状态下的强度特征,此种情况有利于围岩承载状态的形成。

3)增强了松动区域岩石强度。无支护情况下会导致围岩强度的持续下降而最终诱发冒顶、片帮等现象,支护力的存在使得破碎区围岩能维持一定的残余应力,有效控制了破碎区以及塑性区半径,使更大范围的围岩得以处于弹性变形状态,加快了松动圈的产生和稳定过程,这对保持巷道的稳定更为有利。

4)有利于围岩二次稳定形成,具体来说,就是要加固和维护好浅层围岩的张拉破坏区域,充分发挥和利用应力集中压缩区域的承载能力,加速围岩二次稳定过程的形成。

综上所述,与巷道无支护情况下相对比可以看出,破碎围岩体由于得到气囊支护的有效约束,围岩表面张拉破碎区的残余强度虽然较低,但在径向约束条件下,一定的支护阻力阻止了碎胀变形的持续发展,支护阻力通过控制破碎区范围实现对围岩变形的控制作用,改善了围岩内部的受力环境和分布特征,从而不致诱使巷道围岩的变形失稳。

摘要：随着我国浅部资源的逐渐减少和枯竭,矿物资源地下开采的深度越来越大,深埋高应力巷道的围岩变形、破坏规律的研究显得更加重要。采用WYQ1000-I型地下工程综合模拟试验系统,通过三维物理试验模拟巷道的真实开挖、变形和破坏过程,研究有、无支护体系形态下的围岩变形,针对围岩应力以及位移变形破坏特征进行深入分析,得出深埋高应力巷道围岩变形破坏及支护作用机理。

关键词：深埋高应力,三维物理模拟,变形机理

参考文献

[1]靖洪文.深部巷道大松动圈位移分析及应用[M].徐州:中国矿业大学出版社,2001.

[2]孙晓明,何满潮.深部开采软岩巷道耦合支护数值模拟研究[J].中国矿业大学学报,2005,(2):166-169.

[3]孔德森.深部巷道围岩在复合应力场中的稳定性数值模拟分析[J].山东科技大学学报,2001,21(1):68-70.

[4]朱德仁,王金华,康红普,等.巷道煤帮稳定性相似材料模拟试验研究[J].煤炭学报,1998,23(1):42-47.

[5]胡耀青,赵阳升,杨栋,等.带压开采顶板破坏规律的三维相似模拟研究[J].岩石力学与工程学报,2003,22(8):1239-1243.

[6]李世平,吴振业,贺永年,等.岩石力学简明教程[M].北京:煤炭工业出版社,1996.