import numpy as np import matplotlib.pyplot as plt import matplotlib.tri as tri from collections import OrderedDict # 一、定义函数 # 函数1：计算形函数矩阵N def Q4_N(kesi, yita): N1 = (1-kesi)*(1-yita)/4 N2 = (1+kesi)*(1-yita)/4 N3 = (1+kesi)*(1+yita)/4 N4 = (1-kesi)*(1+yita)/4 N = np.array([ [N1, 0, N2, 0, N3, 0, N4, 0], [0, N1, 0, N2, 0, N3, 0, N4] ]) return N # 函数2： 计算形矩阵G，用于计算应变矩阵 def Q4_G(kesi, yita): G = np.array([ [yita/4 - 1/4, 0, 1/4 - yita/4, 0, yita/4 + 1/4, 0, - yita/4 - 1/4, 0], [kesi /4 - 1/4, 0, - kesi/4 - 1/4, 0, kesi/4 + 1/4, 0, 1/4 - kesi/4, 0], [0,yita/4 - 1/4, 0, 1/4 - yita/4, 0, yita/4 + 1/4, 0, - yita/4 - 1/4], [0, kesi /4 - 1/4, 0, - kesi/4 - 1/4, 0, kesi/4 + 1/4, 0, 1/4 - kesi/4] ]) return G # 函数3： 计算雅可比矩阵J def Q4_J(kesi, yita, xy): NN = np.array([ [-(1-yita)/4, (1-yita)/4, (1+yita)/4, -(1+yita)/4], [-(1-kesi)/4, -(1+kesi)/4, (1+kesi)/4, (1-kesi)/4] ]) J = NN @ xy return J # 二、输入信息 # 定义结构材料、截面属性 E = 2.000E+08 # 弹性模量 v = 0.3 # 泊松比 L = 2 # 长度 H = 0.5 # 高度 t = 0.2 # 宽度 # 网格划分 DoE = 0.25 # 网格尺寸 NEofH = int(H/DoE) # 高度方向单元数 NEofL = int(L/DoE) # 长度方向单元数 NodeCoord = np.zeros(((NEofH+1)*(NEofL+1), 2)) # 节点坐标矩阵初始化 EleNode = np.zeros([int(NEofH*NEofL), 4]) # 单元节点矩阵初始化 # 定义节点坐标 for i in range(1, int(NEofL)+2): for j in range(1, int(NEofH)+2): NodeCoord[((i-1)*(int(NEofH)+1)+j)-1, 0] = (L/NEofL)*(i-1) NodeCoord[((i-1)*(int(NEofH)+1)+j)-1, 1] = (H/NEofH)*(j-1) # 定义单元节点 for i in range(1, int(NEofL)+1): for j in range(1, int(NEofH)+1): EleNode[((i-1)*int(NEofH)+j)-1, 0] = (NEofH+1)*(i-1) + 1+(j-1) EleNode[((i-1)*int(NEofH)+j)-1, 1] = (NEofH+1)*i + 1+(j-1) EleNode[((i-1)*int(NEofH)+j)-1, 2] = (NEofH+1)*i + 2+(j-1) EleNode[((i-1)*int(NEofH)+j)-1, 3] = (NEofH+1)*(i-1) + 2+(j-1) # 单元和节点数量 numEle = EleNode.shape[0] # 单元数量 numNode = NodeCoord.shape[0] # 节点数量 # 节点坐标值向量 xx = NodeCoord[:, 0] yy = NodeCoord[:, 1] # 自由度数量 numDOF = 2 * numNode # 约束自由度 ConstrainedDOF = np.zeros([1, 2*(int(NEofH)+1)], dtype='int') for i in range(1, 2*(int(NEofH)+1)+1): ConstrainedDOF[0][i -1] =i # 定义整体位移向量 displacement = np.zeros((1, numDOF)) # 定义整体载荷向量 force = np.zeros((1, numDOF)) p = -100 force[0][int(numDOF)-int(NEofH)-1] = p # 积分点及积分点权重 kesi_yita = np.array([ [-1/np.sqrt(3), -1/np.sqrt(3)], [1/np.sqrt(3), -1/np.sqrt(3)], [1/np.sqrt(3), 1/np.sqrt(3)], [-1/np.sqrt(3), 1/np.sqrt(3)] ]) w = np.array([1, 1]) ww = np.array([w[0]*w[0], w[1]*w[0], w[1]*w[1], w[0]*w[1]]) # 计算分析 # 定义总体刚度矩阵 stiffness_sum = np.zeros((numDOF, numDOF)) # 定义弹性矩阵D D = (E/(1-v**2))*np.array([ [1, v, 0], [v, 1, 0], [0, 0, (1-v)/2] ]) # 遍历每个单元，计算单元刚度矩阵并组装到总体刚度矩阵中 for i in range(numEle): noindex = EleNode[i] xy = np.zeros([4,2]) for j in range(4): xy[j, 0] = xx[int(noindex[j])-1] xy[j, 1] = yy[int(noindex[j])-1] eleK = np.zeros((8, 8)) for m in range(4): kesi = kesi_yita[m][0] yita = kesi_yita[m][1] J = Q4_J(kesi, yita, xy) A = (1/np.linalg.det(J))*np.array([ [J[1,1], -1*J[0,1], 0, 0], [0, 0, -1*J[1,0], J[0,0]], [-1*J[0,1], J[0,0], J[1,1], -1*J[1,0]] ]) G = Q4_G(kesi, yita) B = A @ G eleK =eleK+ww[m]*B.T @ D @ B * t * np.linalg.det(J) # 根据单元节点编号，确定单元自由度 eleDof = np.array([ noindex[0]*2-2, noindex[0]*2-1, noindex[1]*2-2, noindex[1]*2-1, noindex[2]*2-2, noindex[2]*2-1, noindex[3]*2-2, noindex[3]*2-1 ]) # 将单元刚度矩阵组装到总体刚度矩阵中 for i in range(8): for j in range(8): stiffness_sum[int(eleDof[i]), int(eleDof[j])] += eleK[i][j] # 根据给定的约束条件，删除指定被约束自由度对应的行和列 stiffness_sum_sim = np.delete(stiffness_sum, ConstrainedDOF-1, 0) stiffness_sum_sim = np.delete(stiffness_sum_sim, ConstrainedDOF-1, 1) # 解方程，求位移 activeDof = np.delete(np.arange(numDOF) +1, ConstrainedDOF-1) displacement[0][activeDof -1] = np.array(np.asmatrix(stiffness_sum_sim).I * np.asmatrix(force[0][activeDof -1]).T).T # 按节点编号输出两个方向的节点位移分量，格式：节点号，位移1，位移2 Nodeindex = np.arange(1, numNode+1).reshape(numNode, 1) tempdisp =np.hstack((Nodeindex, displacement.reshape(numNode, 2))) # 计算节点应力 stress_Node = np.zeros((numEle * 4, 5)) kesi_yita_Node = np.array([ [-1, -1], [1, -1], [1, 1], [-1, 1] ]) for i in range(numEle): noindex = EleNode[i] xy_node = np.zeros((4,2)) d = np.zeros((8,1)) for j in range(4): xy_node[(j,0)] = xx[int(noindex[j])-1] xy_node[(j,1)] = yy[int(noindex[j])-1] d[int(j*2)] = displacement[0][int(noindex[j]-1)*2] d[int(j*2+1)] = displacement[0][int(noindex[j]-1)*2+1] for m in range(4): kesi = kesi_yita_Node[m][0] yita = kesi_yita_Node[m][1] J = Q4_J(kesi, yita, xy) A = (1/np.linalg.det(J))*np.array([ [J[1,1], -1*J[0,1], 0, 0], [0, 0, -1*J[1,0], J[0,0]], [-1*J[1,0], J[0,0], J[1,1], -1*J[0,1]] ]) G = Q4_G(kesi, yita) B = A @ G sima = D @ B @ d stress_Node[i*4+m] = np.asarray([i+1, noindex[m], sima[0][0], sima[1][0], sima[2][0]], dtype=object) # 计算支座反力 reation =np.zeros((ConstrainedDOF.size, 1)) for i in range(ConstrainedDOF.size): rownum = ConstrainedDOF[0][i] reation[i] = stiffness_sum[rownum -1,] @ displacement.T - force[0][rownum -1] # 四、计算用于绘图相关参数 # 计算个单元最大长度maxlen maxlen = -1 for i in range(numEle): noindex = EleNode[i] deltax = xx[int(noindex[1]-1)] - xx[int(noindex[0]-1)] deltay = yy[int(noindex[1]-1)] - yy[int(noindex[0]-1)] L = np.sqrt(deltax**2 + deltay**2) if L > maxlen: maxlen = L # 计算节点最大绝对位移 maxabsDisp = np.sqrt((displacement[0][0] * displacement[0][0])+(displacement[0][1] * displacement[0][1])) for i in range(numNode): indexU = 2 * (i+1) -2 indexV = 2 * (i+1) -1 dispU = displacement[0][indexU] dispV = displacement[0][indexV] absDisp = np.sqrt(dispU**2 + dispV**2) if absDisp > maxabsDisp: maxabsDisp = absDisp # 计算变形图变形放大系数 factor = 0.1 scalefactor = 0.2 if maxabsDisp > 1e-30: factor = scalefactor * maxlen / maxabsDisp # 五、绘图 fig = plt.figure(figsize=(15, 6)) plt.subplots_adjust(left=0.1, right = 1) # 绘制为变形的结构线条 for i in range(numEle): noindex = np.concatenate([ [EleNode[i][0]], [EleNode[i][1]], [EleNode[i][2]], [EleNode[i][3]], [EleNode[i][0]],]) x = xx[noindex[:].astype('int')-1] y = yy[noindex[:].astype('int')-1] line1 = plt.plot(x, y, label='Undeformed Shape', color='grey', linestyle='--') # 绘制变形后的结构 xx_Deformed = np.empty_like(xx) yy_Deformed = np.empty_like(yy) for i in range(numEle): #获取第i个单元的节点索引，闭合曲线需要重复第一个节点方向索引 noindex = np.concatenate([[EleNode[i][0]], [EleNode[i][1]], [EleNode[i][2]], [EleNode[i][3]], [EleNode[i][0]]]) indexU= 2*noindex[:].astype('int')-2 indexV= 2*noindex[:].astype('int')-1 dispU= displacement[0][indexU] dispV= displacement[0][indexV] xxDeformed = xx[noindex[:].astype('int')-1] + dispU*factor yyDeformed = yy[noindex[:].astype('int')-1] + dispV*factor xx_Deformed[noindex.astype('int')-1] = xxDeformed yy_Deformed[noindex.astype('int')-1] = yyDeformed middlex = 0.5*(xxDeformed[1]+xxDeformed[0]) middley = 0.5*(yyDeformed[1]+yyDeformed[0]) line2 = plt.plot(xxDeformed, yyDeformed, label='Deformed Shape', color='black', marker='o', mfc='w') #标记节点编号 for i in range(numNode): plt.text(xx_Deformed[i]+0.03, yy_Deformed[i], str(i+1), color='red', ha='center', va='center') #设置坐标格式，使得图形以适合的比例进行显示 max_xx=max(max(xx), max(xx_Deformed)) max_yy=max(max(yy), max(yy_Deformed)) min_xx=min(min(xx),min(xx_Deformed)) min_yy=min(min(yy),min(yy_Deformed)) minx= -0.2*(max_xx-min_xx)+min_xx maxx=0.2*(max_xx-min_xx)+max_xx miny= -0.2*(max_yy-min_yy)+min_yy maxy=0.2*(max_yy-min_yy)+max_yy plt.xlim(0.2*minx,0.89*maxx) plt.ylim(miny, maxy) #设置x坐标轴刻度 plt.xticks(np.arange(0, 2.2,0.2)) #设置x坐标轴标签 plt.xlabel('X') #设置y坐标轴标签 plt.ylabel('Y') #绘制应力云图 sima_x= np.zeros((numNode, 1)) for stress_Node_sl in stress_Node: sima_x[int(stress_Node_sl[1])-1]= stress_Node_sl[2]#将方向应力数组赋值给变量2 Z=sima_x #将x坐标数组展平为一维数组#将y坐标数组展平为一维数组 x=xx_Deformed.flatten() y=yy_Deformed.flatten() #将应力数组展平为一维数组 z=Z.flatten() #构建三角剖分对象 triang =tri.Triangulation(x, y) tt=plt.tricontourf(triang,z,cmap="viridis") plt.colorbar(tt) #去掉重复标签 #绘制应力云图的颜色条 handles, labels = plt.gca().get_legend_handles_labels()#获取图例句柄和标签 by_label =OrderedDict(zip(labels,handles))#按照标签和句柄创建有序字典 plt.legend(by_label.values(),by_label.keys(), loc='upper right')#绘制图例，并设置位置为右上角 plt.show()
Q4单元介绍平面应力问题及平面应变问题的单元列式基本一直，差别仅在弹性矩阵 $[D]$ 。基本列式基本方程几何方程对于平面应力问题，正应力 $\sigma_z=0$，且剪应力 $\tau_{zx}=\tau_{xz}=0$，$\tau_{zy}=\tau_{yz}=0$，仅考虑单元平面内的三个平面应力分量 $\sigma_x$、$\sigma_y$ 和 $\tau_{xy}=\tau_{yx}$，根据单元仅有 $\varepsilon_x$、$\varepsilon_y$ 和 $\gamma_{xy}=\gamma_{yx}$ 三个应变分量，根据弹性力学的基本公式可得平面应力问题的几何方程的矩阵形式：$$ \left \{\varepsilon \right \} = \begin{Bmatrix}\varepsilon_x \\ \varepsilon_y \\ \varepsilon_{xy}\end{Bmatrix} = \begin{Bmatrix}\frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \end{Bmatrix} $$平面应变问题的几何方程与此相同。物理方程平面应力问题的应力-应变关系的矩阵形式：$$ \begin{Bmatrix}\sigma_x \\ \sigma_y \\ \tau_{xy}\end{Bmatrix}= = \begin{Bmatrix}\sigma \end{Bmatrix}=\left [ D \right ]\left \{ \varepsilon \right \}=\left [ D \right ]\begin{Bmatrix}\varepsilon_x \\ \varepsilon_y \\ \gamma_{xy}\end{Bmatrix}=\frac{E}{1-\mu^2}\begin{bmatrix}1 & \mu & 0 \\ \mu & 1 & 0 \\ 0 & 0 & \frac{1-\mu}{2}\end{bmatrix}\begin{Bmatrix}\varepsilon_x \\ \varepsilon_y \\ \gamma_{xy}\end{Bmatrix} $$其中 $E$ 为弹性模量， $\mu$ 为泊松比。位移场4 节点矩形单元有 4 个节点和 4 个边，每个节点有两个平动自由度，令单元节点位移向量 $\left \{ a \right \} = \begin{Bmatrix} u_1 & v_1 & u_2 & v_2 & u_3 & v_3 & u_4 & v_4\end{Bmatrix}^T$。单元中任意一点的节点位移 $(u,v)$ 可通过线性拉格朗日插值得到：$$ \begin{cases}u =N_1u_1+N_2u_2+N_3u_3+N_4u_4 \\ v =N_1v_1+N_2v_2+N_3v_3+N_4v_4\end{cases} $$用矩阵表示为：$$ \begin{Bmatrix}u \\ v\end{Bmatrix}=\begin{bmatrix}N_1 & 0 & \cdots & N_4 & 0 \\ 0 & N_1 & \cdots & 0 & N_4\end{bmatrix}\begin{Bmatrix}u_1 \\ v_1 \\ \cdots \\ u_4 \\ v_4\end{Bmatrix}=[N]\left \{ a \right \} $$将单元节点进行编号，并用自然坐标表示，此时单元形函数 $N_i$ 为：$$ \begin{cases}N_1(\xi,\eta) = (1-\xi)(1-\eta)/4 \\ N_2(\xi,\eta) = (1-\xi)(1-\eta)/4 \\ N_3(\xi,\eta) = (1-\xi)(1-\eta)/4 \\ N_4(\xi,\eta) = (1-\xi)(1-\eta)/4\end{cases} $$其中， $\xi\in[-1,1]$，$\eta\in[-1,1]$。可以看出在单元节点处，$N_i=0$。雅可比(Jacobian)矩阵假设一般函数 $f=f(x,y)$，它是 $\xi$ 和 $\eta$ 的隐式函数（变量 $x$ 和 $y$ 均是 $\xi$ 和 $\eta$ 的函数），根据链导法则可得：$$ \begin{cases} \frac{\partial f}{\partial \xi} = \frac{\partial f}{\partial x}\frac{\partial x}{\partial \xi} + \frac{\partial f}{\partial y}\frac{\partial y}{\partial \xi} \\ \frac{\partial f}{\partial \eta} = \frac{\partial f}{\partial x}\frac{\partial x}{\partial \eta} + \frac{\partial f}{\partial y}\frac{\partial y}{\partial \eta}\end{cases} $$写成矩阵形式：$$ \begin{Bmatrix}\frac{\partial f}{\partial \xi} \\ \frac{\partial f}{\partial \eta}\end{Bmatrix} = \begin{bmatrix}\frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta}\end{bmatrix} \begin{Bmatrix}\frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y}\end{Bmatrix} =[J]\begin{Bmatrix}\frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y}\end{Bmatrix} $$其中 $[J]$ 为雅可比矩阵，对于 $Q4$ 单元，可得雅可比矩阵为：$$ [J]=\begin{bmatrix}\frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta}\end{bmatrix}=\begin{bmatrix}\displaystyle \sum_{i=0}^4\frac{\partial N_i}{\partial \xi} x_i & \displaystyle \sum_{i=0}^4\frac{\partial N_i}{\partial \xi} y_i \\ \displaystyle \sum_{i=0}^4\frac{\partial N_i}{\partial \eta} x_i & \displaystyle \sum_{i=0}^4\frac{\partial N_i}{\partial \eta} y_i \end{bmatrix}=\begin{bmatrix}J_{11} & J_{12} \\ J_{21} & J_{22}\end{bmatrix} $$则有：$$ [J]^{-1}=\frac{1}{det([J])}\begin{bmatrix}J_{11} & -J_{12} \\ -J_{21} & J_{22}\end{bmatrix} $$几何方程与应变矩阵根据几何方程，将单元内任意一点的应变 $\left \{ \varepsilon \right \}$ 表示为单元节点位移 $\left \{ a \right \}$ 的函数，即：$$ \left \{ \varepsilon \right \}=\left [ B \right ]\left \{ a \right \} $$式中， $\left\{\varepsilon \right\}=\begin{Bmatrix}\varepsilon_x & \varepsilon_y & \tau_{xy}\end{Bmatrix}^T$，矩阵 $[B]$ 为应变矩阵。上式还可以写成：$$ \begin{Bmatrix}\varepsilon\end{Bmatrix}= \begin{Bmatrix}\varepsilon_x \\ \varepsilon_y \\ \gamma_{xy}\end{Bmatrix}=\begin{Bmatrix}\frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \end{Bmatrix}=\frac{1}{det([J])}\begin{bmatrix}J_{22} & -J_{12} & 0 & 0 \\ 0 & 0 & -J_{21} & J_{11} \\ -J_{21} & J_{11} & J_{22} & -J_{12}\end{bmatrix}\begin{Bmatrix} \frac{\partial u}{\partial \xi} \\ \frac{\partial u}{\partial \eta} \\ \frac{\partial v}{\partial \xi} \\ \frac{\partial v}{\partial \eta} \end{Bmatrix} $$进一步改写如下形式：$$ \begin{Bmatrix}\varepsilon\end{Bmatrix}= \frac{1}{det([J])}\begin{bmatrix}J_{22} & -J_{12} & 0 & 0 \\ 0 & 0 & -J_{21} & J_{11} \\ -J_{21} & J_{11} & J_{22} & -J_{12}\end{bmatrix}\begin{bmatrix}\frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} & 0 \\ \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} & 0 \\ 0 & \frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} \\ 0 & \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} \end{bmatrix}\begin{Bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \\ u_4 \\v_4 \end{Bmatrix} $$为便于计算应变矩阵 $[B]$，令 $[B]=[A][G]$，则有：$$ [A]=\frac{1}{det([J])}\begin{bmatrix}J_{22} & -J_{12} & 0 & 0 \\ 0 & 0 & -J_{21} & J_{11} \\ -J_{21} & J_{11} & J_{22} & -J_{12}\end{bmatrix} $$$$ [G]=\begin{bmatrix}\frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} & 0 \\ \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} & 0 \\ 0 & \frac{\partial N_1}{\partial \xi} & 0 & \frac{\partial N_2}{\partial \xi} & 0 & \frac{\partial N_3}{\partial \xi} & 0 & \frac{\partial N_4}{\partial \xi} \\ 0 & \frac{\partial N_1}{\partial \eta} & 0 & \frac{\partial N_2}{\partial \eta} & 0 & \frac{\partial N_3}{\partial \eta} & 0 & \frac{\partial N_4}{\partial \eta} \end{bmatrix} $$可见，对于 $Q4$ 单元，应变矩阵 $[B]$ 的元素是 $\xi$ 和 $\eta$ 的函数。物理方程与应力矩阵根据物理方程，单元的应力可表示为 $\{ \sigma \}=[D]\{ \varepsilon \}$。结合几何方程，可知任一点的应力可通过节点位移来表示 $\{ \sigma \} =[D][B]\{ a\}=[S]\{a\}$，其中 $[S]=[D][B]$ 为应力矩阵。对于平面应力玩问题，弹性矩阵 $\displaystyle [D]=\frac{E}{1-u^2}\begin{bmatrix} 1 & u & 0 \\\mu & 1 & 0 \\ 0 & 0 & \frac{1-\mu}{2} \end{bmatrix}$，若为平面应变问题，弹性矩阵 $\displaystyle [D]=\frac{E(1-\mu)}{(1+\mu)(1-2\mu)}\begin{bmatrix} 1 & \frac{\mu}{1-\mu} & 0 \\\frac{\mu}{1-\mu} & 1 & 0 \\ 0 & 0 & \frac{1-2\mu}{2(1-\mu)} \end{bmatrix}$，平面问题的弹性矩阵均为 $3\times 3$ 的对称矩阵。单元刚度矩阵单元的应变能可表示为：$$ U = \frac{1}{2}\int_V \left \{ \sigma \right \}^T \left \{ \varepsilon \right \} dV=\frac{1}{2}\left \{a \right \}^T \left[ \int_V [B]^T [D][B]dV \right ]\left \{ a \right \}=\frac{1}{2}\left \{ a \right \}^T\left [ t \int_{-1}^1 \int_{-1}^1 [B]^T[D][B]det([J])d\xi d\eta\right ]\left \{ a \right \} $$由最小势能原理可得单元的刚度矩阵：$$ [k]=t \int_{-1}^1 \int_{-1}^1 [B]^T[D][B]det([J])d\xi d\eta $$其中 $t$ 为单元厚度。对于 $Q4$ 单元，采用 $2\times2$ 个积分点的 Gauss 积分可求得精确结果，则刚度矩阵可表示为：$$ [k]=t\sum_{i=1}^2\sum_{j=1}^2 w_iw_j[B(\xi_i,\eta_j)]^T[D][B(\xi_i,\eta_j)]det[J(\xi_i,\eta_j)] $$其中 $w_i$、$w_j$ 为积分权重。单元载荷列阵及等效节点力单元的载荷主要包括面力（分布在单元边上的面力）、体力（分布在单位体积上的力）及温度作用等。以体力 $\{f_v\}$ 为例，给出其等效节点载荷 $\{f\}$ 的计算方法。单元的体力表示为 $\{ f_v\} = \{ f_{vx}, f_{vy}\}^T$。单元体力 $\{ f_v\}$ 在单元变形 $\{u\}$ 上做的外力功势能可表示为：$$ W=-\int_V\{u\}^T\{f_v\}dV $$由最小势能原理可得，Q4 单元体力的等效节点载荷列阵为：$$ \{f\}=\int_V[N]^T\{f_v\}dV=t\left[ \int_{-1}^1 \int_{-1}^1 [N]^Tdet([J])d\xi d\eta \right]\{f_v\} $$与刚度矩阵类似，上式计算可通过数值积分完成，即$$ \{ f \}=t\sum_{i=1}^2\sum_{j=1}^2w_iw_j\left [ N(\xi_i,\eta_j) \right ]^Tdet[J(\xi_i,\eta_j)]\{f_v\} $$
#桁架单元:桁架结构作用集中力荷载 import numpy as np import matplotlib.pyplot as plt #一、输入信息 #模尺材料参数 E=200000 #弹性模量 A=4532 #截面面积 EA=E*A #抗拉刚度 #节点坐标及单元节点序号 NodeCoord=np.array([[-4500,0],[-1500,0],[1500,0],[4500,0],[-3000,3000],[0,3000],[3000,3000]]) EleNode=np.array([[1,2],[2,3],[3,4],[5,6],[6,7],[1,5],[2,6],[3,7],[2,5],[3,6],[4,7]]) #获得单元数量和单元节点数量 numEle = EleNode.shape[0] # 11 numNode = NodeCoord.shape[0] # 7 #节点坐标值(x，y)列向量 xx= NodeCoord[:,0] yy=NodeCoord[:,1] #被约束自由度 ConstrainedDof = np.array([1,2,8]) #自由度总数量 numDOF = numNode*2 # 14 #整体位移向量定义 displacement = np.zeros((1, numDOF)) # 1*14 矩阵 #整体荷载向量定义 force = np.zeros((1, numDOF)) force[0][(10-1)]=-100000 # -1转换为索引 force[0][(12-1)]=-100000 force[0][(14-1)]=-100000 #二、计算分析 #定义总体刚度矩阵 #创建一个大小为 mumDOFx mumDOF的全零数组，用于存储整体刚度矩阵 stiffness_sum = np.zeros((numDOF,numDOF)) #遍历单元，求解单元刚度矩阵 for i in range(numEle): #求得单元刚度矩阵 noindex = EleNode[i] #获取第i个单元的节点索引 deltax = xx[noindex[1] - 1] - xx[noindex[0] - 1] #计算第i个单元在横向上的坐标差 deltay = yy[noindex[1] - 1] - yy[noindex[0] - 1] #计算第i个单元在纵向上的坐标差 L= np.sqrt(deltax *deltax + deltay * deltay) #计算第个单元的长度L C=deltax/L #计算角度的余弦值C S= deltay/L #计算角度的正弦值S eleK=(E*A/L)*np.array([ [C*C,C*S,-C*C,-C*S], [C*S,S*S,-C*S,-S*S], [-C*C,-C*S,C*C,C*S], [-C*S,-S*S,C*S,S*S] ]) #计算单元刚度矩阵 eleK #获得节点对应自由度 eleDof= np.array([noindex[0]*2-2,noindex[0]*2-1, noindex[1]*2-2,noindex[1]*2-1]) # 将单元刚度矩阵组装到整体刚度矩阵 for i in range(4): for j in range(4): stiffness_sum[int(eleDof[i]),int(eleDof[j])] = stiffness_sum[int(eleDof[i]),int(eleDof[j])] + eleK[i][j] #根据给定的约束条件，删除指定被约束自由度对应的行和列，将整体刚度矩阵进行简化处理 stiffness_sum_sim = np.delete(stiffness_sum, ConstrainedDof-1,0) stiffness_sum_sim = np.delete(stiffness_sum_sim, ConstrainedDof-1,1) #解方程,求位移 activeDof = np.delete((np.arange(numDOF)+1), ConstrainedDof-1) displacement[0][activeDof-1] = np.array(np.asmatrix(stiffness_sum_sim).I * np.asmatrix(force[0][activeDof-1]).T).T #计算单元轴力 axialForce =np.zeros((numEle,1)) for i in range(numEle): noindex = EleNode[i] deltax = xx[noindex[1]-1]-xx[noindex[0]-1] deltay = yy[noindex[1]-1]-yy[noindex[0]-1] L = np.sqrt(deltax * deltax + deltay *deltay) C = deltax/L S = deltay/L eleDof = np.array([noindex[0]*2-1, noindex[0]*2,noindex[1]*2-1,noindex[1]*2]) axialForce[i]=((EA/L)*np.asmatrix([-C,-S,C,S]))*np.asmatrix(displacement[0][eleDof - 1]).T #计算支座反力 reation = np.zeros((len(ConstrainedDof),1)) for i in range (len(ConstrainedDof)): reation[i]= np.asmatrix(stiffness_sum[(ConstrainedDof[i]-1)]) * np.asmatrix(displacement).T #三、计算用于绘图相关参数 #计算各单元最大长度maxlen maxlen =-1 for i in range(numEle): noindex= EleNode[i] deltax = xx[noindex[1]-1]-xx[noindex[0]-1] deltay = yy[noindex[1]-1]-yy[noindex[0]-1] L= np.sqrt(deltax*deltax + deltay * deltay) if L>maxlen: maxlen =L #计算节点最大绝对位移 maxabsDisp =0 for i in range(numNode): indexU=2*(i+1)-1 indexV= 2*(i+1) dispU= displacement[0][indexU-1] dispV= displacement[0][indexV-1] Disp =np.sqrt(dispU*dispU + dispV* dispV) if Disp> maxabsDisp: maxabsDisp = Disp #获取变形图变形放大系数 scalefactor =0.2 factor=1.1 if maxabsDisp> 1e-30: factor = scalefactor * maxlen/maxabsDisp #四、绘制图形 #遍历单元，绘制变形前结构形状 for i in range(numEle): noindex =EleNode[i] x=xx[noindex[:]-1] y=yy[noindex[:]-1] plt.plot(x, y,color='0.6', marker='o', mfc='w', linestyle='--') #遍历单元节点，标注单元节点编号 for i in range(numNode): plt.text(xx[i]+200,yy[i],i+1,color='r', ha='center', va='center') #定义变形后节点坐标值(x，y)列向量 xx_Deformed= np.empty_like(xx) yy_Deformed =np.empty_like(yy) # 日遍单元，绘制变形后结构形状 for i in range(numEle): noindex= EleNode[i] indexU= 2*noindex[:]-1 indexV=2*noindex[:] dispU= displacement[0][indexU-1] dispV= displacement[0][indexV-1] xxDeformed = xx[noindex[:]-1]+ dispU*factor yyDeformed = yy[noindex[:]-1] + dispV*factor xx_Deformed[noindex-1]= xxDeformed yy_Deformed[noindex-1]= yyDeformed plt.plot(xxDeformed, yyDeformed, color='black',marker='o', mfc='w',) #设置坐标格式，使得图形以适合的比例进行显示 max_xx = max(max(xx),max(xx_Deformed)) max_yy = max(max(yy),max(yy_Deformed)) min_xx = min(min(xx),min(xx_Deformed)) min_yy = min(min(yy), min(yy_Deformed)) minx =-0.1*(max_xx-min_xx)+min_xx maxx =0.1*(max_xx-min_xx)+max_xx miny = -0.1*(max_yy-min_yy)+min_yy maxy =0.1*(max_yy-min_yy)+max_yy plt.xlim(minx,maxx) plt.ylim(miny,maxy) plt.show() #绘制内力图 #遍历单元，绘制变形前结构形状 for i in range(numEle): noindex = EleNode[i] x=xx[noindex[:]-1] y=yy[noindex[:]-1] plt.plot(x, y, color='0.6', marker='o', mfc='w', linestyle='--') #定义变形后节点坐标值(x，y)列向量 xx_Deformed= np.empty_like(xx) yy_Deformed = np.empty_like(yy) #将内力值归一化，仅绘图用 absminForce = 1e10 axialForceStd = np.zeros((numEle, 1)) for i in range(numEle): if abs(axialForce[i]) <absminForce: absminForce = axialForce[i] for i in range(numEle): axialForceStd[i] = axialForce[i]/absminForce #遍历单元，绘制变形后结构形状 for i in range(numEle): noindex = EleNode[i] indexU= 2*noindex[:]-1 indexV= 2*noindex[:] dispU= displacement[0][indexU-1] dispV= displacement[0][indexV-1] xxDeformed = xx[noindex[:]-1] + dispU*factor yyDeformed = yy[noindex[:]-1] + dispV*factor xx_Deformed[noindex-1]= xxDeformed yy_Deformed[noindex-1]= yyDeformed if axialForce[i]>0: plt.plot(xxDeformed, yyDeformed, color='green', marker='o', mfc='w', linewidth = 2 * axialForceStd[i][0]) else: plt.plot(xxDeformed, yyDeformed, color='red', marker='o', mfc='w', linewidth = -2 * axialForceStd[i][0]) #遍历单元，标注单元编号 for i in range(numEle): noindex = EleNode[i] x=(xx[noindex[0]-1]+ xx[noindex[1]-1])/2 y=(yy[noindex[0]-1]+yy[noindex[1]-1])/2 plt.text(x, y, i+1, color='b', ha='center', va='center') #设置坐标格式，使得图形以适合的比例进行显示 max_xx = max(max(xx), max(xx_Deformed)) max_yy = max(max(yy), max(yy_Deformed)) min_xx = min(min(xx), min(xx_Deformed)) min_yy = min(min(yy), min(yy_Deformed)) minx= -0.1*(max_xx-min_xx)+min_xx maxx = 0.1*(max_xx-min_xx)+max_xx miny= -0.1*(max_yy-min_yy)+min_yy maxy = 0.1*(max_yy-min_yy)+max_yy plt.xlim(minx,maxx) plt.ylim(miny, maxy) plt.show()
桁架单元介绍桁架单元又称为杆单元，该单元具有以下特性：（1）单元仅承受轴向载荷，不承受弯矩；（2）单元节点仅有平动自由度，没有转动自由度，单元本身只发生轴向变形。桁架单元又可分为平面桁架和空间桁架，平面桁架单元每个节点又 2 个自由度，空间桁架单元每个节点有 3 个自由度。基本方程桁架单元仅在单元端点收到轴向力作用，对于单元中任意长度为 $dx$ 的隔离体，满足以下方程：几何方程：$$ \varepsilon =\frac{du}{dx} $$物理方程：$$ \sigma=E\varepsilon=E\frac{du}{dx} $$平衡方程：$$ F=\sigma A = EA\varepsilon=EA\frac{du}{dx} $$其中， $E$ 为弹性模量， $A$ 为截面面积， $F$ 为杆端轴力，$\varepsilon$ 为轴向应变，$\sigma$ 为轴向应力，$du$ 为轴向力 $F$ 作用下长度为 $dx$ 的隔离体变形。局部坐标系下单元刚度矩阵根据上述几何方程、物理方程和平衡方程，可得：$$ \left \{ f \right \}=\left [ k \right ]\left \{ a \right \} $$其中 $\left\{ k\right\}=\frac{EA}{L}\begin{bmatrix}1&-1 \\ -1 & 1\end{bmatrix}$，即局部坐标下桁架单元的刚度矩阵。单元的坐标转换矩阵整体坐标下单元的刚度矩阵为：$$ \left [ k \right ]=\frac{EA}{L}\begin{bmatrix}c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2\end{bmatrix} $$其中， $c=\cos\theta=\frac{x_2-x_1}{L}$，$s=\sin\theta=\frac{y_2-y_1}{L}$。单元的内力全局坐标下节点位移的单元应力：$$ F=EA\varepsilon=\frac{EA}{L}\begin{Bmatrix}-c & -s & c & s\end{Bmatrix}\begin{Bmatrix}u_1 \\ v_1 \\ u_2 \\ v_2\end{Bmatrix} $$
基本变量(1) 物体内一点的变形由它的三个位移分量来表示：$$ \begin{Bmatrix}u\end{Bmatrix}=\begin{bmatrix}u&v&w\end{bmatrix}^T $$(2) 体力体力也称体积力，就是分布在物体整个体积内部各个质点上的力，例如物体的重力、惯性力、电磁力等。一 般情况下，物体内各点的体力是不同的。物体所受的体积力可用列向量表示：$$ \begin{Bmatrix}f_v\end{Bmatrix}=\begin{bmatrix}f_{vx} & f_{vy} & f_{vz}\end{bmatrix}^T $$(3) 面力面力也称表面力，是分布在物体表面上的力，例如风力、静水压力、物体之间的接触力等。物体在其表面各点所受到的力，一般也是不同的。物体所受的表面力可用列向量表示：$$ \begin{Bmatrix}f_s\end{Bmatrix}=\begin{bmatrix}f_{sx} & f_{sy} & f_{sz}\end{bmatrix}^T $$(4) 外力边界面外法向量可用它的三个方向余弦来表示：$$ \begin{Bmatrix}n\end{Bmatrix}=\begin{bmatrix}n_{x} & n_{y} & n_{z}\end{bmatrix}^T $$(5) 集中力就是作用于物体某一点上的力。不论点是在物体内部还是在表面上，只在点附近很小区域存在作用力，其他区域数值为0，在有限元分析中集中力是常见形式。物体所受的点载荷可用它的三个分量表示：$$ \begin{Bmatrix}P\end{Bmatrix}=\begin{bmatrix}P_{x} & P_{y} & P_{z}\end{bmatrix}^T $$(6) 物体内一点的应力可用列向量表示：$$ \begin{Bmatrix}\sigma\end{Bmatrix}=\begin{bmatrix}\sigma_x & \sigma_y & \sigma_z & \tau_{yz} & \tau_{zx} & \tau_{xy}\end{bmatrix}^T $$正应力分量：每个面上有一个量，只要用一个下标，既可表示作用面，又能区分作用方向，用 $\sigma_x$、 $\sigma_y$、 $\sigma_z$ 分别表示 $x$ 面、$y$ 面、$z$ 面上的正力分量；剪应力分量：每个面上有两个量，共有6个量，要区分作用面及其作用方向，必须用两个下标表示。规定第1个下标表示作用面，第2个下标表示作用方向，如 $\tau_{xy}$ 表示 $x$ 面上与 $y$ 轴平行的剪应力分量。所有剪应力分量为 $\tau_{xy}$、 $\tau_{xz}$、 $\tau_{yx}$、 $\tau_{yz}$、 $\tau_{zx}$、 $\tau_{zy}$。六个剪应力分量之间存在一些互等关系，称为剪应力互等定理：$$ \tau_{xy}=\tau_{yx}\quad\tau_{xz}=\tau_{zx}\quad\tau_{yz}=\tau_{zy} $$说明作用在两个互相垂直面上并且垂直于两面交线的剪应力是互等的，二者大小相等，正负号相同，剪应力的实际方向同时指向或背离两面交线。因此，剪应力符号中的两个下标可以对调，不必再区分哪个下标表示作用面或是作用方向。(7) 物体内一点的应变可用列向量表示：$$ \begin{Bmatrix}\varepsilon\end{Bmatrix}=\begin{bmatrix}\varepsilon_x & \varepsilon_y & \varepsilon_z & \gamma_{yz} & \gamma_{zx} & \gamma_{xy}\end{bmatrix}^T $$所谓形变就是形状的改变，是描述物体受力后发生变形的相对力学量。物体的形变程度一般用六面单元体的三条互相垂直的棱边的伸缩程度以及棱边之间夹角的改变量来描述。棱边的相对伸长或缩短量称为正应变，有时也称线应变，用 $\varepsilon_x$、$\varepsilon_y$、$\varepsilon_z$ 表示，其意义为沿 3 个坐标轴方向的相对伸长量，正应变以伸长为正、缩短为负；棱边之间所夹直角的改变量称为剪应变，记为 $\gamma_{xy}$、$\gamma_{yz}$、$\gamma_{zx}$ ，如 $\gamma_{xy}$ 的意义为 $x$ 与 $y$ 轴夹角的改变量，以弧度表示，剪应变以直角变小为正、增大为负。平衡方程平衡方程是表示应力与外力（体力）之间关系的方程，3D弹性问题单元内任意一点的应力平衡方程表示为：$$ \frac{\partial\sigma_x}{\partial x}+\frac{\partial\tau_{yx}}{\partial y}+\frac{\partial\tau_{zx}}{\partial z}+f_{vx}=0 $$$$ \frac{\partial\tau_{xy}}{\partial x}+\frac{\partial\sigma_{y}}{\partial y}+\frac{\partial\tau_{zy}}{\partial z}+f_{vy}=0 $$$$ \frac{\partial\tau_{xz}}{\partial x}+\frac{\partial\tau_{yz}}{\partial y}+\frac{\partial\sigma_{z}}{\partial z}+f_{vz}=0 $$物理方程物理方程表示应变与应力之间的关系，即广义胡克定律。用应力表示应变的物理方程。$$ \varepsilon_x=\frac{1}{E}\left [ \sigma_x -\mu\left ( \sigma_y + \sigma_z \right ) \right],\quad \gamma_{xy}=\frac{\tau_{xy}}{G} $$$$ \varepsilon_y=\frac{1}{E}\left [ \sigma_y -\mu\left ( \sigma_x + \sigma_z \right ) \right],\quad \gamma_{yz}=\frac{\tau_{yz}}{G} $$$$ \varepsilon_z=\frac{1}{E}\left [ \sigma_z -\mu\left ( \sigma_x + \sigma_y \right ) \right],\quad \gamma_{zx}=\frac{\tau_{zx}}{G} $$式中，$E$、$G$、$\mu$ 分别为材料的弹性模量、剪切模量和泊松比。三者的关系为$$ G=\frac{E}{2(1+\mu)} $$三个方向的正应变与正应力之间存在耦合关系，其影响量为泊松比。剪应变与剪应力之间不存在耦合。物理方程也有用应变表示应力的形式。几何方程几何方程又称柯西方程，表示应变与位移之间的关系。在小变形情况下，应变-位移关系可以用以下公式描述$$ \varepsilon_x=\frac{\partial u}{\partial x},\quad \gamma_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} $$$$ \varepsilon_y=\frac{\partial v}{\partial y},\quad \gamma_{yz}=\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y} $$$$ \varepsilon_z=\frac{\partial w}{\partial z},\quad \gamma_{zx}=\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z} $$可见，剪应变是两个方向变形的协调性问题，如 $\tau_{xy}$ 就是 $x$ 与 $y$ 方向变形的相互影响。
边界条件表示在边界上位移与约束或应力与面力之间的关系。按照边界条件的不同，分为位移边界条件、应力边界条件、混合边界条件。位移边界条件弹性体在所有边界上的位移分量均为已知，在边界上存在 $$ u_s =\overline{u}, \quad v_s = \overline{v}, \quad w_s=\overline{w} $$其中，$u_s$、$v_s$、$w_s$ 表示边界上的位移分量，而 $\overline{u}$、$\overline{v}$、$\overline{w}$ 则表示在边界上是坐标的已知函数。;应力边界条件弹性体在全部边界上所受的面力都是已知的。面力已知的条件可以转换成应力方面的条件，表示为$$ \overline{X}=l\sigma_x+m\tau_{yx}+n\tau_{zx} $$$$ \overline{Y}=l\tau_{xy}+m\sigma_y+n\tau_{zy} $$$$ \overline{Z}=l\tau_{xz}+m\tau_{yz}+n\sigma_z $$式中，$l=\cos(N,x)$，$m=\cos(N,y)$，$n=\cos(N,z)$ 表示边界外法线 $N$ 的方向余弦。上式同样是一个函数方程，表示边界上每一点应力与面力之间的关系，即应力分量的边界值就等于对应的面力分量。但应注意，应力分量与面力分别作用于不同面上，且各有不同的正负号规定。当边界的外法线 $N$ 沿坐标轴正向时，应力与面力符号相同；当边界的外法线N沿坐标轴负向时，二者符号相反。混合边界条件弹性体在一部分位移边界 $S_u$ 上的位移分量为已知，满足位移边界条件式；而在另一部分应力边界 $S_{\tau}$ 上的应力分量为已知，满足应力边界条件式。混合边界条件是实际工程问题中常见的边界形式。弹性力学的基本方程共计15个，即平衡方程3个、几何方程6个、物理方程6个；基本未知量为应力分量6个、应变分量6个、位移分量3个，共计也是15个。基本未知量的数目恰好等于基本方程的数目，在适当的边界条件下，才有可能从基本方程求解这些未知量。说明一点，应用有限元法分析工程问题时，位移边界条件、混合边界条件均可获得很好的结果；而对完全应力边界条件将无法求解，因为应力边界条件只满足力平衡关系，而对刚体位移没有限制。
工程结构分析可以归结为求解弹性力学问题。弹性力学就是研究弹性体或结构由于外部载荷作用或者温度改变等因素所引起的形变、应力和位移。通过力学的平衡条件、变形几何协调性等方法，在外力与应力、形变、位移等各个物理量之间建立基本方程。假设对结构输入某种形式的已知量或者说某个量发生改变之后(通常是外力)，再根据弹性力学的基本方程，寻求采用适当的方法求解这些基本方程，以便得到所需要的相关未知量。弹性力学研究的外力、应力、应变、位移这四个基本量，通过三组基本方程，即应力与外力关系的平衡方程、应力与应变关系的物理方程、应变与位移的几何方程联系起来。应力、应变、位移属于弹性体内部的物理量，严格遵循弹性力学基本方程。一般弹性力学方法是根据弹性力学基本方程，对整个弹性体或结构来分析的。假设应力场能够确定，先根据物理方程求应变，再根据几何方程求位移，几何方程是偏微分方程组，由应变计算位移需要积分，求解起来非常困难，有时是不可能的。如果弹性体内部的位移场能够确定，那么根据几何方程由位移求导计算应变；物理方程是代数方程，则容易由应变确定应力。但问题是计算的应力分量很难满足平衡方程，要得到正确解答同样也非常困难。有限元法放弃了应力与外力之间的平衡方程。根据能量原理，直接在外力与位移场之间建立联系，也就是反映外力改变与弹性体内部位移场变化的关系，依次求解位移、应变、应力，其过程如图所示。现在关键的问题就是对于复杂形体的弹性体或结构中的位移场该如何表达?有限元的基本思想是，在弹性体内选取足够多、有限数量的具有代表性的点，假定这些点的位移已知，再用这些假定的位移量描述其他位置点的位移，就得到了用特定点位移表示的弹性体的位移场，这些选定的有代表性的点称为节点。为了较准确反映原结构特征，通常，尖点、 拐角、截面改变处、位移约束位置、集中载荷作用点等特殊、具有代表性的点都应选为节点。所谓节点数量足够多，就是要保证所得到的数值结果满足精度要求；节点数量有限，就是将连续体中无穷多点的无限自由度转变为有限数量的自由度，以便能够用计算机进行数值计算和处理。在用节点位移分量描述其他点(非节点)的位移分量时，并非考虑所有节点位移的影响，而是先将弹性体按节点分割成很多简单形状的小区域或小块，仅由包含该点的小区域上节点来表示。有限元法是将复杂的连续体分割成有限多个简单形状的小区域或小块，这些小区域或小块称为单元。每个单元形状简单，以位移为基本变量，根据虚位移原理或最小势能原理，求解单元的能量，将每个单元的能量叠加后便得到结构的总能量。有限元法先化整为零、再积零为整。也就是把一个连续体分割成有限数量个单元，首先对每个简单的单元进行分析，然后再将所有单元组合起来进行结构的整体分析。有限元法的实质是将连续体无限自由度问题转化为离散体有限自由度问题，将连续场函数的偏微分方程求解问题转化成有限个参数的代数方程组求解问题。从数学的角度来看，有限元法是将偏微分方程化成代数方程组，然后利用计算机进行求解的方法。
模态分析：是研究结构动力特性一种方法，一般应用在工程振动领域。其中，模态是指机械结构的固有振动特性，每一个模态都有特定的固有频率、阻尼比和模态振型。分析这些模态参数的过程称为模态分析。按计算方法，模态分析可分为计算模态分析和试验模态分析。静力学分析：用来分析结构在给定静力载荷作用下的响应。一般情况下，比较关注的往往是结构的位移、约束反力、应力以及应变等参数。模态向量有一个很重要的特性,即“模态正交性'”。所谓模态分析法，就是利用系统固有的模态的正交性，以系统的各阶模态向量所组成的模态矩阵作为交换矩阵，对通常选取的物理坐标进行线性变换，使得振动系统以物理坐标和物理参数所描述、互相耦合的运动方程组，能够变为一组彼此独立的方程，每个独立方程只含一个独立的模态坐标。这个用模态坐标和模态参数所描述的各个独立方程，称为模态方程。模态分析实质上是一种坐标变换，其目的是为了解除方程的耦合，便于求解由于坐标变换是线性变换，因而系统在原有的物理系统中，对于任意激励的响应，便可视为系统各阶模态的线性组合，故模态分析法，又称模态叠加法。而各阶模态在叠加中所占的比重或加权系数，则取决于各阶的模态坐标响应。一般来说,高阶模态比低阶模态的加权系数要小得多，通常只需要选取前阶模态进行叠加，即可达到足够的精度。由此可知模态分析的主要优点就在于，它能用较少的运动方程或自由度数，直观、简明而又相当精确的去反映一个比较复杂结构的动态特性，从而大大减小测量、分析及计算的工作量。
有限元分析里面的模态分析的阶数
在理解“阶”之前，要先理解与“阶”紧密相连的名词“自由度”。自由度是指用于确定结构空间运动位置所需要的最小、独立的坐标个数。空间上的质点有三个自由度，分别为三个方向的平动自由度；空间上的刚体有六个自由度，分别为三个平动、三个转动自由度。一个连续体实际上有无穷多个自由度，有限元分析时将连续的无穷多个自由度问题离散成为离散的有限多个自由度的问题，此时，结构的自由度也就有限了。因此，可以这样理解，一个自由度对应一阶，连续体有无穷多阶。像弹簧--质量模型为单自由度系统，故对应的频率只有一阶。两自由度系统有两阶。一个具体的系统，每一阶对应着特定的频率、阻尼和模态振型。因此，结构理论上有无穷多阶固有频率，但是我们通常只关注低阶频率，特别是第一阶，也称为基频，这是因为，频率越低，越容易被外界激励起来。模态分析的实质是计算结构振动特征方程的特征值和特征向量。从数学的角度发现模态的灵魂在数学中，结构的频率和振型问题实际就是描述结构的刚度矩阵和质量矩阵相乘得到的矩阵的特征值和特征向量。线性代数的术语来说，振型(mode shape)其实就是特征向量在动力学中的一个物理表现。而我们所提取模态的阶数即对应要获取的方程中特征值的个数。实际的分析对象是无限维的，所以其模态具有无穷阶。但是对于运动起主导作用的只是前面的几阶模态，所以计算时根据需要指定提取前几阶进行计算。复杂的振动一般都可分解为简单振动的组合，而且，这些个简单振动，跟外来的激励样式无关，只跟物体的本身的性质以及边界约束条件有关。求振型的过程，就是把复杂振动“提纯”（数学术语叫做解耦，decoupling）的过程。例如当简支梁受到不同形式的外力时，会有不同的振动样式，再复杂的形式，也不过是前几阶振型的线性组合。由于各阶振型在整个振动中所占的比例不同，在宏观上就表现为振动形态有所不同。找出了振型，就抓住了振动的本质特征，振型是特征向量的一种表现形式。物体按照某一阶固有频率振动时，物体上各个点偏离平衡位置的位移是满足一定的比例关系的，可以用一个向量表示，这个就称之为模态，即振动形态。一阶模态是外力的激励频率与物体固有频率相等的时候出现的，此时物体的振动形态叫做一阶振型或主振型；从工程应用的角度认识模态分析的精髓低阶模态的模态刚度相对比较弱，在同样量级的激励作用下，响应会相对所占的权值大一些，所以，工程上低阶模态比较受关照，理论上低阶模态理论也相对成熟。在实验中，我们就是通过用一定的频率对结构进行激振，观测相应点的位移状况，当观测点的位移达到最大时，此时频率即为固有频率。实际结构的振动形态并不是一个规则的形状，而是各阶振型相叠加的结果。对于没有约束的对象，前6阶为刚体位移模态，频率为0；而对于有约束的对象，则没有刚体模态。约束施加的正确与否，对结构模态分析的影响十分显著，因此对于该问题应十分注意，保证对模型施加的约束与实际情况尽量符合。所以模态分析的目的就是要得到结构的振型和固有频率。所得到的应力、应变、位移值都没有实际量化意义，只能用于定性地考察比较；模态分析的意义在于了解结构的共振区域，为结构设计提供指导，它是开展其它动力学特性分析的基础；为结构系统的振动特性、振动故障诊断以及结构动力特性的优化设计提供依据。总之，模态是一只道行很高的妖，值得一捉。