《工程力学数值方法》课程作业.doc

1、工程力学数值方法工程力学数值方法课程作业目录1 问题概述12 圆孔孔口应力集中问题的理论计算23 圆孔孔口应力集中的ANSYS有限元分析43.1 建立有限元模型43.2 施加荷载并求解43.3 查看分析结果并验证结果的正确性43.4 用理论解对数值解进行验证并得出结论74 圆孔孔口应力集中的边界元分析84.1边界单元间接法的基本原理84.2边界单元间接法的数值方法94.3 孔口应力集中问题115 总结131 问题概述在工程结构设计中，常常会根据需要在构件上进行开孔，由于孔口的存在势必会破坏金属材料的连续性，引起主应力弯曲绕行，应力作用线在孔口附近密集弯曲，造成局部应力增大的现象。由圣维南原理可

2、知，在远离小孔的地方，孔口局部应力集中的影响将消失。孔口附近的应力远大于无孔时的应力，也远大于距孔口较远处的应力，这种现象称为孔口应力集中。应力集中是由于开口后的应力扰动引起的, 最大与最小应力一般发生在孔边上，且应力扰动主要发生在距孔边1.5倍孔口尺寸的范围内。孔口应力集中与孔口形状有关，圆形孔口的应力集中程度比较低，故在工程上被广泛使用。在此，研究小孔口问题，即孔口的尺寸远小于弹性体的尺寸，并且孔边距弹性体的边界比较远(约大于1.5倍的孔口尺寸)，以排除弹性体边界条件对孔口应力分布的影响。孔口应力集中的存在会使材料的脆性加大，构件的强度降低，尤其会严重影响交变应力下构件的持久疲劳极限。因此

3、，有必要对孔口应力集中问题进行深入的研究和探讨，为工程施工提供理论指导和依据。本文首先从理论上对孔口应力集中问题进行分析，然后利用有限元软件ANSYS14.5该问题进行模拟分析,并用理论解对ANSYS计算的数值解进行验证,并将理论解与数值解进行分析比较,总结出孔口应力集中问题的特征。122 圆孔孔口应力集中问题的理论计算图2-1所示一矩形薄板只在左右两边受有均布的拉力，板中有一圆孔，半径为r，矩形薄板的尺寸远大于孔口尺寸，且孔边距矩形板边界较远。 图2-1 带圆孔的矩形薄板 可按“小孔口问题”求解，采用极坐标求解较为方便，坐标原点取在圆孔中心，利用弹性力学的方法可得到该矩形薄板应力状态的基尔斯

4、解答如下:为了研究孔口应力集中情况，在此分情况进行讨论：（1）沿孔边()的环向应力为： 而径向应力，其分布如图2-2（a）所示。（2）沿y轴() 的横截面面上的环向应力为： 其分布规律如图2-2（b）所示。（3）沿x轴() 纵向截面上的环向应力为：其分布规律如图2-2（b）所示。（4）孔边最大应力出现在处，即：最小应力出现在处，即： 小圆孔的应力集中系数为： (a) (b)图2-2 孔口附近应力分布3 圆孔孔口应力集中的ANSYS有限元分析ANSYS是一种大型的通用有限元分析软件，软件主要包括三个部分：前处理模块、分析计算模块和后处理模块。假设矩形板长为400mm，宽为400mm，圆孔半径为2

5、0mm，材质为钢板，采用线性静力分析方法。利用大型通用有限元软件ANSYS14.5进行分析计算，其具体步骤如下:3.1 建立有限元模型(1) 定义单元类型。采用的单元类型为Plane82。(2) 生成几何模型。由于结构和荷载的对称性，取1/4结构进行分析计算，建立一个长为200mm，宽为200mm的正方形，以左小角为圆心建立半径为20mm的1/4圆，进行减运算得到有限元分析模型。(3) 定义材料属性。确定薄板的弹性模量E=2.06105MPa，泊松比为0.3，密度7850kg/m3。(4) 划分网格。每条直线边平均划分为10段，圆弧段平均划分为20段。分别采用三角形单元与四边形单元进行网格划分

6、，如图3-1所示。3.2 施加荷载并求解(1) 定义约束并施加荷载。固定薄板左侧边线，约束左侧边线上所有节点的x方向自由度；在薄板右侧边施加均布荷载，荷载大小为100N/m。(2) 设置分析选项并求解。选择分析类型为静力分析，并进行求解。 a 分析模型 b 三角形单元 c 四边形单元3-1 1/4薄板分析模型及有限元网格划分模型3.3 查看分析结果并验证结果的正确性划分单元时采用了三角形和四边形单元两种划分方法，相应产生两种分析结果，括号前为三角形单元的分析结果，括号内为四边形单元的分析结果。(1) 1/4薄板变形。ANSYS得出的薄板的整体变形中，圆孔沿x方向被拉成椭圆孔，且最大变形为0.1

7、02093mm（0.102132mm）。(2) 1/4薄板x方向的应力云图。ANSYS得出的x方向的最大应力为297.345MPa（254.996MPa），x方向的最小应力为-11.298MPa（-21.776MPa）。1/4椭圆孔左右边缘各有一段压应力区，且沿x轴方向应力由压应力变为拉应力，并逐渐增大到125.877MPa（101.234MPa）。1/4椭圆孔上下边缘均为拉应力，且沿y 轴方向拉应力逐渐由297.345MPa（254.996MPa）迅速减小到125.877MPa（101.234MPa）。可见，ANSYS计算得出的x方向的应力分布与理论计算情况相吻合。(3) 1/4薄板y方向的

8、应力云图。y方向的最大应力为55.244MPa（46.575MPa），y方向的最小应力为-97.425MPa（-70.025MPa）。1/4椭圆孔左右边缘有一段的压应力区，沿着x轴方向压应力急剧变成拉应力，y轴方向拉应力逐渐减小。 a 变形图（三角形单元） b 变形图（四边形单元） c x方向应力云图（三角形单元） d x方向应力云图（四边形单元） e y方向应力云图（三角形单元） f y方向应力云图（四边形单元） g 单元应力云图（三角形单元） h 单元应力云图（四边形单元）图3-2 有限元分析结果 (4) 沿1/4圆孔边环向正应力的ANSYS计算数据如下表3-1： 表3-1 孔边环向正应力

9、的ANSYS计算数据0o30o45o60o95o(MPa)-11.29（-21.77）0.34（0.98）10.33（11.43）21.08（30.21）38.28 （46.57） (5) 沿x轴(从圆孔边缘沿x轴正方向的路径)的环向正应力：在孔边时环向正应力为11.29MPa (-21.77MPa)，沿着x轴方向环向压应力急剧减小为0，随后为环向拉应力逐渐增大，当增大到一定值125.88MPa(101.23MPa)时，又逐渐减小到91.58MPa (70.48MPa)。相应的ANSYS 计算数据如下表3-2：表3-2 沿x轴环向正应力的ANSYS计算数据（mm）r1.5r2r2.5r3r4r

10、5r203040506080100(MPa)-11.29（-21.77）22.30（8.98）40.142（39.73）57.29（54.82）64.13（70.48）125.88（101.23）91.58（70.48） (6) 沿y轴(从圆孔边缘沿y轴正方向的路径)的环向正应力：在孔边为38.28MPa（46.57MPa），沿着y轴方向环向正应力先是急剧减小，然后逐渐减少到4.36MPa（3.21MPa）。相应的ANSYS 计算数据如下表3-3：表3-3 沿y轴环向正应力的ANSYS计算数据（mm）r1.5r2r2.5r3r4r5r203040506080100(MPa)38.28（46.5

11、7）32.66（36.32）26.99（33.62）21.318（20.66）15.32（14.20）9.86（7.70）4.36（3.21）3.4 用理论解对数值解进行验证并得出结论 (1) 沿着孔边的环向正应力: 虽然理论解与数值解之间存在一定的误差，主要是由于ANSYS建模时的诸多因素如网格密度约束条件简化等引起的，但总的来说沿孔边的环向正应力的理论解和ANSYS计算的数值解拟合的很好。三角形单元：从0到90沿孔边应力从最小的11.29MPa骤增到最大的38.28MPa，这是由于应力集中所引起的，应力集中系数为：38.28/11.29=3.39，接近于3.0。四边形单元：从0到90沿孔边

12、应力从最小的21.77MPa骤增到最大的46.57MPa，应力集中系数：46.57/21.77=2.1，误差较大。数据表明三角形单元比四边形单元更加接近理论值。如果不考虑应力集中，则孔边应力应为均匀的100MPa。(2) 沿着y轴的环向正应力:理论解与数值解沿y轴的环向正应力的理论解与ANSYS计算得出数值解拟合的很好，应力集中主要发生在距孔边1.5倍孔口尺寸的范围之内。(3) 沿着x轴的环向正应力:虽然理论解与数值解之间存在一定误差，主要是由于ANSYS建模时的诸多因素如网格密度约束条件简化等引起的，但总的说来沿x轴的环向正应力的理论解与ANSYS计算的数值解拟合的很好。孔边应力为11.29

13、MPa（21.77MPa），但随着远离孔边急剧增大到一定值125.88MPa（101.23MPa），随后逐渐减小到趋近于0MPa。4 圆孔孔口应力集中的边界元分析有限元法计算能使孔口应力集中问题得到较好的解决。近年来，新发展的边界元法只需对边界划分单元，比有限元法输入数据少，过程简便，精度高，总体方程系数矩阵阶数低，是一种强有力的数值计算工具。边界元法可分为直接法、半直接法和间接法，间接法不仅具有以上优点，而且更为简单，易于推广。因此，本文采用边界元法解孔口应力集中问题。4.1边界单元间接法的基本原理假设无限域中一虚拟体I，它与要求解的各向同性均质体I具有相同的几何条件和边界条件，无限域及虚拟

14、体与实体I的材料相同，都是均质的。间接法解弹性力学问题是从Navier方程出发，利用其基本解，此处为开尔文解，以及弹性力学物理方程和几何方程建立虚拟体边界上分布的虚拟力强度和边界面力、位移之间的关系，即边界积分方程，在已知边界条件下，确定虚拟力强度，解虚拟体，从而达到解实体I的目的。Navier方程： （4.1）的基本解，即开尔文解，提供了无限大域内P点的力向量在R点引起的位移： （4.2）利用弹性力学几何方程和物理方程可以获得内力 （4.3）如果Ns(R)是过虚拟体内点R某一斜面S的外法线方向余弦，则此斜面上的应力: 令 则 （4.4）假定边界S上有虚拟分布力则 把内点R移向边界点Q时 （4

15、.5）上式为不考虑体积力时虚拟体的边界积分方程。在一个适定问题中，边界上一部分为已知，另一部分也为已知，从积分方程（4.5）中求得，那么任意点R的应力及位移可以从下式求得 （4.6） 4.2边界单元间接法的数值方法（1） 线性代数方程组的建立计算中，均不考虑体力。将边界离散，另在整个单元内为常数，选取单元中点为节点，则（4.5）为： （4.7） 其中 位移系数 图4-1 边界离散 应力系数 取Q为单元结点，按照给定的边界条件，由（4.7）可以得到线性代数方程组： （4.8）（2） 求位移和应力系数用精确积分的方法求得位移及应力系数，积分在局部坐标中进行，再通过张量转换而获得。（3）对称条件的利

16、用本课题研究的薄板是几何形状及加载条件对称的弹性体，可用对称性大大减小的系数矩阵的规模。由于对称的位移及面力也是对称的，可以证明对称单元虚拟分布强度也对称，因此在满足途中对称条件的情况下，合并绝对值相同的未知数，方程组系数矩阵阶数 图4-2 对称单元可减少为原来的1/4,即只需对1/4边界划分单元。（4）求区域内各点位移及应力在解方程求出后，代入以下公式可以求得区域内任意点的应力及位移。应力： （4.9）其中 位移： （4.10）一般情况下，用三个约束条件可以处理刚体位移，在利用了对称条件后，相应的刚体位移可自动消除。（5） 边界单元间接法计算程序功能及程序流程图输入数据j=1是否为应力边界条

17、件计算应力系数计算位移系数送入系数矩阵j=j+1J n ?解方程求虚拟力强度求边界点及内点应力位移计算主应力等输入结果图4-3 边界元间接法程序流程图4.3 孔口应力集中问题 表4-1 圆孔孔口应力值单元数71631max2.99912.99993.0000计算圆孔内点(0,y)应力值的误差分布，可以看出内点应力值计算结果精度很高，但其误差随y的变化而变化，在y-R=0.010.1附近，误差达到最大，约为0.5%，两侧误差以较快速度减小。图4-4 圆孔内点应力值误差分布采用边界元间接法分析孔口应力集中问题只需对边界划分单元，十分简便，计算结果精度令人满意，线性代数方程组系数矩阵阶数比较低，计算

18、范围广，是一种十分有效的数值计算方法。5 总结小孔口问题的应力集中现象具有如下特征：一是有集中性，孔附近的应力远大于较远处的应力，且最大和最小的应力一般发生在孔边上。二是局部性，由于开孔引起的应力扰动，主要发生在距孔边1.5倍孔口尺寸(如圆孔直径)的范围内;在此区域外，由于开孔引起的应力扰动值一般小于5%，在工程上可以忽略不计。孔口的存在破坏了金属材料的连续性，使得主应力弯曲绕行，应力作用线在孔口附近密集弯曲，造成局部应力增大。孔口附近的应力远大于无孔时的应力，也远大于距孔口较远处的应力，由圣维南原理可知，在远离小孔的地方，孔口局部应力集中的影响将消失。应力集中是由于开口后的应力扰动引起的, 最大与最小应力一般发生在孔边上，且应力扰动主要发生在距孔边1.5倍孔口尺寸的范围内。孔口应力集中程度与孔口形状有关，圆形孔口的应力集中程度较低，应用较为广泛。边界元法只需对边界部分做划分，相比有限元分析方法更为简便易行，解题速度较快。而且在相同离散精度的条件下，边界元解的精度要高于有限元。但是，由于边界元方法形成的线性方程组的系数矩阵是满阵，所以在处理大规模问题时遇到了困难，解题的规模受到限制，仅适合于处理中小规模问题。利用有限元适合于解决大规模问题而边界元适合于解决无限域问题和解的精度高的特点，考虑将二者耦合起来，充分发挥两者的优势，以更好地解决实际问题。