10. Finite element implementation#
The finite element representation of the matrices of [§9] is given below. These matrices appear in expressions to be integrated along the beam. They are therefore calculated at the Gauss points…. \({N}_{i},{N}_{j}\) are the values taken, at the Gauss point in question, by the shape functions relating to the nodes \(i,j\)…
Stiffness matrices relate the increase in displacement of node \(j\) to the increase in internal force at node \(i\) for element \(e\).
10.1. Deformation matrix and internal forces#
The deformation matrix has the continuous expression [éq 6.3-5]:
\(B=\left[\begin{array}{cc}\frac{d}{\text{ds}}1& {x}_{o}\text{'}\wedge \\ 0& \frac{d}{\text{ds}}1\end{array}\right]\).
In finite elements, the contribution of the displacement of node \(i\) to the deformation at the Gauss point in question is obtained by multiplying the 6 components of this displacement by the matrix:
\({B}_{i}^{h}=\left[\begin{array}{cc}{N}_{i}\text{'}1& {N}_{i}{\stackrel{ˆ}{x}}_{o}\text{'}(G)\\ 0& {N}_{i}\text{'}1\end{array}\right]\) eq 10.1-1
The higher index \(h\) indicates that it is a discretized form of the \(B\) matrix.
According to [éq 6.3-4]:
\({F}_{\text{int}i}^{e}={\int }_{e}{B}_{i}^{hT}\left\{\begin{array}{c}f\\ m\end{array}\right\}\text{ds}\) eq 10.1-2
is the interior effort applied to node \(i\) of element \(e\) and due to the generalized constraint field \(\left\{\begin{array}{c}f\\ m\end{array}\right\}\) in the element. This constraint is calculated according to the equation [éq 8.1-2]. But it should be noted that:
on the one hand [éq 6.3-3]: \({R}_{\text{tot}}^{T}\varepsilon ={R}_{\text{tot}}^{T}{\dot{x}}_{o}-{e}_{1};\)
on the other hand, vector :math:`chi` is**updated* at each iteration, as shown in [§10.8].
10.2. Stiffness matrices#
The continuous expression for the material stiffness matrix [éq 9.1-5] is:
\({B}^{T}{\Pi C\Pi }^{T}B\).
For the finite element \(e\), we deduce the matrix relating the increase in displacement of node \(j\) to the increase in internal force at node \(i\):
\({S}_{\text{mat}ij}^{e}={\int }_{e}{B}_{i}^{\text{hT}}{\Pi C\Pi }^{T}{B}_{j}^{h}\text{ds}\text{.}\) eq 10.2-1
We numerically calculate \({\Pi C\Pi }^{T}\) at the current Gauss point.
The continuous expression for the geometric rigidity matrix [éq 9.1-4] is:
\({Y}^{T}EY\)
where:
\(Y=\left[\begin{array}{cc}\frac{d}{\text{ds}}1& 0\\ 0& \frac{d}{\text{ds}}1\\ 0& 1\end{array}\right],\)
and \(E\) is given by the equation [éq 9.1-2].
As for the material stiffness matrix, we deduce:
\({S}_{\text{géom}ij}^{e}={\int }_{e}{Y}_{i}^{\text{hT}}{EY}_{j}^{h}\text{ds},\) eq 10.2-2
where:
\({Y}_{i}^{h}=\left[\begin{array}{cc}{N}_{i}^{\text{'}}1& 0\\ 0& {N}_{i}^{\text{'}}1\\ 0& {N}_{i}^{\text{'}}1\end{array}\right],\)
and where \(E\) is calculated numerically at the current Gauss point.
The continuous expression for the complementary stiffness matrix [éq 9.1-5] is:
\({B}^{T}Z,\)
where \(Z\) is given by the equation [éq 9.1-3].
From this we deduce:
\({S}_{\text{compl}\mathrm{ij}}^{e}={\int }_{e}{B}_{i}^{\text{hT}}{N}_{j}Z\text{ds}\),
where \(Z\) is calculated numerically at the current Gauss point.
10.3. Inertial forces#
According to [éq 7.3-1]:
\({F}_{\mathrm{iner},i}^{e}=-{\int }_{e}\left\{\begin{array}{}{N}_{i}\rho A{\ddot{x}}_{0}\\ {N}_{i}({I}_{\rho }\dot{\omega }+\omega \wedge {I}_{\rho }\omega )\end{array}\right\}\mathrm{ds}\) eq 10.3-1
is the inertial force of element \(e\) at node \(i\).
10.4. Inertia matrix#
Let’s say:
\({I}_{\mathrm{0,}\mathrm{iner}}=\left[\begin{array}{cc}\rho A1& 0\\ 0& {J}_{\rho }\end{array}\right]\) and \({M}_{0\mathrm{iner}\mathrm{ij}}^{e}={\int }_{e}{N}_{i}{N}_{j}{I}_{0\mathrm{iner}}\mathrm{ds}\)
We can see, from [éq 7-4] and [éq 9.2.1-1], that \(\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{iner}}^{e}\) is indeed the inertia matrix of the element \(e\) for the translational movement (diagonal sub-matrix \(\rho A1\) of \({I}_{0\mathrm{iner}}\)). But according to [éq9.2.2.1‑1] and [éq 9.2.2.2-1], it’s not the rotational inertia matrix.
However, matrix \({M}_{o\text{iner}}\), an assembly of \({M}_{o\text{iner}}^{e}\), is used to calculate the initial acceleration of the beam when it leaves its reference position with a zero initial speed \(({\omega }_{o}=0)\). Indeed, we then have, according to [éq 7.3-1]:
\({F}_{\text{ext}}(t=0)=-{F}_{\text{iner}}(t=0)={M}_{o\text{iner}}{\left\{\begin{array}{c}{\ddot{x}}_{o}\\ \dot{w}\end{array}\right\}}_{o},\)
since, in the reference position:
\({I}_{\rho }={J}_{\rho }\)
Let’s call \(J\text{'}\) the sum of the two matrices [éq 9.2.2.1-1] and [éq 9.2.2.2-1] and set:
\({I}_{\text{rot}}=\left[\begin{array}{cc}0& 0\\ 0& J\text{'}\end{array}\right]\) and \({M}_{\text{rot}ij}^{e}={\int }_{e}{N}_{i}{N}_{j}{I}_{\text{rot}}\text{ds}\)
Let’s also ask:
\({{I}_{0}}_{\text{rot}}=\left[\begin{array}{cc}0& 0\\ 0& {J}_{\rho }\end{array}\right]\text{et}{{M}_{0}}_{\text{rot}ij}^{e}={\int }_{e}{N}_{i}{N}_{j}{{I}_{0}}_{\text{rot}}\text{ds}\)
The complete inertia matrix of a beam element is obviously:
\({M}_{\text{iner}}^{e}=\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{iner}}^{e}-\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{rot}}^{e}+{M}_{\text{rot}}^{e}\).
10.5. External forces given#
According to [éq 8.1-3]:
\({F}_{\text{ext}i}^{e}={\int }_{e}\left\{\begin{array}{c}{N}_{i}\overline{f}\\ {N}_{i}\overline{m}\end{array}\right\}\text{ds}\) eq 10.5-1
is the force applied to node \(i\) of element \(e\) which is equivalent to the distributed external forces.
10.6. Linear iteration system#
By discretization in finite elements, equation [éq 9-1] gives:
\({W}^{h}=({F}_{\text{int}}-{F}_{\text{iner}}-{F}_{\text{ext}})\text{.}\delta \phi\).
On the other hand, assuming that the external forces are conservative, we have, according to [§10.2] and [§10.4]:
\({\mathrm{DW}}^{h}\text{.}\Delta \phi =\left[{S}_{\text{mat}}+{S}_{\text{geom}}+{S}_{\text{compl}}+{M}_{\text{rot}}-\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{rot}}+\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{iner}}\right]\Delta \phi \text{.}\delta \phi \text{.}\)
The relationship [éq 9-2], which must be verified by all \(\delta \phi\), therefore leads, at iteration \(n\) from time \(i\), to the following linear system in \(\Delta {\phi }_{i}^{n+1}\):
\(\begin{array}{}\left[{S}_{\text{mat,}i}^{n}+{S}_{\text{geom,}i}^{n}+{S}_{\text{compl,}i}^{n}+{M}_{\text{rot,}i}^{n}-\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{rot}}+\frac{1}{b\Delta {t}^{2}}{M}_{o\text{iner}}\right]\left\{\begin{array}{c}\Delta {x}_{o,i}^{n+1}\\ \Delta {\Psi }_{i}^{n+1}\end{array}\right\}\\ ={F}_{\text{ext,}i}-{F}_{\text{int,}i}^{n}+{F}_{\text{iner,}i}^{n}\text{.}\end{array}\) eq 10.6-1
In the*Code_Aster*, the elementary matrices \({M}_{o\text{iner}}^{e},\), which are independent of the displacement, are assembled only once to constitute the global matrix \({M}_{o\text{iner}},\) which is used in particular to calculate the initial acceleration [§10.4].
The three elementary stiffness matrices \({S}_{\text{mat}}^{e},{S}_{\text{geom}}^{e}\) and \({S}_{\text{compl}}^{e}\), the rotation inertia matrix, \({M}_{\text{rot}}^{e}\), which all four depend on the displacement, and the corrective rotational inertia matrix \(-\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{rot}}^{e},\), which is invariable, are, at each iteration, combined and then assembled to constitute a pseudo-global stiffness matrix \({\tilde{K}}_{i}^{n}\text{.}\)
So the linear system [éq 10.6-1] becomes:
\(\left[{\tilde{K}}_{i}^{n}+\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{iner}}\right]\left\{\begin{array}{c}\Delta {x}_{o,i}^{n+1}\\ \Delta {\Psi }_{i}^{n+1}\end{array}\right\}={F}_{\text{ext ,}i}-{F}_{\text{int,}i}^{n}+{F}_{\text{iner,}i}^{n}\text{.}\) eq 10.6-2
In the case of a static problem, the previous system is simplified to:
\(\left[{K}_{i}^{n}\right]\left\{\begin{array}{c}\Delta {x}_{o,i}^{n+1}\\ \Delta {\Psi }_{i}^{n+1}\end{array}\right\}={F}_{\text{ext,}i}-{F}_{\text{int,}i}^{n}\) eq 10.6-3
where \({K}_{i}^{n}\) is the assembly of the only stiffness matrices at iteration \(n\) of « instant » \(i\) [§8.2]:
\({S}_{\text{mat,}i}^{n}+{S}_{\text{geom,}i}^{n}+{S}_{\text{compl,}i}^{n}\text{.}\)
10.7. Displacement, speed, and acceleration update#
The treatment of translational motion is classical; that of rotational motion is done using quaternions [A8].
10.7.1. Translational movement#
The formulas [éq A3.1-3] and [éq A3.1-4] are applied.
10.7.2. Rotational movement#
The quantities to be updated are:
on the one hand, the rotation vector, the angular speed and the angular acceleration;
on the other hand, for subsequent calculations, the rotation matrix and the increment of the rotation-vector from time \(i-1\) to the current iteration from time \(i\) [§A6].
The update of the rotation vectors is based on the following property of quaternions [§A8.4]: « the quaternion of the product of two rotations is equal to the product of the quaternions of the component rotations ».
So by asking [§An8.6]:
\(\begin{array}{}(\Delta {x}_{o},\Delta x)=Q(\Delta {\Psi }_{i}^{n+1})\\ ({x}_{o},x)=Q({\Psi }_{i}^{n}),\end{array}\)
He comes:
\({\Psi }_{i}^{n+1}={Q}^{-1}\left[(\Delta {x}_{o},\Delta x)°({x}_{o},x)\right]\text{.}\) eq 10.7.2-1
On the other hand, according to equation [éq An6-2], if:
\(({y}_{o},y)=Q({\Psi }_{i-\mathrm{1,}i}^{n}),\)
so:
\({\Psi }_{i-\mathrm{1,}i}^{n+1}={Q}^{-1}\left[(\Delta {x}_{o},\Delta x)°({y}_{o},y)\right]\text{.}\) eq 10.7.2-2
The rotation matrix is updated immediately [§4.2]:
\({R}_{i}^{n+1}=\text{exp}({\stackrel{ˆ}{\Psi }}_{i}^{n+1}),\)
which is calculated according to [éq 4.2-1].
Finally, angular speed and acceleration are updated by the relationships [éq A3.2-3] and [éqA3.2‑4].
10.8. Curvature variation vector update#
The vector \(\chi\), which defines rotational deformation [§6.2], should only be calculated at Gauss points. In the*Code_Aster*, it is treated computationally as an « internal variable ».
According to [éq 6.2-2]:
\({\stackrel{ˆ}{\chi }}_{i}^{n}=\frac{d}{\text{ds}}({R}_{i}^{n})\text{.}{R}_{i}^{\mathrm{nT}}\text{.}\)
And, according to [éq 4.2-6]:
\({\stackrel{ˆ}{\chi }}_{i}^{n+1}=\frac{d}{\text{ds}}\left[\text{exp}(\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}){R}_{i}^{n}\right]{R}_{i}^{\mathrm{nT}}\text{exp}(-\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}})\text{.}\)
Either:
\(\begin{array}{}{\stackrel{ˆ}{\chi }}_{i}^{n+1}=\frac{d}{\text{ds}}\left[\text{exp}(\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}})\right]\text{exp}(-\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}})\\ +\text{exp}(\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}){\stackrel{ˆ}{\chi }}_{i}^{n}\text{exp}(-\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}})\text{.}\end{array}\)
In the previous equation, the matrix of the first member is anti-symmetric by construction; the second matrix of the second member is obviously antisymmetric; therefore the first matrix of the second member is also antisymmetric.
We show in [bib2] Appendix B that the axial vector of this last one matrix is:
\(\begin{array}{cc}\beta =& \frac{\text{sin}\parallel \Delta {\Psi }_{i}^{n+1}\parallel }{\parallel \Delta {\Psi }_{i}^{n+1}\parallel }(\Delta {\Psi }_{i}^{n+1})\text{'}+(1-\frac{\text{sin}\parallel \Delta {\Psi }_{i}^{n+1}\parallel }{\parallel \Delta {\Psi }_{i}^{n+1}\parallel })\frac{\Delta {\Psi }_{i}^{n+1}\text{.}(\Delta {\Psi }_{i}^{n+1})\text{'}}{\parallel \Delta {\Psi }_{i}^{n+1}\parallel }\frac{\Delta {\Psi }_{i}^{n+1}}{\parallel \Delta {\Psi }_{i}^{n+1}\parallel }\\ & +\frac{1}{2}{(\frac{\text{sin}\frac{1}{2}\parallel \Delta {\Psi }_{i}^{n+1}\parallel }{\frac{1}{2}\parallel \Delta {\Psi }_{i}^{n+1}\parallel })}^{2}\Delta {\Psi }_{i}^{n+1}\wedge (\Delta {\Psi }_{i}^{n+1})\text{'}\text{.}\end{array}\)
So:
\({\chi }_{i}^{n+1}=\beta +\text{AXIAL}\left[\text{exp}(\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}){\stackrel{ˆ}{\chi }}_{i}^{n}\text{exp}(-\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}})\right]\text{.}\)
10.9. Initialization before iterations#
In the dynamic case, if the load is constant over time, iterations can only start at time \(i\) if we initialize some of the displacement fields, speed and acceleration to values different from those of time \(i-1\). These initializations are done as follows.
10.9.1. Translational movement#
\({x}_{o,i}^{o}={x}_{o,i-1}\text{.}\)
Then, according to equation [éq An3.1-1]:
\({\ddot{x}}_{o,i}^{o}=-\frac{1}{\beta \Delta t}{\dot{x}}_{o,i-1}+\frac{2\beta -1}{2\beta }{\ddot{x}}_{o,i-1}\text{.}\)
From equation [éq An3.1-2]:
\({\dot{x}}_{o,i}^{o}={\dot{x}}_{o,i-1}+\Delta t\left[(1-\gamma ){\ddot{x}}_{o,i-1}+\gamma {\ddot{x}}_{o,i}^{o}\right]\text{.}\)
10.9.2. Rotational movement#
Expressions similar to the preceding ones are taken:
\({\Psi }_{i}^{o}={\Psi }_{i-1}\) eq 10.9.2-1
\({\dot{\omega }}_{i}^{o}=-\frac{1}{\beta \Delta t}{\omega }_{i-1}+\frac{2\beta -1}{2\beta }{\dot{\omega }}_{i-1}\) eq 10.9.2-2
\({\omega }_{i}^{o}={\omega }_{i-1}+\Delta t\left[(1-\gamma ){\dot{\omega }}_{i-1}+\gamma {\dot{\omega }}_{i}^{o}\right]\text{.}\) eq 10.9.2-3
The second members of [éq 10.9.2-2] and [éq 10.9.2-3] make sense because all the vectors in them are in the vector space tangent to \(\text{SO}(3)\) by \({R}_{i-1}\).
As consequences of equation [éq 10.9.2-1]:
\({R}_{i}^{o}={R}_{i-1}\)
and:
\({\Psi }_{i-\mathrm{1,}i}^{o}=0\text{.}\)
We can see that at the first moment \((i=1)\), the initializations require the knowledge of the initial acceleration \({\left\{\begin{array}{c}{\ddot{x}}_{o}\\ \dot{\omega }\end{array}\right\}}_{o},\) whose calculation is indicated in [§10.4].