一、问题描述

有半径为a中心孔的均匀薄板受到单轴压力，应力为1000MPa，中心孔半径a = 0.5 in., 薄板高2h，宽2w，h = 3 in., w = 6 in., 弹性模量E = 2(10)⁶ psi，泊松比v=0.3，解决平面应力问题，并将有限元的近似解与基于弹性力学理论的精确解进行对比。

二、理论分析

考虑这类中心开孔方板，受到单轴拉力，圆孔圆心和矩形中心重合，处于平面应力状态，使用有限元求解此结构的变形图。

首先对此结构进行单元剖分，确定单元按照有限元的分析流程，要先对此结构进行单元剖分，确定单元与结点编号、以及单元的自由度编号。因为这里是平面应力问题，所以可以采用常应变三角形单元进行网格划分，并且采用的是非结构化的网格。

在该结构中采用的是常应变三角形单元，在整体坐标系中单元刚度矩阵均用下式进行计算：

对于三个节点，有

还可以得到

代入可以得到：

其中A为三角形单元面积，给定三角形节点坐标，可以通过下式得到：

三角形单元应变能为

在三角形单元中，外力做功为

单元总势能为

单元处于平衡状态，总势能最小，有

代入上式得到

其中k即为单元刚度矩阵。

三、圆孔的孔边应力集中理论

 

五、网格划分

六、应力云图

七、对比分析

有限元解(数值解)，最终输出的应力极值为3096MPa；弹性力学书上的理论解为3100MPa，原因是有限元网格划分所存在的误差，导致计算结果存在一定的误差，但由于误差不超过数值的5%，证明有限元仿真结果的准确性。

八、总结

有限元分析的最大特点就是标准化和规范化，这种特点时使大规模分析和计算成为可能。实现有限元分析标准化和规范化的载体就是单元，通过构造具有代表性的单元装配成复杂的结构。在此例中构造三角形单元，先列出小刚度矩阵，再进行装配，最后带入边界条件和外力得出位移、应力、应变的解，并且画出云图。云图中体现了应力集中，即在小孔周围网格密集（即应力大），在远离小孔的地方应力网格稀疏，符合了弹性力学的理论。弹性力学中，孔边的边界条件是极坐标下的正应力与切应力均为零。而通过有限元方法求出的孔边应力大致符合弹性力学理论

附录代码

#%%

import meshpy.triangle as triangle

import numpy as np

import numpy.linalg as la

import pandas as pd

import matplotlib.pyplot as pt

pt.figure(figsize=(12,6))#画布的大小

# Quadratic element demo, by Aravind Alwan

def round_trip_connect(start, end):##划分网格

   return [(i, i + 1) for i in range(start, end)] + [(end, start)]

mesh_points=[]

def main(R):

   JUZHEN = []

   # R = 0.5

   points = [(6, 0), (6, 3), (-6, 3), (-6, -3), (6, -3), (6, 0)]

   facets = round_trip_connect(0, len(points) - 1)

   circ_start = len(points)

   points.extend((R * np.cos(angle), R * np.sin(angle)) for angle in np.linspace(0, 2 * np.pi, 50, endpoint=False))

   facets.extend(round_trip_connect(circ_start, len(points) - 1))