三维刚度矩阵K的详细推导
前置弹性力学矩阵知识
几何变换矩阵S
S=⎣⎢⎢⎢⎢⎢⎢⎢⎡∂x∂00∂y∂0∂z∂0∂y∂0∂x∂∂z∂000∂z∂0∂y∂∂x∂⎦⎥⎥⎥⎥⎥⎥⎥⎤(1.1)
位移矩阵u
u=⎩⎨⎧uvw⎭⎬⎫(1.2)
应变矩阵ϵ
由(1.1)式和(1.2)式导出:
ε=⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧εxεyεzγxyγyzγzx⎭⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎫=⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧∂x∂u∂y∂v∂z∂w∂y∂u+∂x∂v∂z∂v+∂y∂w∂x∂w+∂z∂u⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫=Su(1.3)
节点位移a
aeai=⎩⎪⎪⎨⎪⎪⎧aiajamap⎭⎪⎪⎬⎪⎪⎫=⎩⎨⎧uiviwi⎭⎬⎫etc.(1.4)
形函数N
设任意点位移为
uvw=α1+α2x+α3y+α4z=α5+α6x+α7y+α8z=α9+α10x+α11y+α12z(1.5)
每个节点的空间位置也已知:
Xi=(xi,yi,zi),etc(1.6)
对于u,可以得到
uiujumup=α1+α2xi+α3yi+α4zi=α1+α2xj+α3yj+α4zj=α1+α2xm+α3ym+α4zm=α1+α2xp+α3yp+α4zp(1.7)
写成矩阵形式
X∈R4×4⎣⎢⎢⎡1111xixjxmxpyiyjymypzizjzmzp⎦⎥⎥⎤⋅⎣⎢⎢⎡α1α2α3α4⎦⎥⎥⎤=⎣⎢⎢⎡uiujumup⎦⎥⎥⎤(1.8)
可以解出
uvw=6V1[(ai+bix+ciy+diz)ui+(aj+bjx+cjy+djz)uj+(am+bmx+cmy+dmz)um+(ap+bpx+cpy+dpz)up]=6V1[(ai+bix+ciy+diz)vi+(aj+bjx+cjy+djz)vj+(am+bmx+cmy+dmz)vm+(ap+bpx+cpy+dpz)vp]=6V1[(ai+bix+ciy+diz)wi+(aj+bjx+cjy+djz)wj+(am+bmx+cmy+dmz)wm+(ap+bpx+cpy+dpz)wp](1.9)
其中
6V=det∣∣∣∣∣∣∣∣1111xixjxmxpyiyjymypzizjzmzp∣∣∣∣∣∣∣∣(1.10)
aibicidi=det∣∣∣∣∣∣xjxmxpyjymypzjzmzp∣∣∣∣∣∣=−det∣∣∣∣∣∣111yjymypzjzmzp∣∣∣∣∣∣=−det∣∣∣∣∣∣xjxmxp111zjzmzp∣∣∣∣∣∣=−det∣∣∣∣∣∣xjxmxpyjymyp111∣∣∣∣∣∣(1.11)
其他的系数由更换角标得到,按照右手准则按顺序变换。
可以得到形函数Ni:
Ni=6V1⎣⎡ai+bix+ciy+diz000ai+bix+ciy+diz000ai+bix+ciy+diz⎦⎤(1.12)
形函数插值表达式u=Nae
ε=SNae=Bae=[BiBjBmBp]ae(1.13)
应变-位移矩阵B
由(1.13)式导出,即B=SN
Bi=⎣⎢⎢⎢⎢⎢⎢⎢⎡∂x∂Ni00∂y∂Ni0∂z∂Ni0∂y∂Ni0∂x∂Ni∂z∂Ni000∂z∂Ni0∂y∂Ni∂x∂Ni⎦⎥⎥⎥⎥⎥⎥⎥⎤=6V1⎣⎢⎢⎢⎢⎢⎢⎡bi00ci0di0ci0bidi000di0cibi⎦⎥⎥⎥⎥⎥⎥⎤(1.14)
其他块直接更换角标即可。
各向同性弹性矩阵D
D=(1+ν)(1−2ν)E⎣⎢⎢⎢⎢⎢⎢⎡1−ννν000ν1−νν000νν1−ν00000021−2ν00000021−2ν00000021−2ν⎦⎥⎥⎥⎥⎥⎥⎤(1.15)
刚度矩阵K表达式推导
基于虚功原理的推导
为使节点力在静力上等效于真实边界应力与体力,最简单的方法是对节点施加一个任意的虚位移 δae,然后使该虚位移下的外功等于内功。
1. 虚位移与应变表达式
δu=Nδae,δε=Bδae(2.1)
2. 外力虚功
δaeTqe(2.2)
3. 内力虚功(单位体积)
δεTσ−δuTb(2.3)
代入式 (2.1) 得:
δaeT(BTσ−NTb)(2.4)
4. 平衡条件:外功 = 内功
δaeTqe=δaeT(∫VeBTσdV−∫VeNTbdV)(2.5)
因为该等式对任意虚位移成立,得:
qe=∫VeBTσdV−∫VeNTbdV(2.6)
5. 引入线性本构关系
若采用线性应力-应变关系:
σ=D(ε−ε0)+σ0(2.7)
则有:
qe=Keae+fe(2.8)
其中刚度矩阵为:
Ke=∫VeBTDBdV(2.9a)
等效节点力为:
fe=−∫VeNTbdV−∫VeBTDε0dV+∫VeBTσ0dV(2.9b)
总势能原理的变分推导(最小势能原理)(从虚功原理出发)
1 虚功原理变体
当虚位移 δa、δu、δε 被视为实际量的变分(variation)时,虚功原理可改写为:
δ(aTr+∫VuTbdV+∫AuTtˉdA)=−δW(2.10)
其中,W 是外载荷的势能(potential energy of external loads)。当 r,b,tˉ 与位移无关时(即保守系统),上式成立。
2 内功变分项(应变能)
对于弹性材料,虚功原理中的内力项变为:
δU=∫VδεTσdV(2.11)
其中 U 是应变能(strain energy)。对于线弹性材料,代入本构关系:
σ=D(ε−ε0)+σ0
得应变能表达式:
U=21∫VεTDεdV−∫VεTDε0dV+∫VεTσ0dV(2.12)
当 D 是对称矩阵时,以上变分结果才是正确的,这是单值应变能 U 存在的必要条件。
3 总势能泛函及其驻值条件
因此,虚功原理可以简洁表示为:
δ(U+W)=δ(Π)=0(2.13)
其中,Π 是总势能(total potential energy)。
对于给定位移模式下的有限元系统,令总势能对节点参数 a 的偏导为零,即可求解平衡方程:
∂a∂Π=⎩⎪⎨⎪⎧∂a1∂Π∂a2∂Π⋮⎭⎪⎬⎪⎫=0(2.14)
这就是有限元系统的离散平衡方程。
在稳定弹性系统中,总势能不仅是驻值点(stationary),而且是最小值。
因此,有限元法的本质是:在预设位移模式下,寻找总势能的最小值解。
可以导出和上一节(2.8)的方程
刚度矩阵K具体形式
很容易得到K矩阵的维度为12×12。
随着单元尺寸减小,(2.9a)式可以简化为
Ke=BTDBVe
由于B是分块矩阵,故可以写出K的分块形式:
Ke=Ve⎣⎢⎢⎡BiTDBiBjTDBiBmTDBiBpTDBiBiTDBjBjTDBjBmTDBjBpTDBjBiTDBmBjTDBmBmTDBmBpTDBmBiTDBpBjTDBpBmTDBpBpTDBp⎦⎥⎥⎤=Ve⎣⎢⎢⎡K11K21K31K41K12K22K32K42K13K23K33K43K14K24K34K44⎦⎥⎥⎤(2.15)
又由于D是对称矩阵,因此有:
(BiTDBj)T=BjTDBi(2.16)
K11:
K11=72V(2ν2+ν−1)E⎣⎡2ν(bi2+ci2+di2)−2bi2−ci2−di2−bici−bidi−bici2ν(bi2+ci2+di2)−bi2−2ci2−di2−cidi−bidi−cidi2ν(bi2+ci2+di2)−bi2−ci2−2di2⎦⎤(2.17a)
K22:
K22=72V(2ν2+ν−1)E⎣⎡2ν(bj2+cj2+dj2)−2bj2−cj2−dj2−bjcj−bjdj−bjcj2ν(bj2+cj2+dj2)−bj2−2cj2−dj2−cjdj−bjdj−cjdj2ν(bj2+cj2+dj2)−bj2−cj2−2dj2⎦⎤(2.17b)
K33:
K33=72V(2ν2+ν−1)E⎣⎡2ν(bm2+cm2+dm2)−2bm2−cm2−dm2−bmcm−bmdm−bmcm2ν(bm2+cm2+dm2)−bm2−2cm2−dm2−cmdm−bmdm−cmdm2ν(bm2+cm2+dm2)−bm2−cm2−2dm2⎦⎤(2.17c)
K44:
K44=72V(2ν2+ν−1)E⎣⎡2ν(bp2+cp2+dp2)−2bp2−cp2−dp2−bpcp−bpdp−bpcp2ν(bp2+cp2+dp2)−bp2−2cp2−dp2−cpdp−bpdp−cpdp2ν(bp2+cp2+dp2)−bp2−cp2−2dp2⎦⎤(2.17d)
K12和K21:
K12=K21=−72V(2ν2+ν−1)E⎣⎢⎡2bibj+cicj+didj−2ν(bibj+cicj+didj)bicj−2νbicj+2νbjcibidj−2νbidj+2νbjdibjci+2νbicj−2νbjcibibj+2cicj+didj−2ν(bibj+cicj+didj)cidj−2νcidj+2νcjdibjdi+2νbidj−2νbjdicjdi+2νcidj−2νcjdibibj+cicj+2didj−2ν(bibj+cicj+didj)⎦⎥⎤(2.17e)
K13和K31:
K13=K31=−72V(2ν2+ν−1)E⎣⎢⎡2bibm+cicm+didm−2ν(bibm+cicm+didm)bicm−2νbicm+2νbmcibidm−2νbidm+2νbmdibmci+2νbicm−2νbmcibibm+2cicm+didm−2ν(bibm+cicm+didm)cidm−2νcidm+2νcmdibmdi+2νbidm−2νbmdicmdi+2νcidm−2νcmdibibm+cicm+2didm−2ν(bibm+cicm+didm)⎦⎥⎤(2.17f)
K14 和 K41:
K14=K41=−72V(2ν2+ν−1)E⎣⎢⎡2bibp+cicp+didp−2ν(bibp+cicp+didp)bicp−2νbicp+2νbpcibidp−2νbidp+2νbpdibpci+2νbicp−2νbpcibibp+2cicp+didp−2ν(bibp+cicp+didp)cidp−2νcidp+2νcpdibpdi+2νbidp−2νbpdicpdi+2νcidp−2νcpdibibp+cicp+2didp−2ν(bibp+cicp+didp)⎦⎥⎤(2.17g)
K23 和 K32:
K23=K32=−72V(2ν2+ν−1)E⎣⎢⎡2bjbm+cjcm+djdm−2ν(bjbm+cjcm+djdm)bjcm−2νbjcm+2νbmcjbjdm−2νbjdm+2νbmdjbmcj+2νbjcm−2νbmcjbjbm+2cjcm+djdm−2ν(bjbm+cjcm+djdm)cjdm−2νcjdm+2νcmdjbmdj+2νbjdm−2νbmdjcmdj+2νcjdm−2νcmdjbjbm+cjcm+2djdm−2ν(bjbm+cjcm+djdm)⎦⎥⎤(2.17h)
K24 和 K42:
K24=K42=−72V(2ν2+ν−1)E⎣⎢⎡2bjbp+cjcp+djdp−2ν(bjbp+cjcp+djdp)bjcp−2νbjcp+2νbpcjbjdp−2νbjdp+2νbpdjbpcj+2νbjcp−2νbpcjbjbp+2cjcp+djdp−2ν(bjbp+cjcp+djdp)cjdp−2νcjdp+2νcpdjbpdj+2νbjdp−2νbpdjcpdj+2νcjdp−2νcpdjbjbp+cjcp+2djdp−2ν(bjbp+cjcp+djdp)⎦⎥⎤(2.17i)
K34 和 K43:
K34=K43=−72V(2ν2+ν−1)E⎣⎢⎡2bmbp+cmcp+dmdp−2ν(bmbp+cmcp+dmdp)bmcp−2νbmcp+2νbpcmbmdp−2νbmdp+2νbpdmbpcm+2νbmcp−2νbpcmbmbp+2cmcp+dmdp−2ν(bmbp+cmcp+dmdp)cmdp−2νcmdp+2νcpdmbpdm+2νbmdp−2νbpdmcpdm+2νcmdp−2νcpdmbmbp+cmcp+2dmdp−2ν(bmbp+cmcp+dmdp)⎦⎥⎤(2.17j)
该封面图片由Artur Pawlak在Pixabay上发布