Finite element implementation ================================ The finite element representation of the matrices of [:ref:`§9 <§9>`] is given below. These matrices appear in expressions to be integrated along the beam. They are therefore calculated at the Gauss points.... :math:`{N}_{i},{N}_{j}` are the values taken, at the Gauss point in question, by the shape functions relating to the nodes :math:`i,j`... Stiffness matrices relate the increase in displacement of node :math:`j` to the increase in internal force at node :math:`i` for element :math:`e`. Deformation matrix and internal forces -------------------------------------------- The deformation matrix has the continuous expression [:ref:`éq 6.3-5 <éq 6.3-5>`]: :math:`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 :math:`i` to the deformation at the Gauss point in question is obtained by multiplying the 6 components of this displacement by the matrix: :math:`{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 :math:`h` indicates that it is a discretized form of the :math:`B` matrix. According to [:ref:`éq 6.3-4 <éq 6.3-4>`]: :math:`{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 :math:`i` of element :math:`e` and due to the generalized constraint field :math:`\left\{\begin{array}{c}f\\ m\end{array}\right\}` in the element. This constraint is calculated according to the equation [:ref:`éq 8.1-2 <éq 8.1-2>`]. But it should be noted that: * on the one hand [:ref:`éq 6.3-3 <éq 6.3-3>`]: :math:`{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 [:ref:`§10.8 <§10.8>`]. Stiffness matrices -------------------- The continuous expression for the material stiffness matrix [:ref:`éq 9.1-5 <éq 9.1-5>`] is: :math:`{B}^{T}{\Pi C\Pi }^{T}B`. For the finite element :math:`e`, we deduce the matrix relating the increase in displacement of node :math:`j` to the increase in internal force at node :math:`i`: :math:`{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 :math:`{\Pi C\Pi }^{T}` at the current Gauss point. The continuous expression for the geometric rigidity matrix [:ref:`éq 9.1-4 <éq 9.1-4>`] is: :math:`{Y}^{T}EY` where: :math:`Y=\left[\begin{array}{cc}\frac{d}{\text{ds}}1& 0\\ 0& \frac{d}{\text{ds}}1\\ 0& 1\end{array}\right],` and :math:`E` is given by the equation [:ref:`éq 9.1-2 <éq 9.1-2>`]. As for the material stiffness matrix, we deduce: :math:`{S}_{\text{géom}ij}^{e}={\int }_{e}{Y}_{i}^{\text{hT}}{EY}_{j}^{h}\text{ds},` **eq 10.2-2** where: :math:`{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 :math:`E` is calculated numerically at the current Gauss point. The continuous expression for the complementary stiffness matrix [:ref:`éq 9.1-5 <éq 9.1-5>`] is: :math:`{B}^{T}Z,` where :math:`Z` is given by the equation [:ref:`éq 9.1-3 <éq 9.1-3>`]. From this we deduce: :math:`{S}_{\text{compl}\mathrm{ij}}^{e}={\int }_{e}{B}_{i}^{\text{hT}}{N}_{j}Z\text{ds}`, where :math:`Z` is calculated numerically at the current Gauss point. Inertial forces ---------------- According to [:ref:`éq 7.3-1 <éq 7.3-1>`]: :math:`{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 :math:`e` at node :math:`i`. Inertia matrix ----------------- Let's say: :math:`{I}_{\mathrm{0,}\mathrm{iner}}=\left[\begin{array}{cc}\rho A1& 0\\ 0& {J}_{\rho }\end{array}\right]` and :math:`{M}_{0\mathrm{iner}\mathrm{ij}}^{e}={\int }_{e}{N}_{i}{N}_{j}{I}_{0\mathrm{iner}}\mathrm{ds}` We can see, from [:ref:`éq 7-4 <éq 7-4>`] and [:ref:`éq 9.2.1-1 <éq 9.2.1-1>`], that :math:`\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{iner}}^{e}` is indeed the inertia matrix of the element :math:`e` for the translational movement (diagonal sub-matrix :math:`\rho A1` of :math:`{I}_{0\mathrm{iner}}`). But according to [:ref:`éq9.2.2.1‑1 <éq9.2.2.1‑1>`] and [:ref:`éq 9.2.2.2-1 <éq 9.2.2.2-1>`], **it's not** the rotational inertia matrix. However, matrix :math:`{M}_{o\text{iner}}`, an assembly of :math:`{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** :math:`({\omega }_{o}=0)`. Indeed, we then have, according to [:ref:`éq 7.3-1 <éq 7.3-1>`]: :math:`{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: :math:`{I}_{\rho }={J}_{\rho }` Let's call :math:`J\text{'}` the sum of the two matrices [:ref:`éq 9.2.2.1-1 <éq 9.2.2.1-1>`] and [:ref:`éq 9.2.2.2-1 <éq 9.2.2.2-1>`] and set: :math:`{I}_{\text{rot}}=\left[\begin{array}{cc}0& 0\\ 0& J\text{'}\end{array}\right]` and :math:`{M}_{\text{rot}ij}^{e}={\int }_{e}{N}_{i}{N}_{j}{I}_{\text{rot}}\text{ds}` Let's also ask: :math:`{{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: :math:`{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}`. External forces given -------------------------- According to [:ref:`éq 8.1-3 <éq 8.1-3>`]: :math:`{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 :math:`i` of element :math:`e` which is equivalent to the distributed external forces. Linear iteration system ---------------------------- By discretization in finite elements, equation [:ref:`éq 9-1 <éq 9-1>`] gives: :math:`{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 [:ref:`§10.2 <§10.2>`] and [:ref:`§10.4 <§10.4>`]: :math:`{\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 [:ref:`éq 9-2 <éq 9-2>`], which must be verified by all :math:`\delta \phi`, therefore leads, at iteration :math:`n` from time :math:`i`, to the following linear system in :math:`\Delta {\phi }_{i}^{n+1}`: :math:`\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 :math:`{M}_{o\text{iner}}^{e},`, which are independent of the displacement, are assembled only once to constitute the global matrix :math:`{M}_{o\text{iner}},` which is used in particular to calculate the initial acceleration [:ref:`§10.4 <§10.4>`]. The three elementary stiffness matrices :math:`{S}_{\text{mat}}^{e},{S}_{\text{geom}}^{e}` and :math:`{S}_{\text{compl}}^{e}`, the rotation inertia matrix, :math:`{M}_{\text{rot}}^{e}`, which all four depend on the displacement, and the corrective rotational inertia matrix :math:`-\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 :math:`{\tilde{K}}_{i}^{n}\text{.}` So the linear system [:ref:`éq 10.6-1 <éq 10.6-1>`] becomes: :math:`\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: :math:`\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 :math:`{K}_{i}^{n}` is the assembly of the only stiffness matrices at iteration :math:`n` of "instant" :math:`i` [:ref:`§8.2 <§8.2>`]: :math:`{S}_{\text{mat,}i}^{n}+{S}_{\text{geom,}i}^{n}+{S}_{\text{compl,}i}^{n}\text{.}` Displacement, speed, and acceleration update -------------------------------------------------------------- The treatment of translational motion is classical; that of rotational motion is done using quaternions [:ref:`A8 `]. Translational movement ~~~~~~~~~~~~~~~~~~~~~~~~~~~ The formulas [:ref:`éq A3.1-3 <éq A3.1-3>`] and [:ref:`éq A3.1-4 <éq A3.1-4>`] are applied. 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 :math:`i-1` to the current iteration from time :math:`i` [:ref:`§A6 <§A6>`]. The update of the rotation vectors is based on the following property of quaternions [:ref:`§A8.4 <§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 [:ref:`§An8.6 <§An8.6>`]: :math:`\begin{array}{}(\Delta {x}_{o},\Delta x)=Q(\Delta {\Psi }_{i}^{n+1})\\ ({x}_{o},x)=Q({\Psi }_{i}^{n}),\end{array}` He comes: :math:`{\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 [:ref:`éq An6-2 <éq An6-2>`], if: :math:`({y}_{o},y)=Q({\Psi }_{i-\mathrm{1,}i}^{n}),` so: :math:`{\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 [:ref:`§4.2 <§4.2>`]: :math:`{R}_{i}^{n+1}=\text{exp}({\stackrel{ˆ}{\Psi }}_{i}^{n+1}),` which is calculated according to [:ref:`éq 4.2-1 <éq 4.2-1>`]. Finally, angular speed and acceleration are updated by the relationships [:ref:`éq A3.2-3 <éq A3.2-3>`] and [:ref:`éqA3.2‑4 <éqA3.2‑4>`]. Curvature variation vector update -------------------------------------------- The vector :math:`\chi`, which defines rotational deformation [:ref:`§6.2 <§6.2>`], should only be calculated at Gauss points. In the*Code_Aster*, it is treated computationally as an "internal variable". According to [:ref:`éq 6.2-2 <éq 6.2-2>`]: :math:`{\stackrel{ˆ}{\chi }}_{i}^{n}=\frac{d}{\text{ds}}({R}_{i}^{n})\text{.}{R}_{i}^{\mathrm{nT}}\text{.}` And, according to [:ref:`éq 4.2-6 <éq 4.2-6>`]: :math:`{\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: :math:`\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 [:ref:`bib2 `] Appendix B that the axial vector of this last one matrix is: :math:`\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: :math:`{\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{.}` Initialization before iterations ------------------------ In the dynamic case, if the load is constant over time, iterations can only start at time :math:`i` if we initialize some of the displacement fields, speed and acceleration to values different from those of time :math:`i-1`. These initializations are done as follows. Translational movement ~~~~~~~~~~~~~~~~~~~~~~~~~~~ :math:`{x}_{o,i}^{o}={x}_{o,i-1}\text{.}` Then, according to equation [:ref:`éq An3.1-1 <éq An3.1-1>`]: :math:`{\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 [:ref:`éq An3.1-2 <éq An3.1-2>`]: :math:`{\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{.}` Rotational movement ~~~~~~~~~~~~~~~~~~~~~~~ Expressions similar to the preceding ones are taken: :math:`{\Psi }_{i}^{o}={\Psi }_{i-1}` **eq 10.9.2-1** :math:`{\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** :math:`{\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 [:ref:`éq 10.9.2-2 <éq 10.9.2-2>`] and [:ref:`éq 10.9.2-3 <éq 10.9.2-3>`] make sense because all the vectors in them are in the vector space tangent to :math:`\text{SO}(3)` by :math:`{R}_{i-1}`. As consequences of equation [:ref:`éq 10.9.2-1 <éq 10.9.2-1>`]: :math:`{R}_{i}^{o}={R}_{i-1}` and: :math:`{\Psi }_{i-\mathrm{1,}i}^{o}=0\text{.}` We can see that at the first moment :math:`(i=1)`, the initializations require the knowledge of the initial acceleration :math:`{\left\{\begin{array}{c}{\ddot{x}}_{o}\\ \dot{\omega }\end{array}\right\}}_{o},` whose calculation is indicated in [:ref:`§10.4 <§10.4>`].