有限元分析
单元分析
根据该论文的表述,我们将根据第一部分的推导,构建算法的迭代过程。
可以意识到, L , V , H L,V,H L,V,H实际上是 F 1 , F 2 F_1,F_2 F1,F2的函数,因此可以对 ( 90 ) ∼ ( 92 ) (90)\sim (92) (90)∼(92)进行改写:
Because , F 3 = − F 1 F 4 = − F 2 + w L u T I = F 1 2 + F 2 2 T J = F 3 2 + F 4 2 = F 1 2 + ( − F 2 + w L u ) 2 So that , H = − F 1 ( L u o E A + 1 w ln ( − F 2 + w L u + F 1 2 + ( − F 2 + w L u ) 2 F 1 2 + F 2 2 − F 2 ) ) , V = 1 2 E A w ( ( − F 2 + w L u ) 2 − F 2 2 ) + F 1 2 + ( − F 2 + w L u ) 2 − F 1 2 + F 2 2 w , L = L u o + 1 2 E A w ( ( − F 2 + w L u ) F 1 2 + ( − F 2 + w L u ) 2 + F 2 F 1 2 + F 2 2 + F 1 2 ln ( − F 2 + w L u + F 1 2 + ( − F 2 + w L u ) 2 F 1 2 + F 2 2 − F 2 ) ) . \begin{align} &\textbf{Because},\notag\\ F_3=&-F_1\\ F_4=&-F_2+wL_u\\ T_I=&\sqrt{F_1^2+F_2^2}\\ T_J=&\sqrt{F_3^2+F_4^2}=\sqrt{F_1^2+(-F_2+wL_u)^2}\\ &\textbf{So that},\notag\\ H=&-F_1\left( \frac{L_{uo}}{EA}+\frac{1}{w}\ln \left(\frac{-F_2+wL_u+\sqrt{F_1^2+(-F_2+wL_u)^2}}{\sqrt{F_1^2+F_2^2}-F_2}\right) \right),\\ V=&\frac{1}{2EAw}\left((-F_2+wL_u)^2-F_2^2\right)+\frac{\sqrt{F_1^2+(-F_2+wL_u)^2}-\sqrt{F_1^2+F_2^2}}{w},\\ L=&L_{uo}+\frac{1}{2EAw} \left( (-F_2+wL_u)\sqrt{F_1^2+(-F_2+wL_u)^2}+F_2\sqrt{F_1^2+F_2^2}+\right.\notag\\ &\left.F_1^2\ln\left(\frac{-F_2+wL_u+\sqrt{F_1^2+(-F_2+wL_u)^2}}{\sqrt{F_1^2+F_2^2}-F_2}\right) \right). \end{align} F3=F4=TI=TJ=H=V=L=Because,−F1−F2+wLuF12+F22F32+F42=F12+(−F2+wLu)2So that,−F1(EALuo+w1ln(F12+F22−F2−F2+wLu+F12+(−F2+wLu)2)),2EAw1((−F2+wLu)2−F22)+wF12+(−F2+wLu)2−F12+F22,Luo+2EAw1((−F2+wLu)F12+(−F2+wLu)2+F2F12+F22+F12ln(F12+F22−F2−F2+wLu+F12+(−F2+wLu)2)).
接下来,考虑如图3的迭代过程。
其中 F 1 i F_1^i F1i和 F 2 i F_2^i F2i是迭代过程中第 i i i步的原始力; H i H^i Hi和 V i V^i Vi是由 ( 97 ) , ( 98 ) (97),(98) (97),(98)确定的线段投影;
Δ H i \Delta H^i ΔHi和 Δ V i \Delta V^i ΔVi是闭合差向量 J J ′ JJ^\prime JJ′(misclosure vector)的水平分量和垂直分量。
如果闭合差向量 J J ′ JJ^\prime JJ′超过一定容差值,则对 F 1 i F_1^i F1i和 F 2 i F_2^i F2i进行线性修正,得到步骤 i + 1 i + 1 i+1处的力。
( F 1 i + 1 F 2 i + 1 ) = ( F 1 i F 2 i ) + ( Δ F 1 i Δ F 2 i ) = ( F 1 i F 2 i ) + ( ∂ F 1 i ∂ H i ∂ F 1 i ∂ V i ∂ F 2 i ∂ H i ∂ F 2 i ∂ V i ) ( Δ H i Δ V i ) = k i ( Δ H i Δ V i ) . \begin{align} \begin{pmatrix} F_1^{i+1}\\ F_2^{i+1} \end{pmatrix} = \begin{pmatrix} F_1^{i}\\ F_2^{i} \end{pmatrix} + \begin{pmatrix} \Delta F_1^{i}\\ \Delta F_2^{i} \end{pmatrix} = \begin{pmatrix} F_1^{i}\\ F_2^{i} \end{pmatrix} + \begin{pmatrix} \dfrac{\partial F_1^i}{\partial H^i}&\dfrac{\partial F_1^i}{\partial V^i}\\[8pt] \dfrac{\partial F_2^i}{\partial H^i}&\dfrac{\partial F_2^i}{\partial V^i} \end{pmatrix} \begin{pmatrix} \Delta H^{i}\\ \Delta V^{i} \end{pmatrix} =k^i \begin{pmatrix} \Delta H^{i}\\ \Delta V^{i} \end{pmatrix}. \end{align} (F1i+1F2i+1)=(F1iF2i)+(ΔF1iΔF2i)=(F1iF2i)+ ∂Hi∂F1i∂Hi∂F2i∂Vi∂F1i∂Vi∂F2i (ΔHiΔVi)=ki(ΔHiΔVi).
k i k^i ki为第 i i i步的刚度矩阵,为求解 Δ H i \Delta H^i ΔHi和 Δ V i \Delta V^i ΔVi方便,可以将其用柔度矩阵 f i f^i fi表示为以下形式:
( Δ H i Δ V i ) = ( ∂ F 1 i ∂ H i ∂ F 1 i ∂ V i ∂ F 2 i ∂ H i ∂ F 2 i ∂ V i ) − 1 ( Δ H i Δ V i ) = ( ∂ H i ∂ F 1 i ∂ H i ∂ F 2 i ∂ V i ∂ F 1 i ∂ V i ∂ F 2 i ) ( Δ F 1 i Δ F 2 i ) = f i ( Δ F 1 i Δ F 2 i ) . \begin{align} \begin{pmatrix} \Delta H^{i}\\ \Delta V^{i} \end{pmatrix} = \begin{pmatrix} \dfrac{\partial F_1^i}{\partial H^i}&\dfrac{\partial F_1^i}{\partial V^i}\\[8pt] \dfrac{\partial F_2^i}{\partial H^i}&\dfrac{\partial F_2^i}{\partial V^i} \end{pmatrix}^{-1} \begin{pmatrix} \Delta H^{i}\\ \Delta V^{i} \end{pmatrix} = \begin{pmatrix} \dfrac{\partial H^i}{\partial F_1^i}&\dfrac{\partial H^i}{\partial F_2^i}\\[8pt] \dfrac{\partial V^i}{\partial F_1^i}&\dfrac{\partial V^i}{\partial F_2^i} \end{pmatrix} \begin{pmatrix} \Delta F_1^{i}\\ \Delta F_2^{i} \end{pmatrix} =f^i \begin{pmatrix} \Delta F_1^{i}\\ \Delta F_2^{i} \end{pmatrix}. \end{align} (ΔHiΔVi)= ∂Hi∂F1i∂Hi∂F2i∂Vi∂F1i∂Vi∂F2i −1(ΔHiΔVi)= ∂F1i∂Hi∂F1i∂Vi∂F2i∂Hi∂F2i∂Vi (ΔF1iΔF2i)=fi(ΔF1iΔF2i).
其中:
∂ H i ∂ F 1 i = − ( L u o E A + 1 w ln ( − F 2 i + w L u + F 1 i 2 + ( − F 2 i + w L u ) 2 F 1 i 2 + F 2 i 2 − F 2 i ) ) − F 1 i ( 1 w F 1 i 2 + F 2 i 2 − F 2 i − F 2 i + w L u + F 1 i 2 + ( − F 2 i + w L u ) 2 × 2 F 1 i F 1 i 2 + F 2 i 2 − F 2 i F 1 i 2 + ( − F 2 i + w L u ) 2 − 2 F 1 i − F 2 i + w L u + F 1 i 2 + ( − F 2 i + w L u ) 2 F 1 i 2 + F 2 i 2 − F 2 i ( F 1 i 2 + F 2 i 2 − F 2 i ) 2 ) , = H i F 1 i + 1 w ( F 4 i T J i + F 2 i T I i ) , ∂ H i ∂ F 2 i = F 1 i w ( 1 T J i − 1 T I i ) , ∂ V i ∂ F 1 i = F 1 i w ( 1 T J i − 1 T I i ) , ∂ V i ∂ F 2 i = − L i E A − 1 w ( F 4 i T J i + F 2 i T I i ) . \begin{align} \frac{\partial H^i}{\partial F_1^i}=&-\left( \frac{L_{uo}}{EA}+\frac{1}{w}\ln \left(\frac{-F_2^i+wL_u+\sqrt{F_1^{i2}+(-F_2^i+wL_u)^2}}{\sqrt{F_1^{i2}+F_2^{i2}}-F_2^{i}}\right) \right)\\ &-F_1^i\left( \frac{1}{w} \frac{\sqrt{F_1^{i2}+F_2^{i2}}-F_2^{i}}{-F_2^i+wL_u+\sqrt{F_1^{i2}+(-F_2^i+wL_u)^2}} \right.\notag\\ &\left.\times \frac{ 2F_1^i \dfrac{\sqrt{F_1^{i2}+F_2^{i2}}-F_2^{i}}{\sqrt{F_1^{i2}+(-F_2^i+wL_u)^2}} - 2F_1^i \dfrac{-F_2^i+wL_u+\sqrt{F_1^{i2}+(-F_2^i+wL_u)^2}}{\sqrt{F_1^{i2}+F_2^{i2}}-F_2^{i}} }{\left({\sqrt{F_1^{i2}+F_2^{i2}}-F_2^{i}}\right)^2} \right),\notag\\ =&\frac{H^i}{F_1^i}+\frac{1}{w}\left( \frac{F_4^i}{T_J^i}+\frac{F_2^i}{T_I^i} \right),\\ \frac{\partial H^i}{\partial F_2^i}=&\frac{F_1^i}{w}\left(\frac{1}{T_J^i}-\frac{1}{T_I^i}\right),\\ \frac{\partial V^i}{\partial F_1^i}=&\frac{F_1^i}{w}\left(\frac{1}{T_J^i}-\frac{1}{T_I^i}\right),\\ \frac{\partial V^i}{\partial F_2^i}=&-\frac{L^i}{EA}-\frac{1}{w}\left(\frac{F_4^i}{T_J^i}+\frac{F_2^i}{T_I^i}\right). \end{align} ∂F1i∂Hi==∂F2i∂Hi=∂F1i∂Vi=∂F2i∂Vi=−(EALuo+w1ln(F1i2+F2i2−F2i−F2i+wLu+F1i2+(−F2i+wLu)2))−F1i(w1−F2i+wLu+F1i2+(−F2i+wLu)2F1i2+F2i2−F2i×(F1i2+F2i2−F2i)22F1iF1i2+(−F2i+wLu)2F1i2+F2i2−F2i−2F1iF1i2+F2i2−F2i−F2i+wLu+F1i2+(−F2i+wLu)2 ,F1iHi+w1(TJiF4i+TIiF2i),wF1i(TJi1−TIi1),wF1i(TJi1−TIi1),−EALi−w1(TJiF4i+TIiF2i).
在以上所描述的迭代过程中需要初值 F 1 0 , F 2 0 F_1^0,F_2^0 F10,F20。
接下来将求解 F 1 0 F_1^0 F10,将方程 ( 88 ) (88) (88)进行变形,有:
F
1
0
=
−
w
H
2
λ
0
.
\begin{align} F_1^0=-\frac{wH}{2\lambda^0}. \end{align}
F10=−2λ0wH.
其中
λ
0
\lambda^0
λ0可由
(
87
)
(87)
(87)变换得到,将
sinh
2
λ
λ
2
\frac{\sinh^2 \lambda}{\lambda^2}
λ2sinh2λ进行级数展开:
sinh
2
λ
λ
2
=
1
+
1
3
λ
2
+
o
(
λ
4
)
.
\begin{align} \frac{\sinh^2 \lambda}{\lambda^2}=1+\frac{1}{3}\lambda^2+o(\lambda^4). \end{align}
λ2sinh2λ=1+31λ2+o(λ4).
然后代入到 ( 87 ) (87) (87)式中,即可推得:
λ 0 = 3 ( L u 2 − V 2 H 2 − 1 ) . \begin{align} \lambda^0=\sqrt{3\left(\frac{L_u^2-V^2}{H^2}-1\right)}. \end{align} λ0=3(H2Lu2−V2−1).
文献中还考虑到以下情况。
当对于一根非常紧的索,未拉伸长度 L u L_u Lu小于距离 ∣ I J ∣ \left|IJ\right| ∣IJ∣,无法解出 λ 0 \lambda^0 λ0时,应假设 λ 0 \lambda^0 λ0在垂跨比为 5 % 5\% 5%下的保守值,从而得出 λ 0 \lambda^0 λ0的值为 0.2 0.2 0.2。
其次,当 H H H接近 0 0 0时,应有 λ 0 > 1 0 6 \lambda^0>10^6 λ0>106。
现在求解 F 2 0 F_2^0 F20。为了保证在所有情况下的收敛性,任何点 J J J的纵坐标都应小于点 I I I的纵坐标(即 V < 0 V<0 V<0),所以在程序中算法应对坐标进行修正,在输出结果中进行还原。此处,利用式 ( 89 ) (89) (89)得到 F 2 F_2 F2的严格正起始值。
F 2 0 = w 2 ( − V cosh ( λ ) sinh ( λ ) + L u ) . \begin{align} F_2^0=\frac{w}{2}\left( -V\frac{\cosh\left(\lambda \right)}{\sinh\left(\lambda \right)}+L_u \right). \end{align} F20=2w(−Vsinh(λ)cosh(λ)+Lu).
利用 ( 93 ) ∼ ( 96 ) (93)\sim (96) (93)∼(96),可以得到 F 3 , F 4 , T I , T J , L F_3,F_4,T_I,T_J,L F3,F4,TI,TJ,L的值;用 L u L_u Lu的任意分数替换 L u L_u Lu,用 ( 90 ) ∼ ( 92 ) (90)\sim (92) (90)∼(92)可以确定沿电缆任意个数点的坐标,悬链的几何形状可以存储在阵列 { X C O O R } \{XCOOR\} {XCOOR}中。
至此,我们得到了悬链单元所有信息,可以编写单元的有限元程序。
整体分析
虽然悬链是非线性系统(几何非线性),但是在单步迭代过程中,可以简化为桁架模型。单元的刚度矩阵与端部位移以及由于端部位移引起的端部力的微小变化( α = 1 0 − 3 \alpha =10^{-3} α=10−3)有关。
{ F N } T = ( F I x 0 F I z F J x 0 F J z ) T = ( F 1 0 F 2 F 3 0 F 4 ) T { X N } T = ( X I x 0 X I z X J x 0 X J z ) T { Δ F N } = { K C I } { Δ X N } { K C N } = ( K C I I K C I J K C J I K C J J ) \begin{align} &\{F_N\}^T= \begin{pmatrix} F_{Ix}&0&F_{Iz}&F_{Jx}&0&F_{Jz} \end{pmatrix}^T = \begin{pmatrix} F_1&0&F_2&F_3&0&F_4 \end{pmatrix}^T\\ &\{X_N\}^T= \begin{pmatrix} X_{Ix}&0&X_{Iz}&X_{Jx}&0&X_{Jz} \end{pmatrix}^T\\ &\{\Delta F_N\}=\{KC_I\}\{\Delta X_N\}\\ &\{KC_N\} = \begin{pmatrix} KC_{II}&KC_{IJ}\\ KC_{JI}&KC_{JJ} \end{pmatrix} \end{align} {FN}T=(FIx0FIzFJx0FJz)T=(F10F2F30F4)T{XN}T=(XIx0XIzXJx0XJz)T{ΔFN}={KCI}{ΔXN}{KCN}=(KCIIKCJIKCIJKCJJ)
设$\alpha
是长度
是长度
是长度H
或
或
或V$变化的分数,可以认为是很小的,可以设其为
1
0
(
−
4
)
10^{(-4)}
10(−4)。在前面我们已经算出了索端的支反力
F
3
,
F
4
F_3,F_4
F3,F4,
我们可以令
F
3
′
,
F
4
′
F_3^\prime,F_4^\prime
F3′,F4′为
H
(
1
+
α
)
,
V
H(1+\alpha ),V
H(1+α),V的计算值,令
F
3
′
′
,
F
4
′
′
F_3^{\prime\prime},F_4^{\prime\prime}
F3′′,F4′′为
H
,
V
(
1
+
α
)
H,V(1+\alpha )
H,V(1+α)的计算值。
{
K
C
N
}
=
(
F
3
′
−
F
3
α
H
0
F
3
′
′
−
F
3
α
V
−
F
3
′
−
F
3
α
H
0
−
F
3
′
′
−
F
3
α
V
0
0
0
0
0
0
F
4
′
−
F
4
α
H
0
F
4
′
′
−
F
4
α
V
−
F
4
′
−
F
4
α
H
0
−
F
4
′
′
−
F
4
α
V
−
F
3
′
−
F
3
α
H
0
−
F
3
′
′
−
F
3
α
V
F
3
′
−
F
3
α
H
0
F
3
′
′
−
F
3
α
V
0
0
0
0
0
0
−
F
4
′
−
F
4
α
H
0
−
F
4
′
′
−
F
4
α
V
F
4
′
−
F
4
α
H
0
F
4
′
′
−
F
4
α
V
)
\begin{align} \{KC_N\} = \begin{pmatrix} \dfrac{F_3^{\prime}-F_3}{\alpha H}&0&\dfrac{F_3^{\prime\prime}-F_3}{\alpha V}&-\dfrac{F_3^{\prime}-F_3}{\alpha H}&0&-\dfrac{F_3^{\prime\prime}-F_3}{\alpha V}\\[8pt] 0&0&0&0&0&0\\[8pt] \dfrac{F_4^{\prime}-F_4}{\alpha H}&0&\dfrac{F_4^{\prime\prime}-F_4}{\alpha V}&-\dfrac{F_4^{\prime}-F_4}{\alpha H}&0&-\dfrac{F_4^{\prime\prime}-F_4}{\alpha V}\\[8pt] -\dfrac{F_3^{\prime}-F_3}{\alpha H}&0&-\dfrac{F_3^{\prime\prime}-F_3}{\alpha V}&\dfrac{F_3^{\prime}-F_3}{\alpha H}&0&\dfrac{F_3^{\prime\prime}-F_3}{\alpha V}\\[8pt] 0&0&0&0&0&0\\[8pt] -\dfrac{F_4^{\prime}-F_4}{\alpha H}&0&-\dfrac{F_4^{\prime\prime}-F_4}{\alpha V}&\dfrac{F_4^{\prime}-F_4}{\alpha H}&0&\dfrac{F_4^{\prime\prime}-F_4}{\alpha V}\\[8pt] \end{pmatrix} \end{align}
{KCN}=
αHF3′−F30αHF4′−F4−αHF3′−F30−αHF4′−F4000000αVF3′′−F30αVF4′′−F4−αVF3′′−F30−αVF4′′−F4−αHF3′−F30−αHF4′−F4αHF3′−F30αHF4′−F4000000−αVF3′′−F30−αVF4′′−F4αVF3′′−F30αVF4′′−F4
为了将单元刚度矩阵
{
K
C
N
}
\{KC_N\}
{KCN}整合到整体矩阵中,需要进行坐标变换,转换为整体坐标系中的刚度矩阵
{
K
N
}
\{K_N\}
{KN}。
{
K
N
}
=
{
T
N
}
{
K
C
N
}
{
T
N
}
T
\begin{align} \{K_N\}=\{T_N\}\{KC_N\}\{T_N\}^T \end{align}
{KN}={TN}{KCN}{TN}T
其中,
{
T
N
}
=
(
cos
α
sin
α
0
0
0
0
−
sin
α
cos
α
0
0
0
0
0
0
1
0
0
0
0
0
0
cos
α
sin
α
0
0
0
0
−
sin
α
cos
α
0
0
0
0
0
0
1
)
\begin{align} \{T_N\}= \begin{pmatrix} \cos \alpha &\sin\alpha &0&0&0&0\\ -\sin\alpha &\cos \alpha&0&0&0&0\\ 0&0&1&0&0&0\\ 0&0&0&\cos \alpha &\sin\alpha &0\\ 0&0&0&-\sin\alpha &\cos \alpha&0\\ 0&0&0&0&0&1\\ \end{pmatrix} \end{align}
{TN}=
cosα−sinα0000sinαcosα0000001000000cosα−sinα0000sinαcosα0000001