5. Geometric stiffness - Prestressed structure

This matrix is calculated by the rigi_geom option. It is used to treat buckling problems or vibrations of prestressed structures. In the case of a prestressed structure, which is therefore subject to initial forces (known and independent of time), it is not possible to neglect in the equilibrium equation the terms introduced by the change in geometry from the unstressed state to the prestressed state. This change in geometry only modifies the equilibrium equation by adding a term that is a function of displacements and prestress whose associated matrix is called the geometric rigidity matrix and which is expressed by:

\({W}_{G}=\underset{{V}_{o}}{\int }\frac{\partial {u}_{k}^{\mathrm{3D}}}{\partial {x}_{i}}{\sigma }_{\text{ij}}^{o}\frac{\partial {v}_{k}^{\mathrm{3D}}}{\partial {x}_{j}}\text{dV}\)

where \({\sigma }_{\text{ij}}^{o}\) refers to the prestress tensor. This term appears naturally if we introduce the deformation tensor of GREEN - LAGRANGE into the virtual deformation work:

\(\begin{array}{}{E}_{\text{xx}}={\varepsilon }_{\text{xx}}+{\eta }_{\text{xx}}=\frac{\partial {u}_{x}^{\mathrm{3D}}}{\partial x}+\frac{1}{2}\left[{(\frac{\partial {u}_{x}^{\mathrm{3D}}}{\partial x})}^{2}+{(\frac{\partial {u}_{y}^{\mathrm{3D}}}{\partial x})}^{2}+{(\frac{\partial {u}_{z}^{\mathrm{3D}}}{\partial x})}^{2}\right]\\ 2{E}_{\text{xy}}=2{\varepsilon }_{\text{xy}}+2{\eta }_{\text{xy}}=\frac{\partial {u}_{x}^{\mathrm{3D}}}{\partial y}+\frac{\partial {u}_{y}^{\mathrm{3D}}}{\partial x}+\left[\frac{\partial {u}_{x}^{\mathrm{3D}}}{\partial x}\frac{\partial {u}_{x}^{\mathrm{3D}}}{\partial y}+\frac{\partial {u}_{y}^{\mathrm{3D}}}{\partial x}\frac{\partial {u}_{y}^{\mathrm{3D}}}{\partial y}+\frac{\partial {u}_{z}^{\mathrm{3D}}}{\partial x}\frac{\partial {u}_{z}^{\mathrm{3D}}}{\partial y}\right]\\ 2{E}_{\text{xz}}=2{\varepsilon }_{\text{xz}}+2{\eta }_{\text{xz}}=\frac{\partial {u}_{x}^{\mathrm{3D}}}{\partial z}+\frac{\partial {u}_{z}^{\mathrm{3D}}}{\partial x}+\left[\frac{\partial {u}_{x}^{\mathrm{3D}}}{\partial x}\frac{\partial {u}_{x}^{\mathrm{3D}}}{\partial z}+\frac{\partial {u}_{y}^{\mathrm{3D}}}{\partial x}\frac{\partial {u}_{y}^{\mathrm{3D}}}{\partial z}+\frac{\partial {u}_{z}^{\mathrm{3D}}}{\partial x}\frac{\partial {u}_{z}^{\mathrm{3D}}}{\partial z}\right]\end{array}\)

In the expression of these deformations, quadratic terms \({(\frac{\mathrm{\partial }{u}_{x}^{\mathrm{3D}}}{\mathrm{\partial }x})}^{2}\text{,}\frac{\mathrm{\partial }{u}_{x}^{\mathrm{3D}}}{\mathrm{\partial }x}\frac{\mathrm{\partial }{u}_{x}^{\mathrm{3D}}}{\mathrm{\partial }y}\text{et}\frac{\mathrm{\partial }{u}_{x}^{\mathrm{3D}}}{\mathrm{\partial }x}\frac{\mathrm{\partial }{u}_{x}^{\mathrm{3D}}}{\mathrm{\partial }z}\) are neglected here, following the assumption commonly made by most authors [bib]. For a beam model, the initial stress tensor is reduced in the local axes of the beam to components \({\sigma }_{\text{xx}}\text{,}{\sigma }_{\text{xy}}\) and \({\sigma }_{\text{xz}}\). We use the kinematics introduced in [§2]:

\(\{\begin{array}{}{u}_{x}^{\mathrm{3D}}(x,y,z)={u}_{G}(x)+z{\theta }_{y}(x)-y{\theta }_{z}(x)+\omega (y,z){\theta }_{x,x}(x)\\ {u}_{y}^{\mathrm{3D}}(x,y,z)={v}_{C}(x)-(z-{z}_{c}){\theta }_{x}(x)\\ {u}_{z}^{\mathrm{3D}}(x,y,z)={w}_{C}(x)+(y-{y}_{c}){\theta }_{x}(x)\end{array}\)

and the expression of generalized efforts according to constraints:

\({N}^{0}=\underset{S}{\int }{\sigma }_{\text{xx}}^{o}\text{ds}{V}_{{}_{y}}^{0}=\underset{S}{\int }{\sigma }_{\text{xy}}^{o}\text{ds}{V}_{{}_{z}}^{0}=\underset{S}{\int }{\sigma }_{\text{xz}}^{o}\text{ds}{M}_{{}_{y}}^{0}=\underset{S}{\int }z{\sigma }_{{}_{\text{xx}}}^{0}\text{ds}{M}_{{}_{z}}^{0}=\underset{S}{\int }-y{\sigma }_{{}_{\text{xx}}}^{0}\text{ds}\)

It is also assumed that \({N}^{0}\), \({V}_{{}_{y}}^{0}\), \({V}_{{}_{z}}^{0}\) are constant in the discretized element (which is inaccurate for example for a vertical beam subject to its own weight). The moments are supposed to vary linearly:

\({M}_{y}^{0}=\left({M}_{\mathit{y2}}^{0}-{M}_{\mathit{y1}}^{0}\right)\frac{x}{L}+{M}_{\mathit{y1}}^{0}\) \(\frac{\partial {M}_{y}^{0}}{\partial x}-{V}_{z}^{0}=0\)

\({M}_{z}^{0}=\left({M}_{\mathit{z2}}^{0}-{M}_{\mathit{z1}}^{0}\right)\frac{x}{L}+{M}_{\mathit{z1}}^{0}\) \(\frac{\partial {M}_{z}^{0}}{\partial x}+{V}_{y}^{0}=0\)

These hypotheses make it possible to express \({W}_{G}\) for a right beam with warping in the following way:

\(\begin{array}{c}{W}_{G}={\int }_{0}^{L}{N}^{0}\left({v}_{,x}\delta {v}_{,x}+{w}_{,x}\delta {w}_{,x}\right)+{N}^{0}\left(\frac{{I}_{y}+{I}_{z}}{A}+{y}_{c}^{2}+{z}_{c}^{2}\right){\theta }_{x,x}\delta {\theta }_{x,x}\\ +{z}_{c}{N}^{0}\left({v}_{,x}\delta {\theta }_{x,x}+{\theta }_{x,x}\delta {v}_{,x}\right)-{y}_{c}{N}^{0}\left({w}_{,x}\delta {\theta }_{x,x}+{\theta }_{x,x}\delta {w}_{,x}\right)\\ -{M}_{y}^{0}\left({w}_{,x}\delta {\theta }_{,x}+{\theta }_{,x}\delta {w}_{,x}\right)-{M}_{z}^{0}\left({v}_{,x}\delta {\theta }_{,x}+{\theta }_{,x}\delta {v}_{,x}\right)\\ -\left({y}_{c}{V}_{y}^{0}-{z}_{c}{V}_{z}^{0}\right)\left({\theta }_{x,x}\delta {\theta }_{x}+{\theta }_{x}\delta {\theta }_{x,x}\right)+{V}_{y}^{0}\left({w}_{,x}\delta {\theta }_{x}+{\theta }_{x}\delta {w}_{,x}\right)\\ -{V}_{z}^{0}\left({v}_{,x}\delta {\theta }_{x}+{\theta }_{x}\delta {v}_{,x}\right)+\left(\left({\mathrm{2y}}_{c}-\frac{{I}_{{\mathit{yr}}^{2}}}{{I}_{z}}\right){M}_{z}^{0}+\left(-{\mathrm{2z}}_{c}+\frac{{I}_{{\mathit{zr}}^{2}}}{{I}_{y}}\right){M}_{y}^{0}\right){\theta }_{x,x}\delta {\theta }_{x,x}\\ +\left(-\frac{{I}_{\mathit{yr2}}}{{I}_{z}}\frac{{\text{dM}}_{z}^{0}}{\text{dx}}+\frac{{I}_{\mathit{zr2}}}{{I}_{y}}\frac{{\text{dM}}_{y}^{0}}{\text{dx}}\right)\left({\theta }_{x,x}\delta {\theta }_{x}+{\theta }_{x}\delta {\theta }_{x,x}\right)\end{array}\)

with the terms:

\(\begin{array}{c}{I}_{{\mathit{yr}}^{2}}=\underset{S}{\int }y\left({y}^{2}+{z}^{2}\right)\text{ds}\\ {I}_{{\mathit{zr}}^{2}}=\underset{S}{\int }z\left({y}^{2}+{z}^{2}\right)\text{ds}\end{array}\)

which represent the non-symmetry of the section. In the case where the section has two axes of symmetry (so \(C\) is confused with \(G\)), these terms are zero.

Attention, these terms (which are called IYR2 and IZR2 in the AFFE_CARA_ELEM command) are not currently calculated by MACR_CARA_POUTRE. The user must therefore enter them using tubular values for each type of section (angle, rectangle, etc.).

In addition, to be able to deal with the problems of overturning thin beams, stressed essentially by bending moments and normal forces, it is necessary to add the hypothesis of moderate torsional rotations [bib], [bib].

This results in the following modification of the displacement field (only for calculating the geometric rigidity):

\({u}_{x}^{\mathrm{3D}}(x,y,z)={u}_{G}(x)+z({\theta }_{y}(x)+{\theta }_{x}(x){\theta }_{z}(x))-y({\theta }_{z}(x)-{\theta }_{x}(x){\theta }_{y}(x))+\omega (y,z){\theta }_{x,x}(x)\)

The origin of this expression cannot be detailed here. It is the subject of the thesis of DE VILLE DE GOYET [bib] on the buckling of beams with thin open sections. The hypothesis of moderate torsional rotations (and not infinitesimal) makes it possible to correctly model the torsional movement of a thin section beam (torsional - flexure coupling).

The hypothesis of moderate rotations leads to the addition of the term \({W}_{G}^{1}\) to \({W}_{{}_{G}}^{0}\):

\({W}_{G}^{1}=\frac{1}{2}{\int }_{o}^{L}-{M}_{z}^{o}\delta ({\theta }_{x}{\theta }_{y,x}+{\theta }_{y}{\theta }_{x,x})+{M}_{y}^{o}\delta ({\theta }_{x}{\theta }_{z,x}+{\theta }_{z}{\theta }_{x,x})+{V}_{y}^{o}\delta ({\theta }_{x}{\theta }_{y})+{V}_{z}^{o}\delta ({\theta }_{x}{\theta }_{z})\)

Finally, the geometric stiffness matrix is obtained by discretizing \({W}_{{}_{G}}={W}_{{}_{G}}^{0}+{W}_{{}_{G}}^{1}\) using the same interpolation functions as the stiffness matrix of [§4.4]. After calculating these matrices, you have to change the coordinate system as in [§4.5]. A matrix of geometric rigidity of the form is then obtained:

\({K}_{G}=(\begin{array}{cc}{A}_{1}& {A}_{2}\\ {A}_{2}& {A}_{3}\end{array})\)

The blocks of the matrix are explained below. To simplify expressions, we use:

\(\begin{array}{}{N}_{\text{ey}}^{0}=1\text{.}2\frac{{N}^{o}{e}_{y}}{L}{N}_{\text{ez}}^{0}=1\text{.}2\frac{{N}^{o}{e}_{z}}{L}\\ {\stackrel{ˉ}{M}}_{y}^{0}=\frac{{M}_{\mathrm{y1}}^{0}+{M}_{\mathrm{y2}}^{0}}{L}{\stackrel{ˉ}{M}}_{z}^{0}=\frac{{M}_{\mathrm{z1}}^{0}+{M}_{\mathrm{z2}}^{0}}{L}\\ \Delta {M}_{y}^{0}=\frac{{M}_{\mathrm{y1}}^{0}-{M}_{\mathrm{y2}}^{0}}{2}\Delta {M}_{z}^{0}=\frac{{M}_{\mathrm{z1}}^{0}-{M}_{\mathrm{z2}}^{0}}{2}\\ \tilde{k}={N}^{o}(\frac{{I}_{y}+{I}_{z}}{S}+{e}_{y}^{2}+{e}_{z}^{2})\\ {\tilde{I}}_{y}=-\frac{{I}_{\text{yr}2}}{{I}_{z}}+2{e}_{y}\\ {\tilde{I}}_{z}=\frac{{I}_{\text{zr}2}}{{I}_{y}}-2{e}_{z}\end{array}\)