平面等参四边形单元
欢迎阅读原文：《有限元基础编程百科全书》| 平面等参四边形单元
由于实际问题的复杂性，需要使用一些几何形状不太规整的单元来逼近原问题，特别是在一些复杂的边界上，有时只能采用不规整单元。
但直接研究这些不规整单元则比较困难，如何利用几何规整单元（如三角形单元、矩形单元、正六面体单元）的结果来研究(推导)所对应的几何不规整单元（叫做参数单元(parametric element)）的表达式？这将涉及到几何形状映射、坐标系变换（等参变换、非等参变换）等问题。
有些教材会引导读者先学习矩形单元，然后逐渐引入四边形单元的等参变换思想。本推文为方便起见直接向读者介绍何为等参变换？如何等参变换？等参变换的意义何在？后续所有单元的刚度矩阵均由等参思想推导。
等参变换思想
如下图所示，左侧为等参单元在自然坐标系下的表示，为一个规则的正方形单元(母单元)，边长为2；右侧为物理坐标系下的任意形状四边形单元（子单元）。两者坐标系通过特殊的变换（雅可比矩阵），即可由其中一个坐标系代表另一个坐标系。
以形函数对坐标轴偏导为例：
[ N I , ξ Q 4 N I , η Q 4 ] = [ x , ξ y , ξ x , η y , η ] [ N I , x Q 4 N I , y Q 4 ] \begin{bmatrix}N_{I,\xi}^{Q4}\\N_{I,\eta}^{Q4}\end{bmatrix}=\begin{bmatrix}x_{,\xi}&y_{,\xi}\\x_{,\eta}&y_{,\eta}\end{bmatrix}\begin{bmatrix}N_{I,x}^{Q4}\\N_{I,y}^{Q4}\end{bmatrix} [NI,ξQ4NI,ηQ4]=[x,ξx,ηy,ξy,η][NI,xQ4NI,yQ4]
其中，雅可比矩阵可表示为：
J e = [ x , ξ y , ξ x , η y , η ] \boldsymbol{J}^e=\begin{bmatrix}x_{,\xi}&y_{,\xi}\\x_{,\eta}&y_{,\eta}\end{bmatrix} Je=[x,ξx,ηy,ξy,η]
写到这里，我们可以看到在实际计算过程中为什么要进行等参变换？
它可以将任意形状的单元的形函数都归结于规则单元形状的形函数，大大降低了求形函数的计算成本，因为规则等参单元的形函数是由前人推导而来，是一组确定的量。
通过等参变换可以将已知等参单元的形函数表示任意单元形状的形函数。（用词可能不太严谨，该段解释是我用自己的语言表述，若有错误之处，欢迎批评指正！）
与等参单元概念类似的还有亚参单元(subparametric element)、超参单元(superparametric element)。
- 等参单元：单元内坐标和位移可以通过相同的形函数分别由节点坐标和节点位移插值得到；
- 亚参单元：单元几何映射时使用的节点数小于构建试探解时使用的节点数；
- 超参单元：单元几何映射时使用的节点数大于构建试探解时使用的节点数；
以上概念有个印象就行，目前本文只针对等参单元进行数值编程。
单元介绍
线性应变四边形单元具有4个节点，每个节点两个平动自由度。单元描述如下图所示，节点序号按照逆时针排列（顺序无所谓，但是要一致），单元共有8个自由度，记作 q e \mathbf{q}_e qe，相应的有8个节点力，记作 P e \mathbf{P}_e Pe：
q ( 8 × 1 ) e = [ u 1 v 1 u 2 v 2 u 3 v 3 u 4 v 4 ] T \mathbf{q}_{(8\times1)}^e=\begin{bmatrix}u_1&v_1&u_2&v_2&u_3&v_3&u_4&v_4\end{bmatrix}^T q(8×1)e=[u1v1u2v2u3v3u4v4]T
P ( 8 × 1 ) e = [ P x 1 P y 1 P x 2 P y 2 P x 3 P y 3 P x 4 P y 4 ] T \mathbf{P}_{(8\times1)}^e=\begin{bmatrix}P_{x1}&P_{y1}&P_{x2}&P_{y2}&P_{x3}&P_{y3}&P_{x4}&P_{y4}\end{bmatrix}^T P(8×1)e=[Px1Py1Px2Py2Px3Py3Px4Py4]T
有限元格式
接下来会带着大家再次熟悉有限元基本格式，在以后的有限元问题描述上我都会不断的串讲有限元基本格式，加深印象。
因为格式真的很固定，特别适合计算机处理，当然这不代表有限元简单，只不过很多内容不需要我们去推倒，直接拿来用就行（应用型教学）。
以下的理论读者不要觉得枯燥，是和编程直接相关的理论基础，本书将与编程无关的理论均剔除，在后面的程序中，以下的理论公式均有所体现。
形函数
{ N 1 ( ξ , η ) = 1 4 ( 1 − ξ ) ( 1 − η ) N 2 ( ξ , η ) = 1 4 ( 1 + ξ ) ( 1 − η ) N 3 ( ξ , η ) = 1 4 ( 1 + ξ ) ( 1 + η ) N 4 ( ξ , η ) = 1 4 ( 1 − ξ ) ( 1 + η ) \begin{cases} N_1(\xi ,\eta )=\frac{1}{4}(1-\xi )(1-\eta )\\ N_2(\xi ,\eta )=\frac{1}{4}(1+\xi )(1-\eta )\\ N_3(\xi ,\eta )=\frac{1}{4}(1+\xi )(1+\eta )\\ N_4(\xi ,\eta )=\frac{1}{4}(1-\xi )(1+\eta )\\ \end{cases} ⎩ ⎨ ⎧N1(ξ,η)=41(1−ξ)(1−η)N2(ξ,η)=41(1+ξ)(1−η)N3(ξ,η)=41(1+ξ)(1+η)N4(ξ,η)=41(1−ξ)(1+η)
位移场
单元内某一点的位移，可通过形函数与节点位移描述。
单元内位移场可描述为：
{
u
=
N
1
u
1
+
N
2
u
2
+
N
3
u
3
+
N
4
u
4
v
=
N
1
v
1
+
N
2
v
2
+
N
3
v
3
+
N
4
v
4
\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}
{u=N1u1+N2u2+N3u3+N4u4v=N1v1+N2v2+N3v3+N4v4
矩阵表示法：
{
u
v
}
=
[
N
1
0
N
2
0
N
3
0
N
4
0
0
N
1
0
N
2
0
N
3
0
N
4
]
{
u
1
v
1
u
2
v
2
u
3
v
4
}
\left\{ \begin{array}{c} u\\ v\\ \end{array} \right\} =\left[ \begin{matrix} N_1& 0& N_2& 0& N_3& 0& N_4& 0\\ 0& N_1& 0& N_2& 0& N_3& 0& N_4\\ \end{matrix} \right] \left\{ \begin{array}{c} u_1\\ v_1\\ u_2\\ v_2\\ u_3\\ v_4\\ \end{array} \right\}
{uv}=[N100N1N200N2N300N3N400N4]⎩
⎨
⎧u1v1u2v2u3v4⎭
⎬
⎫
简写： { U } = [ N ] { q e } \{\mathbf{U}\}=[\mathbf{N}]\{\mathbf{q}^e\} {U}=[N]{qe}
形函数
N
i
(
ξ
,
η
)
N_{i}(\xi,\eta)
Ni(ξ,η)也可以反映母单元和子单元的映射关系，设母单元内部一点的坐标为
(
x
,
y
)
(x,y)
(x,y)，可由形函数插值得到：
{
x
=
N
1
x
1
+
N
2
x
2
+
N
3
x
3
+
N
4
x
4
y
=
N
1
y
1
+
N
2
y
2
+
N
3
y
3
+
N
4
y
4
\begin{cases} x=N_1x_1+N_2x_2+N_3x_3+N_4x_4\\ y=N_1y_1+N_2y_2+N_3y_3+N_4y_4\\ \end{cases}
{x=N1x1+N2x2+N3x3+N4x4y=N1y1+N2y2+N3y3+N4y4
雅可比矩阵可表示为：
[ J ] = [ ∂ x ∂ ξ ∂ y ∂ ξ ∂ x ∂ η ∂ y ∂ η ] = [ ∑ i = 1 4 ∂ N i ∂ ξ x i ∑ i = 1 4 ∂ N i ∂ ξ y i ∑ i = 1 4 ∂ N i ∂ η x i ∑ i = 1 4 ∂ N i ∂ η y i ] [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}\sum_{i=1}^4\frac{\partial N_i}{\partial\xi}x_i&\sum_{i=1}^4\frac{\partial N_i}{\partial\xi}y_i\\\sum_{i=1}^4\frac{\partial N_i}{\partial\eta}x_i&\sum_{i=1}^4\frac{\partial N_i}{\partial\eta}y_i\end{bmatrix} [J]=[∂ξ∂x∂η∂x∂ξ∂y∂η∂y]=[∑i=14∂ξ∂Nixi∑i=14∂η∂Nixi∑i=14∂ξ∂Niyi∑i=14∂η∂Niyi]
不要被上面的表达式所吓倒哈！其实在数值编程的过程中，可以把它换成另一种形式：
[
J
]
=
1
4
[
−
(
1
−
η
)
(
1
−
η
)
(
1
+
η
)
−
(
1
+
η
)
−
(
1
−
ξ
)
−
(
1
+
ξ
)
(
1
+
ξ
)
(
1
−
ξ
)
]
[
x
1
y
1
x
2
y
2
x
3
y
3
x
4
y
4
]
[J]=\frac{1}{4}\begin{bmatrix}-(1-\eta)&(1-\eta)&(1+\eta)&-(1+\eta)\\[1ex]-(1-\xi)&-(1+\xi)&(1+\xi)&(1-\xi)\end{bmatrix}\begin{bmatrix}x_1&y_1\\[1ex]x_2&y_2\\[1ex]x_3&y_3\\[1ex]x_4&y_4\end{bmatrix}
[J]=41[−(1−η)−(1−ξ)(1−η)−(1+ξ)(1+η)(1+ξ)−(1+η)(1−ξ)]
x1x2x3x4y1y2y3y4
如此以来是不是简单许多，左侧项为形函数的偏导矩阵，右侧项为节点坐标。
function [JacobianMatrix,XYDerivatives] = Jacobian(nodeCoordinates,naturalDerivatives)
JacobianMatrix = nodeCoordinates'*naturalDerivatives;
XYDerivatives = naturalDerivatives/JacobianMatrix;
end
几何方程
单元内某一点的应变，可通过几何矩阵 B B B与节点位移描述。
ε
(
3
×
1
)
(
x
,
y
)
=
[
∂
]
(
3
×
2
)
N
(
x
,
y
)
⋅
q
e
=
B
(
3
×
8
)
(
x
,
y
)
⋅
q
(
8
×
1
)
e
\boldsymbol{\underset{(3\times1)}{\varepsilon}}(x,y)=\begin{bmatrix}\partial\end{bmatrix}_{(3\times2)}\mathbf{N}\begin{pmatrix}x,y\end{pmatrix}\cdotp\mathbf{q}^e=\mathbf{B}_{(3\times8)}\begin{pmatrix}x,y\end{pmatrix}\cdotp\mathbf{q}_{(8\times1)}^e
(3×1)ε(x,y)=[∂](3×2)N(x,y)⋅qe=B(3×8)(x,y)⋅q(8×1)e
其中，
B
(
3
×
8
)
(
x
,
y
)
=
[
∂
]
(
3
×
2
)
N
(
2
×
8
)
=
[
∂
∂
x
0
0
∂
∂
y
∂
∂
y
∂
∂
x
]
[
N
1
0
N
2
0
N
3
0
N
4
0
0
N
1
0
N
2
0
N
3
0
N
4
]
\mathbf{B}_{(3\times8)}\big(x,y\big)=\begin{bmatrix}\partial\end{bmatrix}_{(3\times2)}\mathbf{N}_{(2\times8)}=\begin{bmatrix}\frac{\partial}{\partial x}&0\\\\0&\frac{\partial}{\partial y}\\\\\frac{\partial}{\partial y}&\frac{\partial}{\partial x}\end{bmatrix}\begin{bmatrix}N_1&0&N_2&0&N_3&0&N_4&0\\0&N_1&0&N_2&0&N_3&0&N_4\end{bmatrix}
B(3×8)(x,y)=[∂](3×2)N(2×8)=
∂x∂0∂y∂0∂y∂∂x∂
[N100N1N200N2N300N3N400N4]
再次展开：
[
B
]
=
[
∂
N
1
∂
x
0
∂
N
2
∂
x
0
∂
N
3
∂
x
0
∂
N
4
∂
x
0
0
∂
N
1
∂
y
0
∂
N
2
∂
y
0
∂
N
3
∂
y
0
∂
N
4
∂
y
∂
N
1
∂
y
∂
N
1
∂
x
∂
N
2
∂
y
∂
N
2
∂
x
∂
N
3
∂
y
∂
N
3
∂
x
∂
N
4
∂
y
∂
N
4
∂
x
]
[B]=\left[ \begin{matrix} \frac{\partial N_1}{\partial x}& 0& \frac{\partial N_2}{\partial x}& 0& \frac{\partial N_3}{\partial x}& 0& \frac{\partial N_4}{\partial x}& 0\\ 0& \frac{\partial N_1}{\partial y}& 0& \frac{\partial N_2}{\partial y}& 0& \frac{\partial N_3}{\partial y}& 0& \frac{\partial N_4}{\partial y}\\ \frac{\partial N_1}{\partial y}& \frac{\partial N_1}{\partial x}& \frac{\partial N_2}{\partial y}& \frac{\partial N_2}{\partial x}& \frac{\partial N_3}{\partial y}& \frac{\partial N_3}{\partial x}& \frac{\partial N_4}{\partial y}& \frac{\partial N_4}{\partial x}\\ \end{matrix} \right]
[B]=
∂x∂N10∂y∂N10∂y∂N1∂x∂N1∂x∂N20∂y∂N20∂y∂N2∂x∂N2∂x∂N30∂y∂N30∂y∂N3∂x∂N3∂x∂N40∂y∂N40∂y∂N4∂x∂N4
以上的公式中体现的是形函数对物理坐标轴的偏导，这时候就可以用等参变换关系，用母单元的形函数偏导来表示。
{ ∂ N i ∂ ξ ∂ N i ∂ η } = [ ∂ x ∂ ξ ∂ y ∂ ξ ∂ x ∂ η ∂ y ∂ η ] { ∂ N i ∂ x ∂ N i ∂ y } \begin{Bmatrix}\frac{\partial N_i}{\partial\xi}\\\frac{\partial N_i}{\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 N_i}{\partial x}\\\frac{\partial N_i}{\partial y}\end{Bmatrix} {∂ξ∂Ni∂η∂Ni}=[∂ξ∂x∂η∂x∂ξ∂y∂η∂y]{∂x∂Ni∂y∂Ni}
进一步变换：
{ ∂ N i ∂ x ∂ N i ∂ y } = [ J ] − 1 [ ∂ x ∂ ξ ∂ y ∂ ξ ∂ x ∂ η ∂ y ∂ η ] \begin{Bmatrix}\frac{\partial N_i}{\partial x}\\\frac{\partial N_i}{\partial y}\end{Bmatrix}=[J]^{-1}\begin{bmatrix}\frac{\partial x}{\partial\xi}&\frac{\partial y}{\partial\xi}\\\frac{\partial x}{\partial\eta}&\frac{\partial y}{\partial\eta}\end{bmatrix} {∂x∂Ni∂y∂Ni}=[J]−1[∂ξ∂x∂η∂x∂ξ∂y∂η∂y]
本构方程
针对线弹性各向同性材料，平面应力状态下：
σ ( 3 × 1 ) ( x , y ) = [ σ x x σ y y τ x y ] = E 1 − μ 2 [ 1 μ 0 μ 1 0 0 0 1 − μ 2 ] [ ε x x ε y y γ x y ] = D ( 3 × 3 ) ⋅ ε ( 3 × 1 ) \boldsymbol{\sigma}_{(3\times1)}\left(x,y\right)=\begin{bmatrix}\sigma_{xx}\\\sigma_{yy}\\\tau_{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_{xx}\\\varepsilon_{yy}\\\gamma_{xy}\end{bmatrix}=\mathbf{D}_{(3\times3)}\cdot\mathbf{\varepsilon}_{(3\times1)} σ(3×1)(x,y)= σxxσyyτxy =1−μ2E 1μ0μ100021−μ εxxεyyγxy =D(3×3)⋅ε(3×1)
平面应变状态：将上式中 ( E , μ ) (E,\mu) (E,μ)的换为 ( E 1 − μ 2 , μ 1 − μ ) (\frac E{1-\mu^2},\frac\mu{1-\mu}) (1−μ2E,1−μμ)即可。
单元刚度方程
[ K ] = [ ∫ A e [ B ] T [ D ] [ B ] t d A ] [K]=\left[\int_{A_e}[B]^\mathrm{T}[D][B]t\mathrm{d}A\right] [K]=[∫Ae[B]T[D][B]tdA]
在对单元区域内进行积分时，可采用高斯—勒让德数值积分技术，至于为什么选用这个积分方法，我的猜想是与勒让德多项式积分区间在[-1,1]之内有关，个人猜想哈。
K = t ∫ − 1 + 1 ∫ − 1 + 1 [ B ( ξ , η ) ] T [ D ] [ B ( ξ , η ) ] det [ J ( ξ , η ) ] d η d ξ = t ∑ i = 1 n g p ∑ j = 1 n g p W i W j [ B ( ξ i , η j ] T [ D ] [ B ( ξ i , η j ) ] det [ J ( ξ i , η j ) ] \begin{aligned} K&=t\int_{-1}^{+1}{\int_{-1}^{+1}{\left[ B\left( \xi ,\eta \right) \right] ^{\mathrm{T}}[D][B(\xi ,\eta )]\det\mathrm{[}J(\xi ,\eta )]\mathrm{d}\eta \mathrm{d}\xi}}\\ &=t\sum_{i=1}^{\mathrm{ngp}}{\sum_{j=1}^{\mathrm{ngp}}{W_i}}W_j[B(\xi _i,\eta _j]^{\mathrm{T}}[D][B(\xi _i,\eta _j)]\det\mathrm{[}J(\xi _i,\eta _j)]\\ \end{aligned} K=t∫−1+1∫−1+1[B(ξ,η)]T[D][B(ξ,η)]det[J(ξ,η)]dηdξ=ti=1∑ngpj=1∑ngpWiWj[B(ξi,ηj]T[D][B(ξi,ηj)]det[J(ξi,ηj)]
上式中，t为单元厚度，ngp为高斯积分点的数量，W为高斯点的权重系数。常用高斯积分点位置及权重系数汇总见下表：
for n=1:size(gaussWeights,1)
% 高斯积分点位置
xi_Gauss=gaussLocations_cols(n,1);
eta_Gauss=gaussLocations_cols(n,2);
% 形函数在自然坐标系下的偏导
[shapeFunction,naturalDerivatives] = shapeFunctionQuad(xi_Gauss,eta_Gauss,elemType);
% 计算雅可比矩阵及形函数在物理坐标系下的偏导
[Jacob,XYderivatives]=Jacobian(elemNodeCoordinate,naturalDerivatives);
% 计算几何矩阵B
B(1,1:2:end) = XYderivatives(:,1)';
B(2,2:2:end) = XYderivatives(:,2)';
B(3,1:2:end) = XYderivatives(:,2)';
B(3,2:2:end) = XYderivatives(:,1)';
% 积分点循环计算单元刚度矩阵
k = k+B.'*D*t*B*gaussWeights(n)*det(Jacob);
end
绕节点平均应力技术和面力、体力施加与前一节的常应变三角形单元完全一样，可以照搬，这里就不详细展开了。
积分点应力外插至节点
可参考历史推文：有限元计算过程中积分点应力如何外插至节点处？【公式推导篇】和单元积分点应力如何外插至节点上 | 数值实现篇。
核心理念：坐标系的转换。
假设
(
ξ
,
η
)
\left( \xi ,\eta \right)
(ξ,η)是母单元的自然坐标系，
(
ξ
^
,
η
^
)
\left( \hat{\xi},\hat{\eta} \right)
(ξ^,η^)是由高斯积分点控制的坐标系（术语可能不专业），假设高斯积分方案为
2
×
2
2 \times2
2×2。坐标系转换关系：
(
ξ
,
η
)
=
(
ξ
^
,
η
^
)
/
3
o
r
(
ξ
^
,
η
^
)
=
(
ξ
,
η
)
3
(\xi,\eta)=(\hat{\xi},\hat{\eta})/\sqrt{3}\quad\mathrm{or}\quad(\hat{\xi},\hat{\eta})=(\xi,\eta)\sqrt{3}
(ξ,η)=(ξ^,η^)/3or(ξ^,η^)=(ξ,η)3
单元内任一点的应力 σ P \sigma _P σP，由4个高斯积分点应力 σ G i \sigma _{Gi} σGi进行插值时，可表示为
σ P = ∑ i = 1 4 N i σ G i = [ N 1 ( ξ ^ , η ^ ) N 2 ( ξ ^ , η ^ ) N 3 ( ξ ^ , η ^ ) N 4 ( ξ ^ , η ^ ) ] [ σ G 1 σ G 2 σ G 3 σ G 4 ] \sigma_{P}=\sum_{i=1}^{4}N_{i}\sigma_{Gi}=\begin{bmatrix}N_{1}(\hat{\xi},\hat{\eta})N_{2}(\hat{\xi},\hat{\eta})N_{3}(\hat{\xi},\hat{\eta})N_{4}(\hat{\xi},\hat{\eta})\end{bmatrix}\begin{bmatrix}\sigma_{G1}\\\sigma_{G2}\\\sigma_{G3}\\\sigma_{G4}\end{bmatrix} σP=i=1∑4NiσGi=[N1(ξ^,η^)N2(ξ^,η^)N3(ξ^,η^)N4(ξ^,η^)] σG1σG2σG3σG4
其中， N ( ξ ^ , η ^ ) N(\hat{\xi},\hat{\eta}) N(ξ^,η^)是基于高斯积分点的形函数，第一个积分点的坐标在母单元坐标系下为（-1，-1），根据上述的坐标系转换的方式，在高斯积分点的坐标系下，第一个单元节点在高斯积分点坐标系下坐标为 ( − 3 , − 3 ) \left( -\sqrt{3},-\sqrt{3} \right) (−3,−3)，将此坐标值代入第一个形函数，得 1 + 3 2 1+\frac{\sqrt{3}}{2} 1+23，相同的道理，可推导至四个节点在4个形函数下的 4 × 4 4 \times 4 4×4外插矩阵：
[ σ 1 σ 2 σ 3 σ 4 ] = [ 1 + 0.5 3 − 0.5 1 − 0.5 3 − 0.5 − 0.5 1 + 0.5 3 − 0.5 1 − 0.5 3 1 − 0.5 3 − 0.5 1 + 0.5 3 − 0.5 − 0.5 1 − 0.5 3 − 0.5 1 + 0.5 3 ] [ σ G 1 σ G 2 σ G 3 σ G 4 ] \begin{bmatrix}\sigma_1\\\sigma_2\\\sigma_3\\\sigma_4\end{bmatrix}=\begin{bmatrix}1+0.5\sqrt{3}&-0.5&1-0.5\sqrt{3}&-0.5\\-0.5&1+0.5\sqrt{3}&-0.5&1-0.5\sqrt{3}\\1-0.5\sqrt{3}&-0.5&1+0.5\sqrt{3}&-0.5\\-0.5&1-0.5\sqrt{3}&-0.5&1+0.5\sqrt{3}\end{bmatrix}\begin{bmatrix}\sigma_{G1}\\\sigma_{G2}\\\sigma_{G3}\\\sigma_{G4}\end{bmatrix} σ1σ2σ3σ4 = 1+0.53−0.51−0.53−0.5−0.51+0.53−0.51−0.531−0.53−0.51+0.53−0.5−0.51−0.53−0.51+0.53 σG1σG2σG3σG4
for i = 1:3
StressElem(e,:,i) = [1+0.5*sqrt(3) -0.5 1-0.5*sqrt(3) -0.5;
-0.5 1+0.5*sqrt(3) -0.5 1-0.5*sqrt(3);
1-0.5*sqrt(3) -0.5 1+0.5*sqrt(3) -0.5;
-0.5 1-0.5*sqrt(3) -0.5 1+0.5*sqrt(3)]*...
[stress(e,1,i);stress(e,2,i);stress(e,3,i);stress(e,4,i)];
end
数值案例
先考虑如下的单侧受拉模型，如下图所示，模型左侧节点的自由度被固定（UX=UY=0），右侧的节点承受大小为1向右的位移，材料属性见前处理文件：
*Material
1e+06,0.3,1.0,1
弹性模量，泊松比，厚度，平面应力/应变标识
场输出变量有：
- 节点位移：Umag、U1、U2
- 节点应力: Sx、Sy、Sxy
- 节点应力: Von Mises
精度验证
位移场对比
应力场 S x x S_{xx} Sxx对比
应力场 S y y S_{yy} Syy对比
应力场 S x y S_{xy} Sxy对比
文档精心排版的PDF缩略图
9719
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
