Home > Theory > Numerical Issues > Algorithms and Convergence Criteria

Algorithms and convergence criteria

The slender structure is modelled by a series of nodes connected altogether by elements of different kinds. Each nodes is given a position vector \(\vec{X_0}\) and a rotation vector \(\vec{\Theta_x}\) . Hereafter the vector \(\vec{X}\) stands for the positions as well as for the rotations.

In static, the total virtual work \(G\) consists in external and internal efforts contributions : $$ G(\vec{X})=\vec{Q_{int}}.\delta\vec{X}+\vec{Q_{ext}}.\delta\vec{X} $$

From initial values the problem is to find the positions for which \(G(\vec{X})=0\).

Newton-Raphson algorithm

The iterative Newton-Raphson scheme is used to deals with differential non- linearities (large displacements and rotations in 3D, material laws, contact forces quantification...).

For a given set of positions, \(G(\vec{X})\) is evaluated. If \(G(\vec{X})=0\) is not reached, a new set is calculated from the preceeding using the stiffness matrix which gives the variation \(\Delta G(\vec{X})\) of \(G(\vec{X})\) with a displacement \(\Delta\vec{X}\).

Different options are available and may be chosen with the keyword *NUMINT.

Matrix storage :

First, the way the matrix \(\Delta G(\vec{X})\) is stored may influence the calculations. Indeed two storage possibilities are offered : a band type and a profile type. The first one is the default . For slender structures, a node is usually in relation with its nearest neighbours. The numerical translation of this property is a sparsed matrix. The effective terms of the matrix are not far from one another and it is possible to define a band length.

• Lband = 12 for binodal beams,

• Lband = 18 for trinodal beams,

• Lband = 6 for bars and cables.

The goal is to store the effective terms only so as to prevent useless operations. In a band storage mode, the highest band length is evaluated and Lband colons are stored from the first diagonal. In a profile storage mode, only the effective terms are stored which requires a working vector linking the real matrix and the stored one.

The second mode is very usefull everytime special connections are created between different seabed connections. In that case, the band length is greatly increased for only a few lines of the matrix. The band storage leads to a much higher cputime consumption.

Additionnal matrix :

In the default option, an additionnal matrix is added to the calculated one. This prevents nul terms when the inversion is done. This is especially the case with loose cables on the seabed. However with beam elements, this option is not recommended.

Fixed point method

The Newton-Raphson algorithm can not deal with some non-linearities which may not be derived. This actually happens with riser clashes and the wake effect... Indeed the detection of contacts as well as the wake flow evaluation are treated with a fixed point method. This means that some iterative procedures have been implemented above the Newton-Raphson algorithm.

For instance, in case of risers clashes, the aera of neighbourhood are first determined. According to them, contact elements are created with their own normal vector. At that stage, the virtual work principle is solved taking into account contact reactions. Their modulus are adjusted but the contact elements are fixed as well as the contact directions even if the structures move. Thats why getting out of the Newton-Raphson algorithm, the location of contact elements is up-dated and the system of equations is solved again.

Finally, convergence is reached when only one iteration is enough to get the equilibrium. For practical reasons, the number of loops for the fixed point method is limited with the keyword *CONTACT. The final convergence state is supposed to be acceptable.

Newmark's algorithm :

In dynamic, The behaviour of the structure is a series of equilibrium states between the external, internal and inertia efforts contributions : $$ G(\vec{X},\dot{\vec{X}},\ddot{\vec{X}})=\vec{Q_{int}}(\vec{X}).\delta\vec{X}+\vec{Q_{ext}}(\vec{X},\dot{\vec{X}},\ddot{\vec{X}}).\delta\vec{X}+\vec{Q_{inner}}(\vec{X},\dot{\vec{X}},\ddot{\vec{X}}).\delta\vec{X} $$

Newmark's algorithm allows to express the unknowns \(\dot{\vec{X}}\) and \(\ddot{\vec{X}}\) as functions of \(\vec{X}\) and by so doing to come back to a single unkown problem.

Let \(\vec{X}^*_{t+dt}\) and \(\dot{\vec{X}}^*_{t+dt}\), be the predictors at time t+dt calculated with the following formula :

\[ \vec{X}^*_{t+dt}=\vec{X}_{t}+\Delta t\dot{\vec{X}}_{t}+(1-2\beta)\frac{\Delta t^2}{2}\ddot{\vec{X}}_{t} \]

\[ \dot{\vec{X}}^*_{t+dt}=\dot{\vec{X}}_{t}+(1-\gamma)\Delta t\ddot{\vec{X}}_{t} \]

For a given position vector \(\vec{X}_{t+dt}\) at t+dt, speeds and accelerations come :

\[ \dot{\vec{X}}_{t+dt}=\dot{\vec{X}}^*_{t+dt}+\frac{\gamma}{\beta\Delta t}(\vec{X}_{t+dt}-\vec{X}_{t+dt}^*)=fv(\vec{X}_{t+dt}) \]

\[ \ddot{\vec{X}}_{t+dt}=\frac{1}{\beta\Delta t^2}(\vec{X}_{t+dt}-\vec{X}_{t+dt}^*)=fa(\vec{X}_{t+dt}) \]

With \(\beta\) and \(\gamma\) parameters of the numerical scheme.

The virtual work principle then shrinks to a single unknown function like in static :

\[ G(\vec{X}_{t+dt},\dot{\vec{X}}_{t+dt},\ddot{\vec{X}}_{t+dt})=G(\vec{X}_{t+dt}) \]

A new stiffness may be calculated \(\tilde {K}\):

\[ \Delta G=K\Delta\vec{X}_{t+dt}+M\Delta\ddot{\vec{X}}_{t+dt}+C\Delta\dot{\vec{X}}_{t+dt}=\tilde{K}\Delta\vec{X}_{t+dt} \]

with the Newmark's formula :

\[ \Delta\dot{\vec{X}}_{t+dt}=\frac{\gamma}{\beta\Delta t}\Delta \vec{X}_{t+dt} \quad;\quad \Delta\ddot{\vec{X}}_{t+dt}=\frac{\gamma}{\beta\Delta t^2}\Delta \vec{X}_{t+dt} \]

and \(\Delta\ddot{\vec{X}}_{t+dt}=\frac{1}{\beta\Delta t^2}\Delta \vec{X}_{t+dt}\)

\[ \tilde{K}=K+\frac{1}{\beta\Delta t^2}M+\frac{\gamma}{\beta\Delta t}C \]

At each time step, the Newton-Raphson scheme is then used so as to find :

\[ G(\vec{X}_{t+dt})=0 \]

At the beginning of each time step, the choice of initial values may help the convergence of the Newton-Raphson algorithm. Four possibilities are offered :

1) Select the predicted values as initial ones :

\[ \begin{cases} \vec{X}_{t+dt}^0 &=\vec{X}_{t+dt}^* \\ \dot{\vec{X}}^0_{t+dt} &=\dot{\vec{X}}_{t+dt}^* \\ \ddot{\vec{X}}^0_{t+dt} &= \vec0 \end{cases} \]

2) Keep the positions of the previous time step :

\[ \begin{cases} \vec{X}_{t+dt}^0 &=\vec{X}_{t} \\ \dot{\vec{X}^0}_{t+dt} &=\dot{\vec{X}}_{t+dt}^*+\frac{\gamma}{\beta\Delta t}(\vec{X}_{t}-\vec{X}_{t+dt}^*) \\ \ddot{\vec{X}^0}_{t+dt} &=\frac{1}{\beta\Delta t^2}(\vec{X}_{t}-\vec{X}_{t+dt}^*) \\ \end{cases} \]

3) Keep the accelerations of the previous time step :

\[ \begin{cases} \vec{X}_{t+dt}^0 &=\vec{X}_{t+dt}^*+\dot{\vec{X}}_t\Delta t+\frac{\Delta t^2}{2}\ddot{\vec{X}}_{t} \\ \dot{\vec{X}^0}_{t+dt} &=\dot{\vec{X}}_{t+dt}^*+\ddot{\vec{X}}_{t}\gamma\Delta t \\ \ddot{\vec{X}^0}_{t+dt} &=\ddot{\vec{X}}_{t} \end{cases} \]

4) Last option, equivalent to 1) if \(\beta=\frac{1}{2}\):

\[ \begin{cases} \vec{X}_{t+dt}^0 &=\vec{X}_{t}+\dot{\vec{X}}_{t}\Delta t \\ \dot{\vec{X}}^0_{t+dt} &=\dot{\vec{X}}_{t+dt}^*+\ddot{\vec{X}}_{t}(1-\frac{1}{2\beta})\gamma \\ \ddot{\vec{X}}^0_{t+dt} &= \ddot{\vec{X}}_{t}(1-\frac{1}{2\beta}) \end{cases} \]

These options are available with the keyword *DYNAMIC. The option by default is the first one but the acceleration is the mean value between the actual acceleration and the previous one. This introduces a numerical damping.