三维刚度矩阵K的详细推导

前置弹性力学矩阵知识

几何变换矩阵S

S=[x000y000zyx00zyz0x](1.1)\mathbf{S} = \begin{bmatrix} \frac{\partial}{\partial x} & 0 & 0 \\ 0 & \frac{\partial}{\partial y} & 0 \\ 0 & 0 & \frac{\partial}{\partial z} \\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \\ 0 & \frac{\partial}{\partial z} & \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} & 0 & \frac{\partial}{\partial x} \end{bmatrix}\tag{1.1}

位移矩阵u

u={uvw}(1.2)\mathbf{u} = \begin{Bmatrix} u \\ v \\ w \end{Bmatrix}\tag{1.2}

应变矩阵ϵ\epsilon

由(1.1)式和(1.2)式导出:

ε={εxεyεzγxyγyzγzx}={uxvywzuy+vxvz+wywx+uz}=Su(1.3)\boldsymbol{\varepsilon} = \left\{ \begin{array}{c} \varepsilon_x \\ \varepsilon_y \\ \varepsilon_z \\ \gamma_{xy} \\ \gamma_{yz} \\ \gamma_{zx} \end{array} \right\} = \left\{ \begin{array}{c} \dfrac{\partial u}{\partial x} \\ \dfrac{\partial v}{\partial y} \\ \dfrac{\partial w}{\partial z} \\ \dfrac{\partial u}{\partial y} + \dfrac{\partial v}{\partial x} \\ \dfrac{\partial v}{\partial z} + \dfrac{\partial w}{\partial y} \\ \dfrac{\partial w}{\partial x} + \dfrac{\partial u}{\partial z} \end{array} \right\} = \mathbf{S}\mathbf{u} \tag{1.3}

节点位移a

ae={aiajamap}ai={uiviwi}etc.(1.4)\begin{aligned} \mathbf{a}^e &= \left\{ \begin{array}{c} \mathbf{a}_i \\ \mathbf{a}_j \\ \mathbf{a}_m \\ \mathbf{a}_p \\ \end{array} \right\} \\ \mathbf{a}_i &= \left\{ \begin{array}{c} u_i \\ v_i \\ w_i \\ \end{array} \right\} \quad \text{etc.} \end{aligned}\tag{1.4}

形函数N

设任意点位移为

u=α1+α2x+α3y+α4zv=α5+α6x+α7y+α8zw=α9+α10x+α11y+α12z(1.5)\begin{aligned} u &= \alpha_1 + \alpha_2 x + \alpha_3 y + \alpha_4 z \\ v &= \alpha_5 + \alpha_6 x + \alpha_7 y + \alpha_8 z \\ w &= \alpha_9 + \alpha_{10} x + \alpha_{11} y + \alpha_{12} z \end{aligned}\tag{1.5}

每个节点的空间位置也已知:

Xi=(xi,yi,zi),etc(1.6)X_i=(x_i,y_i,z_i), etc\tag{1.6}

对于uu,可以得到

ui=α1+α2xi+α3yi+α4ziuj=α1+α2xj+α3yj+α4zjum=α1+α2xm+α3ym+α4zmup=α1+α2xp+α3yp+α4zp(1.7)\begin{aligned} u_i &= \alpha_1 + \alpha_2 x_i + \alpha_3 y_i + \alpha_4 z_i \\ u_j &= \alpha_1 + \alpha_2 x_j + \alpha_3 y_j + \alpha_4 z_j \\ u_m &= \alpha_1 + \alpha_2 x_m + \alpha_3 y_m + \alpha_4 z_m \\ u_p &= \alpha_1 + \alpha_2 x_p + \alpha_3 y_p + \alpha_4 z_p \end{aligned}\tag{1.7}

写成矩阵形式

[1xiyizi1xjyjzj1xmymzm1xpypzp]XR4×4[α1α2α3α4]=[uiujumup](1.8)\underbrace{ \begin{bmatrix} 1 & x_i & y_i & z_i \\ 1 & x_j & y_j & z_j \\ 1 & x_m & y_m & z_m \\ 1 & x_p & y_p & z_p \end{bmatrix}}_{\mathbf{X} \in \mathbb{R}^{4 \times 4}} \cdot \begin{bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \end{bmatrix} = \begin{bmatrix} u_i \\ u_j \\ u_m \\ u_p \end{bmatrix}\tag{1.8}

可以解出

u=16V[(ai+bix+ciy+diz)ui+(aj+bjx+cjy+djz)uj+(am+bmx+cmy+dmz)um+(ap+bpx+cpy+dpz)up]v=16V[(ai+bix+ciy+diz)vi+(aj+bjx+cjy+djz)vj+(am+bmx+cmy+dmz)vm+(ap+bpx+cpy+dpz)vp]w=16V[(ai+bix+ciy+diz)wi+(aj+bjx+cjy+djz)wj+(am+bmx+cmy+dmz)wm+(ap+bpx+cpy+dpz)wp](1.9)\begin{aligned} u &= \frac{1}{6V} \Big[ (a_i + b_i x + c_i y + d_i z) u_i + (a_j + b_j x + c_j y + d_j z) u_j \\ &\quad + (a_m + b_m x + c_m y + d_m z) u_m + (a_p + b_p x + c_p y + d_p z) u_p \Big] \\ \\ v &= \frac{1}{6V} \Big[ (a_i + b_i x + c_i y + d_i z) v_i + (a_j + b_j x + c_j y + d_j z) v_j \\ &\quad + (a_m + b_m x + c_m y + d_m z) v_m + (a_p + b_p x + c_p y + d_p z) v_p \Big] \\ \\ w &= \frac{1}{6V} \Big[ (a_i + b_i x + c_i y + d_i z) w_i + (a_j + b_j x + c_j y + d_j z) w_j \\ &\quad + (a_m + b_m x + c_m y + d_m z) w_m + (a_p + b_p x + c_p y + d_p z) w_p \Big] \end{aligned}\tag{1.9}

其中

6V=det1xiyizi1xjyjzj1xmymzm1xpypzp(1.10)6V = \det \begin{vmatrix} 1 & x_i & y_i & z_i \\ 1 & x_j & y_j & z_j \\ 1 & x_m & y_m & z_m \\ 1 & x_p & y_p & z_p \\ \end{vmatrix}\tag{1.10}

ai=detxjyjzjxmymzmxpypzpbi=det1yjzj1ymzm1ypzpci=detxj1zjxm1zmxp1zpdi=detxjyj1xmym1xpyp1(1.11)\begin{aligned} a_i &= \det \begin{vmatrix} x_j & y_j & z_j \\ x_m & y_m & z_m \\ x_p & y_p & z_p \\ \end{vmatrix} \\ b_i &= -\det \begin{vmatrix} 1 & y_j & z_j \\ 1 & y_m & z_m \\ 1 & y_p & z_p \\ \end{vmatrix} \\ c_i &= -\det \begin{vmatrix} x_j & 1 & z_j \\ x_m & 1 & z_m \\ x_p & 1 & z_p \\ \end{vmatrix} \\ d_i &= -\det \begin{vmatrix} x_j & y_j & 1 \\ x_m & y_m & 1 \\ x_p & y_p & 1 \\ \end{vmatrix} \end{aligned}\tag{1.11}

其他的系数由更换角标得到,按照右手准则按顺序变换。

可以得到形函数NiN_i:

Ni=16V[ai+bix+ciy+diz000ai+bix+ciy+diz000ai+bix+ciy+diz](1.12)N_i=\frac{1}{6V} \begin{bmatrix} a_i + b_i x + c_i y + d_i z & 0 & 0 \\ 0 & a_i + b_i x + c_i y + d_i z & 0 \\ 0 & 0 & a_i + b_i x + c_i y + d_i z \\ \end{bmatrix}\tag{1.12}

形函数插值表达式u=Naeu=Na^e

ε=SNae=Bae=[BiBjBmBp]ae(1.13)\boldsymbol{\varepsilon} = \mathbf{S} \mathbf{N} \mathbf{a}^e = \mathbf{B} \mathbf{a}^e = \begin{bmatrix} \mathbf{B}_i & \mathbf{B}_j & \mathbf{B}_m & \mathbf{B}_p \end{bmatrix} \mathbf{a}^e\tag{1.13}

应变-位移矩阵B

由(1.13)式导出,即B=SNB=SN

Bi=[Nix000Niy000NizNiyNix00NizNiyNiz0Nix]=16V[bi000ci000dicibi00dicidi0bi](1.14)\mathbf{B}_i = \begin{bmatrix} \frac{\partial N_i}{\partial x} & 0 & 0 \\ 0 & \frac{\partial N_i}{\partial y} & 0 \\ 0 & 0 & \frac{\partial N_i}{\partial z} \\ \frac{\partial N_i}{\partial y} & \frac{\partial N_i}{\partial x} & 0 \\ 0 & \frac{\partial N_i}{\partial z} & \frac{\partial N_i}{\partial y} \\ \frac{\partial N_i}{\partial z} & 0 & \frac{\partial N_i}{\partial x} \end{bmatrix} = \frac{1}{6V} \begin{bmatrix} b_i & 0 & 0 \\ 0 & c_i & 0 \\ 0 & 0 & d_i \\ c_i & b_i & 0 \\ 0 & d_i & c_i \\ d_i & 0 & b_i \end{bmatrix}\tag{1.14}

其他块直接更换角标即可。

各向同性弹性矩阵D

D=E(1+ν)(12ν)[1ννν000ν1νν000νν1ν00000012ν200000012ν200000012ν2](1.15)\mathbf{D} = \frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & \nu & 0 & 0 & 0 \\ \nu & 1-\nu & \nu & 0 & 0 & 0 \\ \nu & \nu & 1-\nu & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} \end{bmatrix} \tag{1.15}

刚度矩阵K表达式推导

基于虚功原理的推导

为使节点力在静力上等效于真实边界应力与体力,最简单的方法是对节点施加一个任意的虚位移 δae\delta \mathbf{a}^e,然后使该虚位移下的外功等于内功

1. 虚位移与应变表达式

δu=Nδae,δε=Bδae(2.1)\delta \mathbf{u} = \mathbf{N} \, \delta \mathbf{a}^e, \quad \delta \boldsymbol{\varepsilon} = \mathbf{B} \, \delta \mathbf{a}^e \tag{2.1}

2. 外力虚功

δaeTqe(2.2)\delta \mathbf{a}^{eT} \mathbf{q}^e \tag{2.2}

3. 内力虚功(单位体积)

δεTσδuTb(2.3)\delta \boldsymbol{\varepsilon}^T \boldsymbol{\sigma} - \delta \mathbf{u}^T \mathbf{b} \tag{2.3}

代入式 (2.1) 得:

δaeT(BTσNTb)(2.4)\delta \mathbf{a}^{eT} \left( \mathbf{B}^T \boldsymbol{\sigma} - \mathbf{N}^T \mathbf{b} \right) \tag{2.4}

4. 平衡条件:外功 = 内功

δaeTqe=δaeT(VeBTσdVVeNTbdV)(2.5)\delta \mathbf{a}^{eT} \mathbf{q}^e = \delta \mathbf{a}^{eT} \left( \int_{V^e} \mathbf{B}^T \boldsymbol{\sigma} \, dV - \int_{V^e} \mathbf{N}^T \mathbf{b} \, dV \right) \tag{2.5}

因为该等式对任意虚位移成立,得:

qe=VeBTσdVVeNTbdV(2.6)\boxed{\mathbf{q}^e = \int_{V^e} \mathbf{B}^T \boldsymbol{\sigma} \, dV - \int_{V^e} \mathbf{N}^T \mathbf{b} \, dV } \tag{2.6}

5. 引入线性本构关系

若采用线性应力-应变关系:

σ=D(εε0)+σ0(2.7)\boldsymbol{\sigma} = \mathbf{D}(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}_0) + \boldsymbol{\sigma}_0 \tag{2.7}

则有:

qe=Keae+fe(2.8)\boxed{ \mathbf{q}^e = \mathbf{K}^e \mathbf{a}^e + \mathbf{f}^e } \tag{2.8}

其中刚度矩阵为:

Ke=VeBTDBdV(2.9a)\boxed{ \mathbf{K}^e = \int_{V^e} \mathbf{B}^T \mathbf{D} \mathbf{B} \, dV } \tag{2.9a}

等效节点力为:

fe=VeNTbdVVeBTDε0dV+VeBTσ0dV(2.9b)\boxed{ \mathbf{f}^e = - \int_{V^e} \mathbf{N}^T \mathbf{b} \, dV - \int_{V^e} \mathbf{B}^T \mathbf{D} \boldsymbol{\varepsilon}_0 \, dV + \int_{V^e} \mathbf{B}^T \boldsymbol{\sigma}_0 \, dV } \tag{2.9b}

总势能原理的变分推导(最小势能原理)(从虚功原理出发)

1 虚功原理变体

当虚位移 δa\delta \mathbf{a}δu\delta \mathbf{u}δε\delta \boldsymbol{\varepsilon} 被视为实际量的变分(variation)时,虚功原理可改写为:

δ(aTr+VuTbdV+AuTtˉdA)=δW(2.10)\delta \left( \mathbf{a}^T \mathbf{r} + \int_V \mathbf{u}^T \mathbf{b} \, dV + \int_A \mathbf{u}^T \bar{\mathbf{t}} \, dA \right) = - \delta W \tag{2.10}

其中,WW外载荷的势能(potential energy of external loads)。当 r,b,tˉ\mathbf{r}, \mathbf{b}, \bar{\mathbf{t}} 与位移无关时(即保守系统),上式成立。

2 内功变分项(应变能)

对于弹性材料,虚功原理中的内力项变为:

δU=VδεTσdV(2.11)\delta U = \int_V \delta \boldsymbol{\varepsilon}^T \boldsymbol{\sigma} \, dV \tag{2.11}

其中 UU 是应变能(strain energy)。对于线弹性材料,代入本构关系:

σ=D(εε0)+σ0\boldsymbol{\sigma} = \mathbf{D} (\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}_0) + \boldsymbol{\sigma}_0

得应变能表达式:

U=12VεTDεdVVεTDε0dV+VεTσ0dV(2.12)U = \frac{1}{2} \int_V \boldsymbol{\varepsilon}^T \mathbf{D} \boldsymbol{\varepsilon} \, dV - \int_V \boldsymbol{\varepsilon}^T \mathbf{D} \boldsymbol{\varepsilon}_0 \, dV + \int_V \boldsymbol{\varepsilon}^T \boldsymbol{\sigma}_0 \, dV \tag{2.12}

D\mathbf{D} 是对称矩阵时,以上变分结果才是正确的,这是单值应变能 UU 存在的必要条件。

3 总势能泛函及其驻值条件

因此,虚功原理可以简洁表示为:

δ(U+W)=δ(Π)=0(2.13)\delta (U + W) = \delta (\Pi) = 0 \tag{2.13}

其中,Π\Pi总势能(total potential energy)。

对于给定位移模式下的有限元系统,令总势能对节点参数 a\mathbf{a} 的偏导为零,即可求解平衡方程:

Πa={Πa1Πa2}=0(2.14)\frac{\partial \Pi}{\partial \mathbf{a}} = \begin{Bmatrix} \frac{\partial \Pi}{\partial a_1} \\ \frac{\partial \Pi}{\partial a_2} \\ \vdots \end{Bmatrix} = \mathbf{0} \tag{2.14}

这就是有限元系统的离散平衡方程。

在稳定弹性系统中,总势能不仅是驻值点(stationary),而且是最小值。
因此,有限元法的本质是:在预设位移模式下,寻找总势能的最小值解。

可以导出和上一节(2.8)的方程

刚度矩阵K具体形式

很容易得到K矩阵的维度为12×1212\times12

随着单元尺寸减小,(2.9a)式可以简化为

Ke=BTDBVe\mathbf{K}^e = \mathbf{B}^\mathrm{T} \mathbf{D} \mathbf{B} V^e

由于B是分块矩阵,故可以写出K的分块形式:

Ke=Ve[BiTDBiBiTDBjBiTDBmBiTDBpBjTDBiBjTDBjBjTDBmBjTDBpBmTDBiBmTDBjBmTDBmBmTDBpBpTDBiBpTDBjBpTDBmBpTDBp]=Ve[K11K12K13K14K21K22K23K24K31K32K33K34K41K42K43K44](2.15)\mathbf{K}^e=V^e \begin{bmatrix} \mathbf{B}_i^\mathrm{T} \mathbf{D} \mathbf{B}_i & \mathbf{B}_i^\mathrm{T} \mathbf{D} \mathbf{B}_j & \mathbf{B}_i^\mathrm{T} \mathbf{D} \mathbf{B}_m & \mathbf{B}_i^\mathrm{T} \mathbf{D} \mathbf{B}_p \\ \mathbf{B}_j^\mathrm{T} \mathbf{D} \mathbf{B}_i & \mathbf{B}_j^\mathrm{T} \mathbf{D} \mathbf{B}_j & \mathbf{B}_j^\mathrm{T} \mathbf{D} \mathbf{B}_m & \mathbf{B}_j^\mathrm{T} \mathbf{D} \mathbf{B}_p \\ \mathbf{B}_m^\mathrm{T} \mathbf{D} \mathbf{B}_i & \mathbf{B}_m^\mathrm{T} \mathbf{D} \mathbf{B}_j & \mathbf{B}_m^\mathrm{T} \mathbf{D} \mathbf{B}_m & \mathbf{B}_m^\mathrm{T} \mathbf{D} \mathbf{B}_p \\ \mathbf{B}_p^\mathrm{T} \mathbf{D} \mathbf{B}_i & \mathbf{B}_p^\mathrm{T} \mathbf{D} \mathbf{B}_j & \mathbf{B}_p^\mathrm{T} \mathbf{D} \mathbf{B}_m & \mathbf{B}_p^\mathrm{T} \mathbf{D} \mathbf{B}_p \end{bmatrix} =V^e \begin{bmatrix} \mathbf{K}_{11} & \mathbf{K}_{12} & \mathbf{K}_{13} & \mathbf{K}_{14}\\ \mathbf{K}_{21} & \mathbf{K}_{22} & \mathbf{K}_{23} & \mathbf{K}_{24}\\ \mathbf{K}_{31} & \mathbf{K}_{32} & \mathbf{K}_{33} & \mathbf{K}_{34}\\ \mathbf{K}_{41} & \mathbf{K}_{42} & \mathbf{K}_{43} & \mathbf{K}_{44}\\ \end{bmatrix} \tag{2.15}

又由于D是对称矩阵,因此有:

(BiTDBj)T=BjTDBi(2.16)(\mathbf{B_i}^\mathrm{T} \mathbf{D} \mathbf{B_j})^\mathrm{T}=\mathbf{B_j}^\mathrm{T} \mathbf{D} \mathbf{B_i} \tag{2.16}

K11\mathbf{K}_{11}:

K11=E72V(2ν2+ν1)[2ν(bi2+ci2+di2)2bi2ci2di2bicibidibici2ν(bi2+ci2+di2)bi22ci2di2cidibidicidi2ν(bi2+ci2+di2)bi2ci22di2](2.17a)\mathbf{K}_{11} = \frac{E}{72 V (2 \nu^2 + \nu - 1)} \begin{bmatrix} 2\nu (b_i^2 + c_i^2 + d_i^2) - 2 b_i^2 - c_i^2 - d_i^2 & - b_i c_i & - b_i d_i \\ - b_i c_i & 2 \nu (b_i^2 + c_i^2 + d_i^2) - b_i^2 - 2 c_i^2 - d_i^2 & - c_i d_i \\ - b_i d_i & - c_i d_i & 2 \nu (b_i^2 + c_i^2 + d_i^2) - b_i^2 - c_i^2 - 2 d_i^2 \end{bmatrix}\tag{2.17a}

K22\mathbf{K}_{22}:

K22=E72V(2ν2+ν1)[2ν(bj2+cj2+dj2)2bj2cj2dj2bjcjbjdjbjcj2ν(bj2+cj2+dj2)bj22cj2dj2cjdjbjdjcjdj2ν(bj2+cj2+dj2)bj2cj22dj2](2.17b)\mathbf{K}_{22} = \frac{E}{72 V (2 \nu^2 + \nu - 1)} \begin{bmatrix} 2\nu (b_j^2 + c_j^2 + d_j^2) - 2 b_j^2 - c_j^2 - d_j^2 & - b_j c_j & - b_j d_j \\ - b_j c_j & 2 \nu (b_j^2 + c_j^2 + d_j^2) - b_j^2 - 2 c_j^2 - d_j^2 & - c_j d_j \\ - b_j d_j & - c_j d_j & 2 \nu (b_j^2 + c_j^2 + d_j^2) - b_j^2 - c_j^2 - 2 d_j^2 \end{bmatrix}\tag{2.17b}

K33\mathbf{K}_{33}:

K33=E72V(2ν2+ν1)[2ν(bm2+cm2+dm2)2bm2cm2dm2bmcmbmdmbmcm2ν(bm2+cm2+dm2)bm22cm2dm2cmdmbmdmcmdm2ν(bm2+cm2+dm2)bm2cm22dm2](2.17c)\mathbf{K}_{33} = \frac{E}{72 V (2 \nu^2 + \nu - 1)} \begin{bmatrix} 2\nu (b_m^2 + c_m^2 + d_m^2) - 2 b_m^2 - c_m^2 - d_m^2 & - b_m c_m & - b_m d_m \\ - b_m c_m & 2 \nu (b_m^2 + c_m^2 + d_m^2) - b_m^2 - 2 c_m^2 - d_m^2 & - c_m d_m \\ - b_m d_m & - c_m d_m & 2 \nu (b_m^2 + c_m^2 + d_m^2) - b_m^2 - c_m^2 - 2 d_m^2 \end{bmatrix}\tag{2.17c}

K44\mathbf{K}_{44}:

K44=E72V(2ν2+ν1)[2ν(bp2+cp2+dp2)2bp2cp2dp2bpcpbpdpbpcp2ν(bp2+cp2+dp2)bp22cp2dp2cpdpbpdpcpdp2ν(bp2+cp2+dp2)bp2cp22dp2](2.17d)\mathbf{K}_{44} = \frac{E}{72 V (2 \nu^2 + \nu - 1)} \begin{bmatrix} 2\nu (b_p^2 + c_p^2 + d_p^2) - 2 b_p^2 - c_p^2 - d_p^2 & - b_p c_p & - b_p d_p \\ - b_p c_p & 2 \nu (b_p^2 + c_p^2 + d_p^2) - b_p^2 - 2 c_p^2 - d_p^2 & - c_p d_p \\ - b_p d_p & - c_p d_p & 2 \nu (b_p^2 + c_p^2 + d_p^2) - b_p^2 - c_p^2 - 2 d_p^2 \end{bmatrix}\tag{2.17d}

K12\mathbf{K}_{12}K21\mathbf{K}_{21}:

K12=K21=E72V(2ν2+ν1)[2bibj+cicj+didj2ν(bibj+cicj+didj)bjci+2νbicj2νbjcibjdi+2νbidj2νbjdibicj2νbicj+2νbjcibibj+2cicj+didj2ν(bibj+cicj+didj)cjdi+2νcidj2νcjdibidj2νbidj+2νbjdicidj2νcidj+2νcjdibibj+cicj+2didj2ν(bibj+cicj+didj)](2.17e)\mathbf{K}_{12} = \mathbf{K}_{21} = - \frac{E}{72 V (2 \nu^2 + \nu - 1)} \begin{bmatrix} 2 b_i b_j + c_i c_j + d_i d_j - 2 \nu (b_i b_j + c_i c_j + d_i d_j) & b_j c_i + 2 \nu b_i c_j - 2 \nu b_j c_i & b_j d_i + 2 \nu b_i d_j - 2 \nu b_j d_i \\[1ex] b_i c_j - 2 \nu b_i c_j + 2 \nu b_j c_i & b_i b_j + 2 c_i c_j + d_i d_j - 2 \nu (b_i b_j + c_i c_j + d_i d_j) & c_j d_i + 2 \nu c_i d_j - 2 \nu c_j d_i \\[1ex] b_i d_j - 2 \nu b_i d_j + 2 \nu b_j d_i & c_i d_j - 2 \nu c_i d_j + 2 \nu c_j d_i & b_i b_j + c_i c_j + 2 d_i d_j - 2 \nu (b_i b_j + c_i c_j + d_i d_j) \end{bmatrix}\tag{2.17e}

K13\mathbf{K}_{13}K31\mathbf{K}_{31}:

K13=K31=E72V(2ν2+ν1)[2bibm+cicm+didm2ν(bibm+cicm+didm)bmci+2νbicm2νbmcibmdi+2νbidm2νbmdibicm2νbicm+2νbmcibibm+2cicm+didm2ν(bibm+cicm+didm)cmdi+2νcidm2νcmdibidm2νbidm+2νbmdicidm2νcidm+2νcmdibibm+cicm+2didm2ν(bibm+cicm+didm)](2.17f)\mathbf{K}_{13} = \mathbf{K}_{31} = - \frac{E}{72 V (2 \nu^2 + \nu - 1)} \begin{bmatrix} 2 b_i b_m + c_i c_m + d_i d_m - 2 \nu (b_i b_m + c_i c_m + d_i d_m) & b_m c_i + 2 \nu b_i c_m - 2 \nu b_m c_i & b_m d_i + 2 \nu b_i d_m - 2 \nu b_m d_i \\[1ex] b_i c_m - 2 \nu b_i c_m + 2 \nu b_m c_i & b_i b_m + 2 c_i c_m + d_i d_m - 2 \nu (b_i b_m + c_i c_m + d_i d_m) & c_m d_i + 2 \nu c_i d_m - 2 \nu c_m d_i \\[1ex] b_i d_m - 2 \nu b_i d_m + 2 \nu b_m d_i & c_i d_m - 2 \nu c_i d_m + 2 \nu c_m d_i & b_i b_m + c_i c_m + 2 d_i d_m - 2 \nu (b_i b_m + c_i c_m + d_i d_m) \end{bmatrix}\tag{2.17f}

K14\mathbf{K}_{14}K41\mathbf{K}_{41}:

K14=K41=E72V(2ν2+ν1)[2bibp+cicp+didp2ν(bibp+cicp+didp)bpci+2νbicp2νbpcibpdi+2νbidp2νbpdibicp2νbicp+2νbpcibibp+2cicp+didp2ν(bibp+cicp+didp)cpdi+2νcidp2νcpdibidp2νbidp+2νbpdicidp2νcidp+2νcpdibibp+cicp+2didp2ν(bibp+cicp+didp)](2.17g)\mathbf{K}_{14} = \mathbf{K}_{41} = - \frac{E}{72 V (2 \nu^2 + \nu - 1)} \begin{bmatrix} 2 b_i b_p + c_i c_p + d_i d_p - 2 \nu (b_i b_p + c_i c_p + d_i d_p) & b_p c_i + 2 \nu b_i c_p - 2 \nu b_p c_i & b_p d_i + 2 \nu b_i d_p - 2 \nu b_p d_i \\[1ex] b_i c_p - 2 \nu b_i c_p + 2 \nu b_p c_i & b_i b_p + 2 c_i c_p + d_i d_p - 2 \nu (b_i b_p + c_i c_p + d_i d_p) & c_p d_i + 2 \nu c_i d_p - 2 \nu c_p d_i \\[1ex] b_i d_p - 2 \nu b_i d_p + 2 \nu b_p d_i & c_i d_p - 2 \nu c_i d_p + 2 \nu c_p d_i & b_i b_p + c_i c_p + 2 d_i d_p - 2 \nu (b_i b_p + c_i c_p + d_i d_p) \end{bmatrix}\tag{2.17g}

K23\mathbf{K}_{23}K32\mathbf{K}_{32}:

K23=K32=E72V(2ν2+ν1)[2bjbm+cjcm+djdm2ν(bjbm+cjcm+djdm)bmcj+2νbjcm2νbmcjbmdj+2νbjdm2νbmdjbjcm2νbjcm+2νbmcjbjbm+2cjcm+djdm2ν(bjbm+cjcm+djdm)cmdj+2νcjdm2νcmdjbjdm2νbjdm+2νbmdjcjdm2νcjdm+2νcmdjbjbm+cjcm+2djdm2ν(bjbm+cjcm+djdm)](2.17h)\mathbf{K}_{23} = \mathbf{K}_{32} = - \frac{E}{72 V (2 \nu^2 + \nu - 1)} \begin{bmatrix} 2 b_j b_m + c_j c_m + d_j d_m - 2 \nu (b_j b_m + c_j c_m + d_j d_m) & b_m c_j + 2 \nu b_j c_m - 2 \nu b_m c_j & b_m d_j + 2 \nu b_j d_m - 2 \nu b_m d_j \\[1ex] b_j c_m - 2 \nu b_j c_m + 2 \nu b_m c_j & b_j b_m + 2 c_j c_m + d_j d_m - 2 \nu (b_j b_m + c_j c_m + d_j d_m) & c_m d_j + 2 \nu c_j d_m - 2 \nu c_m d_j \\[1ex] b_j d_m - 2 \nu b_j d_m + 2 \nu b_m d_j & c_j d_m - 2 \nu c_j d_m + 2 \nu c_m d_j & b_j b_m + c_j c_m + 2 d_j d_m - 2 \nu (b_j b_m + c_j c_m + d_j d_m) \end{bmatrix}\tag{2.17h}

K24\mathbf{K}_{24}K42\mathbf{K}_{42}:

K24=K42=E72V(2ν2+ν1)[2bjbp+cjcp+djdp2ν(bjbp+cjcp+djdp)bpcj+2νbjcp2νbpcjbpdj+2νbjdp2νbpdjbjcp2νbjcp+2νbpcjbjbp+2cjcp+djdp2ν(bjbp+cjcp+djdp)cpdj+2νcjdp2νcpdjbjdp2νbjdp+2νbpdjcjdp2νcjdp+2νcpdjbjbp+cjcp+2djdp2ν(bjbp+cjcp+djdp)](2.17i)\mathbf{K}_{24} = \mathbf{K}_{42} = - \frac{E}{72 V (2 \nu^2 + \nu - 1)} \begin{bmatrix} 2 b_j b_p + c_j c_p + d_j d_p - 2 \nu (b_j b_p + c_j c_p + d_j d_p) & b_p c_j + 2 \nu b_j c_p - 2 \nu b_p c_j & b_p d_j + 2 \nu b_j d_p - 2 \nu b_p d_j \\[1ex] b_j c_p - 2 \nu b_j c_p + 2 \nu b_p c_j & b_j b_p + 2 c_j c_p + d_j d_p - 2 \nu (b_j b_p + c_j c_p + d_j d_p) & c_p d_j + 2 \nu c_j d_p - 2 \nu c_p d_j \\[1ex] b_j d_p - 2 \nu b_j d_p + 2 \nu b_p d_j & c_j d_p - 2 \nu c_j d_p + 2 \nu c_p d_j & b_j b_p + c_j c_p + 2 d_j d_p - 2 \nu (b_j b_p + c_j c_p + d_j d_p) \end{bmatrix}\tag{2.17i}

K34\mathbf{K}_{34}K43\mathbf{K}_{43}:

K34=K43=E72V(2ν2+ν1)[2bmbp+cmcp+dmdp2ν(bmbp+cmcp+dmdp)bpcm+2νbmcp2νbpcmbpdm+2νbmdp2νbpdmbmcp2νbmcp+2νbpcmbmbp+2cmcp+dmdp2ν(bmbp+cmcp+dmdp)cpdm+2νcmdp2νcpdmbmdp2νbmdp+2νbpdmcmdp2νcmdp+2νcpdmbmbp+cmcp+2dmdp2ν(bmbp+cmcp+dmdp)](2.17j)\mathbf{K}_{34} = \mathbf{K}_{43} = - \frac{E}{72 V (2 \nu^2 + \nu - 1)} \begin{bmatrix} 2 b_m b_p + c_m c_p + d_m d_p - 2 \nu (b_m b_p + c_m c_p + d_m d_p) & b_p c_m + 2 \nu b_m c_p - 2 \nu b_p c_m & b_p d_m + 2 \nu b_m d_p - 2 \nu b_p d_m \\[1ex] b_m c_p - 2 \nu b_m c_p + 2 \nu b_p c_m & b_m b_p + 2 c_m c_p + d_m d_p - 2 \nu (b_m b_p + c_m c_p + d_m d_p) & c_p d_m + 2 \nu c_m d_p - 2 \nu c_p d_m \\[1ex] b_m d_p - 2 \nu b_m d_p + 2 \nu b_p d_m & c_m d_p - 2 \nu c_m d_p + 2 \nu c_p d_m & b_m b_p + c_m c_p + 2 d_m d_p - 2 \nu (b_m b_p + c_m c_p + d_m d_p) \end{bmatrix}\tag{2.17j}

该封面图片由Artur PawlakPixabay上发布