1.1.1 Finite Element Discretization
The general energy equation below has already been derived by virtual work principle in accordance with Updated Lagrange formulation here. The next step is to transform it into finite element form for truss element. \({}_{t} \boldsymbol{e}_{ij}\) and \({}_{t} \boldsymbol{\eta}_{ij}\) are the linear and nonlinear decomposition of the Green-Lagrange strain tensor.
\(\displaystyle \begin{aligned}
&\int\limits_{{}^{t}V} {}_{t} C_{ijkl} \: {}_{t} e_{kl} \: \delta {}_{t} e_{ij} \: d^t V
+ \int\limits_{{}^{t}V} {}_{t} C_{ijkl} \: {}_{t} \eta_{kl} \: \delta {}_{t} e_{ij} \: d^t V
+ \int\limits_{{}^{t}V} {}_{t} C_{ijkl} \: {}_{t} e_{kl} \: \delta {}_{t} \eta_{ij} \: d^t V \\
&+ \int\limits_{{}^{t}V} {}_{t} C_{ijkl} \: {}_{t} \eta_{kl} \: \delta {}_{t} \eta_{ij} \: d^t V
+ \int\limits_{{}^{t}V} {}^{t} \sigma_{ij} \: \delta {}_{t} \eta_{ij} \: d^t V
+ \int\limits_{{}^{t}V} {}^{t} \sigma_{ij} \: \delta {}_{t} e_{ij} \: d^t V
= {}^{t+ \Delta t} R
\end{aligned}\)
Since there is no resistance against shear or moment in truss element, we consider only the axial strain effect. The Green-Lagrange Strain tensor and its decompositions become:
\(\displaystyle E_{11} = \frac{\partial u}{\partial x} + \frac{1}{2}\left[\left(\frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial v}{\partial x}\right)^2\right]\)
\(\displaystyle {e}_{11} = \frac{\partial u}{\partial x} \)
\(\displaystyle {\eta}_{11} = \frac{1}{2}\left[\left(\frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial v}{\partial x}\right)^2\right]\)
The first variations of linear and nonlinear strains:
\(\displaystyle \delta{e}_{11} = \frac{\partial \delta u}{\partial x} \)
\(\displaystyle \delta \eta_{11} = \frac{\partial u}{\partial x}\frac{\partial \delta u}{\partial x} + \frac{\partial v}{\partial x}\frac{\partial \delta v}{\partial x}\)
Let’s define elastic modulus for the axial direction as E. Since reference axis is regularly changed in each iteration, it is necessary to update elastic module with the equation given below (see derivation). Also, consider the cross-sectional area is constant during nonlinear process and called A.
\(\displaystyle {}_{t} C_{ijkl} = {}_{}^{t}E = {}_{}^{0}E\left(\frac{_{}^{t}L}{_{}^{0}L}\right)^3\)
\(\displaystyle dV = A\:dx\)
Since the area (A) is constant, from Cauchy’s stress tensor \(({}_{}^{t}\sigma_{ij})\):
\({}_{}^{t}F_1 = \int\limits_A{}_{}^{t}\sigma_{11}dA\:\:\rightarrow \:\: F_1 = \sigma_{11}A\)
It’s the total axial force along truss at time t or \((C_1)\) configuration.
For a truss element, we can define the displacement field by nodal displacements and shape functions:
\(\displaystyle u = \sum_iN_iu_i \) and \(\displaystyle v = \sum_iN_iv_i \)
\(\displaystyle \left\{\begin{array}{cc} N_1(x) \\ N_2(x) \\ \end{array}\right\} = \left\{\begin{array}{cc} 1 – \frac{x}{L} \\ \frac{x}{L} \\ \end{array}\right\}\)
\(\displaystyle \left\{\begin{array}{cc} \frac{dN_1}{dx} \\ \frac{dN_2}{dx} \\ \end{array}\right\} = \left\{\begin{array}{cc} – \frac{1}{L} \\ \frac{1}{L} \\ \end{array}\right\}\)
\(\displaystyle \frac{\partial u}{\partial x} = \sum_i\frac{dN_i}{dx}u_i = \left[\begin{array}{cc} \frac{dN_1}{dx} & 0 & \frac{dN_2}{dx} & 0 \\ \end{array}\right]\left\{\begin{array}{cc} u_A \\ v_A \\ u_B \\ v_B \\ \end{array}\right\} = \left[B_1\right] \left[u\right]\)
\(\displaystyle \frac{\partial v}{\partial x} = \sum_i\frac{dN_i}{dx}v_i = \left[\begin{array}{center center} 0 & \frac{dN_1}{dx} & 0 & \frac{dN_2}{dx} \\ \end{array}\right]\left\{\begin{array}{center center} u_A \\ v_A \\ u_B \\ v_B \\ \end{array}\right\} = \left[B_2\right]\left[u\right]\)
1.1.2 Stiffness Matrices
Let’s examine each term in the above equation by one by:
\(\begin{aligned}
\int\limits_{{}^{t}V} {}_{t} C_{ijkl} \ {}_{t} e_{kl} \ \delta {}_{t} e_{ij} \ d^tV
&= \int\limits_0^L E \left[B_1\right] \left[u\right] \left[B_1\right] \left[\delta u\right] A \ dx \\
&= u^TEA\int\limits_0^L\left[B_1\right]^T\left[B_1\right]dx \: \delta u\\
&=\delta u^T\frac{EA}{L}\left[\begin{array}{cc} 1 & 0 & – 1 & 0 \\ 0 & 0 & 0 & 0 \\ – 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ \end{array}\right]u\\
&=\delta u^T\left[k_e\right]u
\end{aligned}\)
\(\begin{aligned}
\int\limits_{{}^{t}V} {}_{t} C_{ijkl} \: {}_{t} \eta_{kl} \: \delta {}_{t} e_{ij} \: d^t V
&= \int\limits_0^LE\frac{1}{2}\left\{\left(\left[B_1\right]\left[u\right]\right)^2 + \left(\left[B_2\right]\left[u\right]\right)^2\right\}\left[B_1\right]\left[\delta u\right]A\:dx \\
&= u^T\frac{EA}{2}\int\limits_0^L\left(\left[B_1\right]^T\left[B_1\right] + \left[B_2\right]^T\left[B_2\right]\right)\left[u\right]\left[B_1\right]dx\:\delta u\\
&=\delta u^T\frac{EA}{2L^2}\left[\begin{array}{center center} \Delta u & \Delta v & – \Delta u & – \Delta v \\ 0 & 0 & 0 & 0 \\ – \Delta u & – \Delta v & \Delta u & \Delta v \\ 0 & 0 & 0 & 0 \\ \end{array}\right]u\\
&=\delta u^T\left[k_1\right]u
\end{aligned}\)
\(\begin{aligned}
\int\limits_{{}^{t}V} {}_{t} C_{ijkl} \: {}_{t} e_{kl} \: \delta {}_{t} \eta_{ij} \: d^t V
&= \int\limits_0^LE\left[B_1\right]\left[u\right]\left\{\left[B_1\right]\left[u\right]\left[B_1\right]\left[\delta u\right] + \left[B_2\right]\left[u\right]\left[B_2\right]\left[\delta u\right]\right\}A\:dx\\
&= u^TEA\int\limits_0^L\left[B_1\right]^T\left[u\right]^T\left(\left[B_1\right]^T\left[B_1\right] + \left[B_2\right]^T\left[B_2\right]\right)dx\:\delta u\\
&=\delta u^T\frac{EA}{2L^2}\left[\begin{array}{center center} 2\Delta u & 0 & – 2\Delta u & 0 \\ 2\Delta v & 0 & – 2\Delta u & 0 \\ – 2\Delta u & 0 & 2\Delta u & 0 \\ – 2\Delta v & 0 & 2\Delta v & 0 \\ \end{array}\right]u\\
&=\delta u^T\left[k_2\right]u
\end{aligned}\)
\(\begin{aligned}
\int\limits_{{}^{t}V} {}_{t} C_{ijkl} \: {}_{t} \eta_{kl} \: \delta {}_{t} \eta_{ij} \: d^t V
&= \int\limits_0^LE\frac{1}{2}\left\{\left(\left[B_1\right]\left[u\right]\right)^2 + \left(\left[B_2\right]\left[u\right]\right)^2\right\}\left\{\left[B_1\right]\left[u\right]\left[B_1\right]\left[\delta u\right] + \left[B_2\right]\left[u\right]\left[B_2\right]\left[\delta u\right]\right\}A\:dx\\
&= u^T\frac{EA}{2}\int\limits_0^L\left\{\left[B_1\right]^T\left[B_1\right] + \left[B_2\right]^T\left[B_2\right]\right\}\left[u\right]\left[u\right]^T\left\{\left[B_1\right]^T\left[B_1\right] + \left[B_2\right]^T\left[B_2\right]\right\}dx\:\delta u\\
&=\delta u^T\frac{EA}{2L^3}\left[\begin{array}{cc} \Delta u^2 & \Delta u\Delta v & – \Delta u^2 & – \Delta u\Delta v \\ \Delta u\Delta v & \Delta v^2 & – \Delta u\Delta v & – \Delta v^2 \\ – \Delta u^2 & – \Delta u\Delta v & \Delta u^2 & \Delta u\Delta v \\ – \Delta u\Delta v & – \Delta v^2 & \Delta u\Delta v & \Delta v^2 \\ \end{array}\right]u\\
&=\delta u^T\left[k_3\right]u
\end{aligned}\)
\(\begin{aligned}
\int\limits_{{}^{t}V} {}^{t} \sigma_{ij} \: \delta {}_{t} \eta_{ij} \: d^t V
&= \int\limits_0^L\sigma_{11}\left\{\left[B_1\right]\left[u\right]\left[B_1\right]\left[\delta u\right] + \left[B_2\right]\left[u\right]\left[B_2\right]\left[\delta u\right]\right\}A\:dx\\
&= u^TF_1\int\limits_0^L\left(\left[B_1\right]^T\left[B_1\right] + \left[B_2\right]^T\left[B_2\right]\right)dx\:\delta u\\
&=\delta u^T\frac{F_1}{L}\left[\begin{array}{cc} 1 & 0 & – 1 & 0 \\ 0 & 1 & 0 & – 1 \\ – 1 & 0 & 1 & 0 \\ 0 & – 1 & 0 & 1 \\ \end{array}\right]u\\
&=\delta u^T\left[k_g\right]u
\end{aligned}\)
\(\begin{aligned}
\int\limits_{{}^{t}V} {}^{t} \sigma_{ij} \: \delta {}_{t} e_{ij} \: d^t V
&= \int\limits_0^L\sigma_{11}\left[B_1\right]\left[\delta u\right]A\:dx\\
&= F_1\int\limits_0^L\left[B_1\right]dx\:\delta u\\
&=\delta u^T\left[\begin{array}{cc} – F_1 \\ 0 \\ F_1 \\ 0 \\ \end{array}\right]\\
&=\delta u^T\left[_{}^{1}f\right]
\end{aligned}\)
Similarly,
\(\begin{aligned}
{}^{t+ \Delta t} R
&= \int\limits_{{}^{t+ \Delta t}V} {}^{t+ \Delta t} \sigma_{ij} \: \delta {}_{t+ \Delta t} e_{ij} \: d{}^{t+ \Delta t}V\\
&= \delta u^T\left[_{}^{2}f\right]
\end{aligned}\)
\([k_e]\) and \([k_g]\) represent elastic and geometric stiffness matrix while \([k_1]\), \([k_2]\) and \([k_3]\) higher-order stiffness matrices. \({}^{1} f\) and \({}^{2} f\) are the internal force vectors at \((C_1)\) and \((C_2)\) configurations.
where
\(\Delta u = u_B \:-\: u_A\)
\(\Delta v = v_B \:-\: v_A\)
\(\left[_{}^1f\right]^T = \left[\begin{array}{cc} _{}^{1}F_{xa} & _{}^{1}F_{ya} & _{}^{1}F_{xb} & _{}^{1}F_{yb} \\ \end{array}\right]\)
\(\left[_{}^2f\right]^T = \left[\begin{array}{cc} _{}^{2}F_{xa} & _{}^{2}F_{ya} & _{}^{2}F_{xb} & _{}^{2}F_{yb} \\ \end{array}\right]\)
\(F_{xa}\) and \(F_{xb}\) have reversed signs while \(F_{ya}\) and \(F_{yb}\) are zero for both configurations since no shear force is resisted by truss element. Also, note that:
\(\displaystyle F_1={}_{}^1F_{xb}\)
1.1.3 Incremental Stiffness Equation
Substituting each term derived above into the virtual work equation yields:
\(\displaystyle \delta u^T\left(\left[k_e\right] + \left[k_g\right] + \left[k_1\right] + \left[k_2\right] + \left[k_3\right]\right)u + \delta u^T\left[_{}^{1}f\right] = \delta u^T\left[_{}^{2}f\right]\)
\(\boxed{\left(\left[k_e\right] + \left[k_g\right] + \left[k_1\right] + \left[k_2\right] + \left[k_3\right]\right)u + \left[_{}^{1}f\right] = \left[_{}^{2}f\right]}\)
1.1.4 Rotation Matrix
Let \(\theta \) is the positive angle it makes with the x-axis:
\(\displaystyle \left[R\right] = \left[\begin{array}{cc} \cos \theta & sin\theta & 0 & 0 \\ – sin\theta & \cos \theta & 0 & 0 \\ 0 & 0 & \cos \theta & sin\theta \\ 0 & 0 & – sin\theta & \cos \theta \\ \end{array}\right]\)
It’s used to transform local stiffness matrix \([k]\) to global stiffness matrix \([K]\).
\(\begin{aligned}
\left[K\right]
&= \left[R\right]^T\left(\left[k_e\right] + \left[k_g\right] + \left[k_1\right] + \left[k_2\right] + \left[k_3\right]\right)\left[R\right]\\
&= \left[R\right]^T\left[k\right]\left[R\right]
\end{aligned}\)
1.1.5 Geometry Update
After solving the incremental stiffness equation, global displacement matrices \((u_i^G)\) and \((v_i^G)\) for each member are obtained. These are also to be used for geometry update for the next increment in nonlinear algorithm.
\(\displaystyle {}^{t+ \Delta t} x_{i} = {}^{t} x_{i} + u_i^G\)
\(\displaystyle {}^{t+ \Delta t} y_{i} = {}^{t} y_{i} + v_i^G\)
For the new nodal coordinates, \(_{}^{2}L\) and \(_{}^{2}\theta\) is calculated for \((C_2)\) configuration by:
\(\displaystyle _{}^{2}L = \sqrt{({}^{t+ \Delta t} x_{2} \:-\: {}^{t+ \Delta t} x_{1})^2 + ({}^{t+ \Delta t} y_{2} \:-\: {}^{t+ \Delta t} y_{1})^2}\)
\(\displaystyle _{}^{2}\theta = arctan\left(\frac{{}^{t+ \Delta t} y_{2} \:-\: {}^{t+ \Delta t} y_{1}}{{}^{t+ \Delta t} x_{2} \:-\:{}^{t+ \Delta t} x_{1}}\right)\)
Be carefull that \(_{}^{2}\theta\) is acute angle, so it’s necessary to detect in which quadrant the truss is aligned on coordinate frame.
1.1.6 Force Recovery
The another role of global displacements \((u^G)\) is to calculate internal force along truss axis; however, we need to transform in into local displacements by rotation matrix \(R(_{}^{1}\theta)\). Here, the rotation matrix must be calculated with respect to the \((C_1)\) configuration because the local displacements are along the local axis of truss; therefore, the local force vector generated by this displacement vector is parallel to it like shown in the figure below.
\(\displaystyle [u^L] = \left[R(_{}^{1}\theta)\right][u^G]\)
By using \([u^L]\), let’s find the local internal forces generated by each stiffness matrix during \((C_1)\) to \((C_2)\):
\(\displaystyle [k_e][u^L] = \frac{EA}{L}\left[\begin{array}{cc} 1 & 0 & – 1 & 0 \\ 0 & 0 & 0 & 0 \\ – 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ \end{array}\right] \left\{\begin{array}{center center} u_A \\ v_A \\ u_B \\ v_B \\ \end{array}\right\} = \left[\begin{array}{cc} – F_{xl} \\ 0 \\ F_{xl} \\ 0 \\ \end{array}\right]\)
\(\displaystyle [k_g][u^L] = \frac{F_1}{L}\left[\begin{array}{cc} 1 & 0 & – 1 & 0 \\ 0 & 1 & 0 & -1 \\ – 1 & 0 & 1 & 0 \\ 0 & -1 & 0 & 1 \\ \end{array}\right] \left\{\begin{array}{center center} u_A \\ v_A \\ u_B \\ v_B \\ \end{array}\right\} = \left[\begin{array}{center center} – F_1\frac{\Delta u}{L} \\ – F_1\frac{\Delta v}{L} \\ F_1\frac{\Delta u}{L} \\ F_1\frac{\Delta v}{L} \\ \end{array}\right]\)
\(\displaystyle [k_1][u^L] = \frac{EA}{2L^2}\left[\begin{array}{cc} \Delta u & \Delta v & – \Delta u & – \Delta v \\ 0 & 0 & 0 & 0 \\ – \Delta u & – \Delta v & \Delta u & \Delta v \\ 0 & 0 & 0 & 0 \\ \end{array}\right] \left\{\begin{array}{center center} u_A \\ v_A \\ u_B \\ v_B \\ \end{array}\right\} = \left[\begin{array}{cc} – F_{xn} \\ 0 \\ F_{xn} \\ 0 \\ \end{array}\right]\)
\(\displaystyle [k_2][u^L] = \frac{EA}{2L^2}\left[\begin{array}{cc} 2\Delta u & 0 & – 2\Delta u & 0 \\ 2\Delta v & 0 & – 2\Delta u & 0 \\ – 2\Delta u & 0 & 2\Delta u & 0 \\ – 2\Delta v & 0 & 2\Delta v & 0 \\ \end{array}\right] \left\{\begin{array}{center center} u_A \\ v_A \\ u_B \\ v_B \\ \end{array}\right\} = \left[\begin{array}{center center} – F_{xl}\frac{\Delta u}{L} \\ – F_{xl}\frac{\Delta v}{L} \\ F_{xl}\frac{\Delta u}{L} \\ F_{xl}\frac{\Delta v}{L} \\ \end{array}\right]\)
\(\displaystyle [k_3][u^L] = \frac{EA}{2L^3}\left[\begin{array}{cc} \Delta u^2 & \Delta u\Delta v & – \Delta u^2 & – \Delta u\Delta v \\ \Delta u\Delta v & \Delta v^2 & – \Delta u\Delta v & – \Delta v^2 \\ – \Delta u^2 & – \Delta u\Delta v & \Delta u^2 & \Delta u\Delta v \\ – \Delta u\Delta v & – \Delta v^2 & \Delta u\Delta v & \Delta v^2 \\ \end{array}\right] \left\{\begin{array}{center center} u_A \\ v_A \\ u_B \\ v_B \\ \end{array}\right\} = \left[\begin{array}{center center} – F_{xn}\frac{\Delta u}{L} \\ – F_{xn}\frac{\Delta v}{L} \\ F_{xn}\frac{\Delta u}{L} \\ F_{xn}\frac{\Delta v}{L} \\ \end{array}\right]\)
where
\(\displaystyle F_{xl} = EA\frac{\Delta u}{L}\)
\(\displaystyle F_{xn} = \frac{EA}{2}\left(\frac{\Delta u^2}{L^2} + \frac{\Delta v^2}{L^2}\right)\)
The forces generated by \([k_2] \) are in effect to transform force increment \(F_{xl}\) generated by \([k_e] \):
\(\displaystyle \left[\begin{array}{cc} – F_{xl}\left(1 + \frac{\Delta u}{L}\right) & – F_{xl}\frac{\Delta v}{L} & F_{xl}\left(1 + \frac{\Delta u}{L}\right) & F_{xl}\frac{\Delta v}{L} \\ \end{array}\right]^T\)
The forces generated by \([k_3] \) are in effect to transform force increment \(F_{xn}\) generated by \([k_1] \):
\(\displaystyle \left[\begin{array}{cc} – F_{xn}\left(1 + \frac{\Delta u}{L}\right) & – F_{xn}\frac{\Delta v}{L} & F_{xn}\left(1 + \frac{\Delta u}{L}\right) & F_{xn}\frac{\Delta v}{L} \\ \end{array}\right]^T\)
The forces generated by \([k_g] \) are in effect to transform force increment \(F_1\) generated by \(\left[_{}^{1}f\right] \):
\(\displaystyle \left[\begin{array}{cc} – F_1\left(1 + \frac{\Delta u}{L}\right) & – F_1\frac{\Delta v}{L} & F_1\left(1 + \frac{\Delta u}{L}\right) & F_1\frac{\Delta v}{L} \\ \end{array}\right]^T\)
These forces are considered as the local forces shown with blue arrows in the figure above, so the summation of resultant forces for 3 vectors gives the axial force at \((C_2)\) configuration.
\(\begin{aligned}
_{2}^{}F_{xl}
&= \sqrt{\left[F_{xl}\left(1 + \frac{\Delta u}{L}\right)\right]^2 + \left[F_{xl}\frac{\Delta v}{L}\right]^2}\\
&=\frac{F_{xl}}{L}\sqrt{\left(\Delta u + L\right)^2 + \Delta v^2}\\
&= F_{xl}\frac{_{}^{2}L}{L}\\
&= EA\frac{\Delta u}{L}\frac{\left(_{}^{2}L\right)}{L}
\end{aligned}\)
\(\begin{aligned}
_{2}^{}F_{xn}
&=\frac{F_{xn}}{L}\sqrt{\left(\Delta u + L\right)^2 + \Delta v^2}\\
&= F_{xn}\frac{_{}^{2}L}{L}\\
&= \frac{EA}{2}\left(\frac{\Delta u^2}{L^2} + \frac{\Delta v^2}{L^2}\right)\frac{\left(_{}^{2}L\right)}{L}
\end{aligned}\)
\(\begin{aligned}
_{2}^{}F_1
&=\frac{F_1}{L}\sqrt{\left(\Delta u + L\right)^2 + \Delta v^2}\\
&= F_1\frac{_{}^{2}L}{L}
\end{aligned}\)
\(\boxed{{}_{2}^{}F_{xl} + {}_{2}^{}F_{xn} + {}_{2}^{}F_1 = \left(EA\left(\frac{\Delta u}{L} + \frac{1}{2}\left\{\left(\frac{\Delta u}{L}\right)^2 + \left(\frac{\Delta v}{L}\right)^2\right\}\right) + F_1\right)\frac{_{}^{2}L}{L} = {}_{}^{2}F_{xb}}\)
\(\displaystyle {}_{}^{2}F_{xb}\) is nothing but the total axial force at \((C_2)\) configuration. Since there is no shear, we can place it in \(\left[_{}^2f\right]^T\):
\(\left[_{}^2f\right]^T = \left[\begin{array}{cc} -_{}^{2}F_{xb} & 0 & _{}^{2}F_{xb} & 0 \\ \end{array}\right]\)
The final procedure is to find total global resisting forces \(({}_{}^2F_r)\) at \((C_2)\). Be carefull that, \({}_{}^2f\) is now at \((C_2)\) configuration; therefore, new angle \(_{}^{2}\theta\) must be used.
\(\boxed{{}_{}^2F_r = \left[R(_{}^2\theta )\right]^T \left[_{}^2f\right]}\)
References:
[1] Leu, L.-J., & Yang, Y.-B. (1990). Effects of rigid body and stretching on nonlinear analysis of trusses. Journal of Structural Engineering, 116(10), 2582-2598.
[2] Torkamani, M. A. M., & Shieh, J.-H. (2011). Higher-order stiffness matrices in nonlinear finite element analysis of plane truss structures. Engineering Structures, 33(12), 3516-3526.
[3] Yang YB, Kuo SR. (1994) Theory and analysis of nonlinear framed structures.