基于python语言分析分析圆板中心开孔的应力集中问题

Using Python/C/C++/Fortran, write a finite element code to evaluate the response of plane system.

1) Your program should read a file named input that contains any useful information.

2) Your program should produce a file named output that contains the following results:

Using your completed code, solve the plane stress problem shown in Figure 1. A uniform thin plate with a central hole of radius a is subjected to uniaxial stress . Use the finite element method to determine the stress concentration factor given the physical data psi, a = 0.5 in., h = 3 in., w = 6 in., E = 2(10)⁶ psi, and Poisson’s ratio = 0.3. In addition, compare the approximate solution of finite element method with the exact solution based on the elasticity theory, and discuss the difference between them.

Figure 1

3) Your codes should be easy enough to follow. Here are some ways to improve the readability of your codes:

ØPut comments as much as you can!

ØPerform different parts of your finite element analysis in different modules and then call the modules in the main body of your code. For example, write separate modules for reading the input file, initialization of equations, assembly of different vectors and matrices, solving the equation systems, and writing the output file.

4) Email me (jingyh@hit.edu.cn) your code named Student_ID-Name-Project2.zip for testing. Also, prepare a short report on what different modules of your code do. Attach to your report the input file for the example system as well as the output file you will obtain from running the code for that example. Hand in the hardcopy of your report in class.

一、问题描述

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

二、理论分析

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

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

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

对于三个节点，有

还可以得到

将用表示：

代入，可以得到：

插值函数有如下形式：

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

三角形单元应变能为

其中：

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

单元总势能为

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

代入上式得到

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

. ，得到

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

四、程序与分析

#%%          

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))          

    def needs_refinement(vertices, area):          
        bary = np.sum(np.array(vertices), axis=0) / 3          
        max_area = 0.1 + (la.norm(bary, np.inf) - 1) * 0.1          
        return bool(area > max_area)          

    info = triangle.MeshInfo()          
    info.set_points(points)          
    info.set_holes([(0, 0)])          
    info.set_facets(facets)          

    mesh = triangle.build(info, refinement_func=needs_refinement)          

    mesh_points = np.array(mesh.points)+[6,3]          
    mesh_tris = np.array(mesh.elements)          
    print(mesh_points.shape)          
    print(mesh_tris.shape)          
    # print(mesh_points)          
    # print(mesh_tris)          
    pt.figure(figsize=(12,6))          
    pt.triplot(mesh_points[:, 0], mesh_points[:, 1], mesh_tris)          
    for i in range(mesh_points.shape[0]):     #对节点进行编号          
        pt.text (mesh_points[i][0], mesh_points[i][1],str(i), fontsize=6, color = "r", style = "italic")          


    for j in range(len(mesh_tris)):      #对单元进行编号          
        x1=mesh_points[mesh_tris[j][0]][0]          
        x2=mesh_points[mesh_tris[j][1]][0]          
        x3=mesh_points[mesh_tris[j][2]][0]          
        y1=mesh_points[mesh_tris[j][0]][1]          
        y2=mesh_points[mesh_tris[j][1]][1]          
        y3=mesh_points[mesh_tris[j][2]][1]          

        JUZHEN.append([x1, x2, x3, y1, y2, y3, mesh_tris[j][0], mesh_tris[j][1], mesh_tris[j][2]]) #三个节点横坐标，三个节点的纵坐标，节点编号，          
        pt.text ((x1+x2+x3)/3,(y1+y2+y3)/3,str(j), fontsize=6, color = "y", style = "italic")          
    return mesh_points, mesh_tris, JUZHEN          
    # print('\n', mesh_tris)          
    # print('\n', mesh_tris)          
mesh_points, mesh_tris, JUZHEN = main(0.5)          
# mesh_points1, mesh_tris1 = main(1)          

#%%          

print(np.shape(JUZHEN))          
JUZHEN          

#%%          

Z=pd.read_excel('./input.xlsx', sheet_name='Sheet1', header=None) #input原始数据          
Z = np.array(Z)  #将格式转换为array          
print(Z)          
R=Z[0,0]   #导入各项参数          
R          
xc=Z[1,0]          
yc=Z[1,1]          
xmin=Z[2,0]          
xmax=Z[3,0]          
ymin=Z[4,0]          
ymax=Z[5,0]          
m=0.1          
n=0.3          
E=Z[0,2]          
v=Z[1,2]          
#%%          
# for i in range(len(left_bounary)):          
#     print("nall")          
# print(left_bounary)          
# print(right_bounary)          
A1 = []          
for j in range(np.shape(JUZHEN)[0]):          
    x1 = JUZHEN[j][0]          
    x2 = JUZHEN[j][1]          
    x3 = JUZHEN[j][2]          
    y1 = JUZHEN[j][3]          
    y2 = JUZHEN[j][4]          
    y3 = JUZHEN[j][5]          
    n_array = np.array([[x1,y1,1],[x2,y2,1],[x3,y3,1]])          
    det = np.linalg.det(n_array) #计算方阵行列式的值          
    A=abs(1/2*det)   #三角形单元面积          
    JUZHEN[j].append(A)          
    A1.append(A)          
A = A1          
A          
# JUZHEN.append(A)          
print(A)          
#%%          
JUZHEN[:][:] #输出矩阵的全部行和列          
#%%          
x1 = []          
y1 = []          
x2 = []          
y2 = []          
x3 = []          
y3 = []          
for i in range(np.array(JUZHEN).shape[0]):          
    x1.append(JUZHEN[i][0])          
    x2.append(JUZHEN[i][1])          
    x3.append(JUZHEN[i][2])          
    y1.append(JUZHEN[i][3])          
    y2.append(JUZHEN[i][4])          
    y3.append(JUZHEN[i][5])          
x1 = np.array(x1) #将矩阵格式转成array格式进行运算          
x2 = np.array(x2)          
x3 = np.array(x3)          
y1 = np.array(y1)          
y2 = np.array(y2)          
y3 = np.array(y3)          
alpha1=x2*y3-x3*y2          
alpha2=x3*y1-x1*y3          
alpha3=x1*y2-x2*y1          
beta1=y2-y3          
beta2=y3-y1          
beta3=y1-y2          
gamma1=x3-x2          
gamma2=x1-x3          
gamma3=x2-x1          
n=321  #节点个数          
m=554  #单元个数          
F = np.zeros([2*n, 1])          
mesh_points          
#%%          
from scipy import integrate          
F = np.zeros([2*n, 1])          
for i in range(len(x1)):          
    if (x1[i]==12 and x2[i]==12):          
        # print(i)          
        # print(x1[i], y1[i], x2[i], y2[i])          
        # # 1节点位置          
        dian1 = JUZHEN[i][6]          
        # print(JUZHEN[i][6])          
        # # 2节点位置          
        dian2 = JUZHEN[i][7]          
        # print(JUZHEN[i][7])          
        def function(y):          
            return 1/(2*A[i])*(alpha1[i]+beta1[i]*12+gamma1[i]*y)          
        def function1(y):          
            return 1/(2*A[i])*(alpha2[i]+beta2[i]*12+gamma2[i]*y)          
        F[2*dian1,0]=F[2*dian1,0]+integrate.quad(function, y1[i],y2[i])[0]*1000          
        F[2*dian1+1,0]=0          
        F[2*dian2,0]=F[2*dian2,0]+integrate.quad(function1, y1[i],y2[i])[0]*1000          
        F[2*dian2+1,0]=0          
    elif (x1[i]==12 and x3[i]==12):          
        dian1 = JUZHEN[i][6]          
        dian3 = JUZHEN[i][8]          
        def function(y):          
            return 1/(2*A[i])*(alpha1[i]+beta1[i]*12+gamma1[i]*y)          
        def function2(y):          
            return 1/(2*A[i])*(alpha3[i]+beta3[i]*12+gamma3[i]*y)          
        F[2*dian1,0]=F[2*dian1,0]+integrate.quad(function, y1[i],y3[i])[0]*1000          
        F[2*dian3,0]=F[2*dian3,0]+integrate.quad(function2, y1[i],y3[i])[0]*1000          
        # print(i)          
    #     print(x1[i], y1[i], x3[i], y3[i])          
    elif (x2[i]==12 and x3[i]==12):          
        dian2 = JUZHEN[i][7]          
        dian3 = JUZHEN[i][8]          
        def function1(y):          
            return 1/(2*A[i])*(alpha1[i]+beta1[i]*12+gamma1[i]*y)          
        def function2(y):          
            return 1/(2*A[i])*(alpha3[i]+beta3[i]*12+gamma3[i]*y)          
        F[2*dian2,0]=F[2*dian2,0]+integrate.quad(function1, y2[i],y3[i])[0]*1000          
        F[2*dian3,0]=F[2*dian3,0]+integrate.quad(function2, y2[i],y3[i])[0]*1000          
    #     print(i)          
    #     print(x2[i], y2[i], x3[i], y3[i])          
    elif (x1[i]==0 and x2[i]==0):          
        print(i)          
        print(x2[i], y2[i], x3[i], y3[i])          
print(F)          
#%%          
# n_array = np.array([[x1,y1,1],[x2,y2,1],[x3,y3,1]])          
# det = np.linalg.det(n_array) #计算方阵行列式的值          
# A=abs(1/2*det)   #三角形单元面积          
# for i in range(np.shape(JUZHEN)[0]):          
A=1          
T=1          
Ve=A*T          
# print(Ve)          
D1=np.array([[1,v,0],[v,1,0],[0,0,(1-v)/2]])          
det = np.linalg.det(D1)          
# print(det)          
# print(D1)          
D=[[0,0,0],[0,0,0],[0,0,0]]          
for i in range(D1.shape[0]):          
    for j in range(D1.shape[1]):          
        if D1[i][j]==0:          
            D[i][j]=0          
        else:          
             D[i][j]=E/((1-v**2)*D1[i][j])          
# print(D)          
# where_are_nan = np.isnan(D)          
# where_are_inf = np.isinf(D)          
# D[where_are_nan] = 0          
# D[where_are_inf] = 0          
# print(D)          

# print(np.array([list(beta1), [0 for i in range(len(beta1))]]))          
# B1 = np.array([[list(beta1), list(beta2), list(beta3), [0 for i in range(len(beta1))],  [0 for i in range(len(beta1))],          
#                 [0 for i in range(len(beta1))]],          
#               [list(gamma1),list(gamma2),list(gamma3),[0 for i in range(len(beta1))],  [0 for i in range(len(beta1))],          
#                [0 for i in range(len(beta1))]],          
#               [list(gamma1), list(gamma2), list(gamma3),list(beta1), list(beta2),          
#                 list(beta3)]])          
# B1.shape          
# B1          
# B1=np.array([beta1,0,beta2,0,beta3,0,],[0,gamma1,0,gamma2,0,gamma3],[gamma1,beta1,gamma2,beta2,gamma3,beta3]])          
# print(B1)          
# det = np.linalg.det(B1)          
# B=1/(2*A)*B1          
# print(B)          
#%%          
# 建立总纲          
K = np.zeros([2*n, 2*n])          
for j in range(np.array(JUZHEN).shape[0]):          
    cc=[2*JUZHEN[j][6]-1,2*JUZHEN[j][6],2*JUZHEN[j][7]-1,2*JUZHEN[j][7],2*JUZHEN[j][8]-1,2*JUZHEN[j][8]]          
    B1=np.array([[beta1[j],0,beta2[0],0,beta3[j],0],[0,gamma1[j],0,gamma2[j],0,gamma3[j]],[gamma1[j],beta1[j],gamma2[j],beta2[j],gamma3[j],beta3[j]]])          
    B2=B1.T          
    B2=np.array(B2)          
    kk=np.dot(B2,D)          
    kk1=np.dot(kk,B1)          
    for II in range(6):          
        for JJ in range(6):          
            K[cc[II],cc[JJ]]=K[cc[II],cc[JJ]]+kk1[II,JJ]          
K          
# xlswrite('output2.xlsx',K)          
# import xlwt          
# import xlrd          
# # ws = wb.add_sheet('test')          
#          
#          
# K1 = np.linalg.inv(K)          
#%%          
%求每个单元的应力应变          
X=np.zero(mesh_points.shape[0],3)          
Y=np.zero(mesh_points.shape[0],3)          
UU=np.zero(mesh_points.shape[0],3)          
Epsilon=np.zero(mesh_points.shape[0],3)          
Epsilone=np.zero(mesh_points.shape[0],1)          
Sigma=np.zero(mesh_points.shape[0],3)          
Sigma_e=np.zero(mesh_points.shape[0],1)          
#%%          
# for j=1:N(1)          
# l=t(j,1);  m=t(j,2);  n=t(j,3);%节点全局标号          
#     x1=p(t(j,1),1);  x2=p(t(j,2),1);  x3=p(t(j,3),1);%节点横坐标          
#     y1=p(t(j,1),2);  y2=p(t(j,2),2);  y3=p(t(j,3),2);%节点纵坐标          
#     A=abs(1/2*det([x1 y1 1;x2 y2 1;x3 y3 1]));%三角形单元面积          
#     alpha1=x2*y3-x3*y2;          
#     alpha2=x3*y1-x1*y3;          
#     alpha3=x1*y2-x2*y1;          
#     beta1=y2-y3;          
#     beta2=y3-y1;          
#     beta3=y1-y2;          
#     gamma1=x3-x2;          
#     gamma2=x1-x3;          
#     gamma3=x2-x1;\          
#          
#     Ve=A*T;          
#     D=E/(1-v^2)*[1,v,0;v,1,0;0,0,(1-v)/2];          
#     B=1/(2*A)*[beta1,0,beta2,0,beta3,0;0,gamma1,0,gamma2,0,gamma3;gamma1,beta1,gamma2,beta2,gamma3,beta3];          
for j in range(np.shape(JUZHEN)[0]):          
    x1 = JUZHEN[j][0]          
    x2 = JUZHEN[j][1]          
    x3 = JUZHEN[j][2]          
    y1 = JUZHEN[j][3]          
    y2 = JUZHEN[j][4]          
    y3 = JUZHEN[j][5]          
    n_array = np.array([[x1,y1,1],[x2,y2,1],[x3,y3,1]])          
    det = np.linalg.det(n_array) #计算方阵行列式的值          
    A=abs(1/2*det)   #三角形单元面积          
    JUZHEN[j].append(A)          
    A1.append(A)          
    x1 = []          
    y1 = []          
    x2 = []          
    y2 = []          
    x3 = []          
    y3 = []          
# for i in range(np.array(JUZHEN).shape[0]):          
    x1.append(JUZHEN[j][0])          

    x2.append(JUZHEN[j][1])          

    x3.append(JUZHEN[j][2])          

    y1.append(JUZHEN[j][3])          

    y2.append(JUZHEN[j][4])          

    y3.append(JUZHEN[j][5])          

    x1 = np.array(x1) #将矩阵格式转成array格式进行运算          
    x2 = np.array(x2)          
    x3 = np.array(x3)          
    y1 = np.array(y1)          
    y2 = np.array(y2)          
    y3 = np.array(y3)          
    alpha1=x2*y3-x3*y2          
    alpha2=x3*y1-x1*y3          
    alpha3=x1*y2-x2*y1          
    beta1=y2-y3          
    beta2=y3-y1          
    beta3=y1-y2          
    gamma1=x3-x2          
    gamma2=x1-x3          
    gamma3=x2-x1          

f3=figure          
title('应力云图','Fontsize',24)          
figure(f3)          
hold on          
axis equal          
axis off          
for j=1:N(1)          
    patch(X(j,:)',Y(j,:)',Sigma_e(j,1))          
    colormap(bone)          
end

五、网格划分

六、应力云图

七、对比分析

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

八、总结

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