计算力学|采用python进行有限元模拟

从abaqus输出的inp文件中读取节点和单元信息

import meshio

mesh = meshio.read('Job-3.inp')

coords = mesh.points###coords即为各个节点的坐标

Edof = mesh.cells_dict['triangle']#Edof为三角形单元的节点号

1.单元刚度矩阵

def element_stiffness(n1,coords,E,v,t):

    node1 = coords[n1[0]]

    node2 = coords[n1[1]]

    node3 = coords[n1[2]]

    xi = node1[0]

    yi = node1[1]

    xj = node2[0]

    yj = node2[1]

    xm = node3[0]

    ym = node3[1]

    A=abs((xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))*0.5)

    bi=yj-ym

    bj=ym-yi

    bm=yi-yj

    ci=-xi+xm

    cj=-xm+xi

    cm=-xi+xj

B=np.array([[bi,0,bj,0,bm,0],[0,ci,0,cj,0,cm],[ci,bi,cj,bj,cm,bm]])/(2*A)

    BT=B.transpose()

    D=np.array([[1,v,0],[v,1,0],[0,0,(1-v)/2]])*E/(1-v*v)

    Ke=BT.dot(D).dot(B)*t*A

    return Ke

2.总体刚度矩阵

def assemble_stiffness(coords,edof,E,v,t):

    NoN = np.size(coords, 0)

    NoE = np.size(edof, 0)

    PD = np.size(coords, 1)

    NPE = np.size(edof, 1)

    K=np.zeros([NoN*PD,NoN*PD])

    for i in range(len(edof)):

        n1 = edof[i, 0:(NPE)]

        Ke=element_stiffness(i,n1,coords, E,v,t)

        for r in range(NPE):

            for p in range(NPE):

                K[2*n1[r] , 2*n1[p] ] = K[2*n1[r] - 0, 2*n1[p]-0] + Ke[2*r, 2*p]

                K[2*n1[r] , 2*n1[p]+1] = K[2*n1[r] , 2*n1[p]+1] + Ke[2*r, 2*p+1]

                K[2*n1[r] +1, 2*n1[p] ] = K[2*n1[r]+1, 2*n1[p] ] + Ke[2*r+1, 2*p]

                K[2*n1[r] +1, 2*n1[p]+1] = K[2*n1[r]+1 , 2*n1[p]+1] + Ke[2*r+1, 2*p+1]

    return K

3.施加载荷

在薄板的左右两端施加均布载荷，首先根据单元节点横坐标判断在左右两端，判断在左右两端的单元个数，然后在单元上设置近似的均布载荷。

def F(coords,NoN,PD,b,c):

    G = np.zeros([NoN , PD])

    U = np.zeros([NoN , PD])

    # 施加左侧载荷

    m=0

    for i in range(NoN):

        if int(coords[i,:][0]*10**2) == -(b/2)*10**2+1.:

            m=m+1

    p=c/m*689.5

    k = 0

    for i in range(NoN):

        if int(coords[i][0]*10**2) == -(b/2)*10**2+1.:

            print(i)

            k = k + 1           

            G[i][0]= (-p)

    # 施加右侧载荷

    n=0

    for i in range(NoN):

        if int(coords[i,0:][0]*10**2) == (b/2)*10**2-1.:

            n = n + 1

    p = c/ n*689.5

    for i in range(NoN):

        if int(coords[i,0:][0]*10**2) == (b/2)*10**2-1.:

            G[i][0] = p

    return G

4.求解位移

U = np.zeros([2*NoN , PD-1])

T= np.linalg.inv(K)#T为K（总体刚度矩阵）逆矩阵

U=np.dot(T,GT)#U=[K]-1[G]

U=U.reshape(2*NoN,PD-1)#把位移转化为一维数组

5.计算应力应变

epsilonx = np.zeros(NoE)

epsilony = np.zeros(NoE)

epsilonxy = np.zeros(NoE)

A = np.zeros([3,3])

for i in range(NoE):

      dot_id = Edof[i]

      dot1 = np.array(coords[int(Edof[i, 0])])

      dot2 = np.array(coords[int(Edof[i, 1])])

      dot3 = np.array(coords[int(Edof[i, 2])])

      beta1 = dot2[1] - dot3[1]

      beta2 = dot3[1] - dot1[1]

      beta3 = dot1[1] - dot2[1]

      gamma1 = -dot2[0] + dot3[0]

      gamma2 = -dot3[0] + dot1[0]

      gamma3 = -dot1[0] + dot2[0]

      B=np.array([[1, dot1[0], dot1[1]], [1, dot2[0], dot2[1]],[1, dot3[0], dot3[1]]]) / 2

      A = np.linalg.det(B)

      epsilonx[i - 1] = (beta1 * U[dot_id[0]] + beta2 * U[dot_id[1]] + beta3 * U[dot_id[2]]) / (2 * A)

      epsilony[i - 1] = (gamma1 * U[NoN + dot_id[0]] +

            gamma2 * U[NoN + dot_id[1]] + gamma3 * U[NoN + dot_id[2]]) / (2 * A)

      epsilonxy[i - 1] = (beta1 * U[NoN + dot_id[0]] +

             beta2 * U[NoN + dot_id[1]] + beta3 * U[NoN + dot_id[2]]

             + gamma1 * U[dot_id[0]] + gamma2 * U[dot_id[1]] + gamma3 * U[dot_id[2]]) / (2 *A)

      sigmax = (E / (1 - v ** 2)) * (epsilonx + v * epsilony)

      sigmay = (E / (1 - v ** 2)) * (epsilony + v * epsilonx)

      sigmaxy = (E / (2 * (1 + v))) * epsilonxy

五、结果及分析

编程结果：

X方向位移（Abaqus）

y方向位移（Abaqus）

应力：

数值比较：最大x方向位移

数值比较：最大y方向位移

数值比较：最大x方向应力

分析：在用Abaqus求解时，将中间整个圆环边界位移设置为零，相当于减少了薄板的位移量，因此出现了编程位移数值大于Abaqus数值的问题。而在使用相同模型时，位移的大小决定了应变的大小，应力呈现与应变相同的变化趋势，所以编程的应力大于Abaqus的应力，但结果在同一个数量级，可以基本确定是以上原因。