计算科学 - 4节点四边形单元的刚度矩阵计算 - 吾爱随笔录

我正在为传热（标量场问题）编写有限元代码，并从简单的 4 节点四边形元素开始。我尝试在物理坐标系中计算电导（刚度）矩阵，并将答案与等参系统进行比较。

元素看起来像：

 |Y
 |
4 _____________3
 |             |
 |             |
 |             |
 |             |
 |             |
 |_____________|  ---->X
 1             2

（边1-2不是直的。）节点的坐标1,2,3,4是：

1: (0,0)    
2: (2,0.5)    
3: (2,1)    
4: (0,1)

然后我对双单元进行等参数变换(e,n)：

       |n
       |
4 _____________3
 |     |       |
 |     |       |
 |     |------ | ----->e
 |             |
 |             |
 |_____________|  
 1             2

因此，节点的坐标1,2,3,4变为：

1: (-1,-1)    
2: (1,-1)    
3: (1,1)    
4: (-1,1)

我创建了一个 python simpy 脚本来计算刚度矩阵。

import sympy as sp
import numpy as np
from sympy import init_printing

n,e = sp.symbols('n e')
# shape functions in isoparametric
N1_iso = (1/4)*(e-1)*(n-1)
N2_iso = (1/4)*(e+1)*(-n+1)
N3_iso = (1/4)*(e+1)*(n+1)
N4_iso = (1/4)*(-e+1)*(n+1)

#gradient of shape function matrix in isoparametric
GN_mat = sp.Matrix([ [sp.diff(N1_iso,e),sp.diff(N2_iso,e), sp.diff(N3_iso,e), sp.diff(N4_iso,e)],
         [sp.diff(N1_iso,n),sp.diff(N2_iso,n), sp.diff(N3_iso,n), sp.diff(N4_iso,n)]])

#point matrix
P_mat = sp.Matrix([[0.0,0.0],[2.0,0.5],[2.0,1.0],[0.0,1.0]])
#Jacobian matrix
J_mat = (1.0/4.0)*(GN_mat*P_mat)

#Jacobian matrix determinant and inverse 
J_det=J_mat.det()
J_inv=J_mat.inv()

#"B" matrix, assuming unit conductivity matrix 
B_mat=J_inv*GN_mat

#conductance (stiffness) matrix
expr1 = B_mat.T*B_mat*J_det

#Double integrating conductance (stiffness) matrix
expr2 = sp.integrate(expr1,(n,-1,1))
#Stiffness (conductance) matrix using isoparametric transformation
K_mat=sp.integrate( (expr2),(e,-1,1))
print('K_mat = \n {}'.format(K_mat ))

#Now in physical coordinate system
x,y = sp.symbols('x y')
#Area of the element = 2
A=1.5

#shape function
N1 = (1/A)* (x-2)*(y-1)
N2 = (1/A)* (x-0)*(1-y)
N3 = (1/A)* (x-0)*(y-0.5)
N4 = (1/A)* (2-x)*(y-0)

B_mat2 = sp.Matrix([[sp.diff(N1,x), sp.diff(N2,x), sp.diff(N3,x), sp.diff(N4,x)],
                    [sp.diff(N1,y), sp.diff(N2,y), sp.diff(N3,y), sp.diff(N4,y)]])

expr10 = sp.integrate(B_mat2.T*B_mat2,(y,0,1))
#Stiffness (conductance) matrix using physical coordinates
K_mat2 = sp.integrate(expr10,(x,0,2))
print('     K_mat2 = \n {}'.format(K_mat2))

如果我提出第 2 点：(2,0.0)，域变为矩形，我得到相同的答案。但是，如果不是这样，我会得到不同的答案K_mat和K_mat2。有什么建议吗？