Time domain simulation

A simulation in the time domain is performed by solving for the dynamic equation of motion for each time step. At a time step
$$n$$
, the dynamic equation of motion is given by
$$M\ddot{d_n} + C\dot{d_n}+f^{int}_n=f^{ext}_n$$


where 
 
$$M$$
 is the mass matrix
  
$$C$$
 is the damping matrix
  
$$d$$
 is the vector containing the displacement at each node of the system
  
$$f^{int}$$
 is the vector containing the internal forces (i.e. response forces) at each node of the system
 and   
$$f^{ext}$$
 is the vector containing the external forces (i.e. excitation forces) of the system

It is common to talk about the left hand side of the equation (referring to the terms
$$M\ddot{d_n} + C\dot{d_n}+f^{int}_n$$
) and the right-hand side of the equation (referring to the term
$$f^{ext}_n$$
). You might see this wording in the litterature or elsewhere in this User Manual. For practical purposes however, here we will rearrange the equation such that it becomes

$$M\ddot{d_n} + C\dot{d_n}+f^{int}_n-f^{ext}_n=0$$


If the system is assumed to be known at time step 
$$n-1$$
, it is possible to substract the equation at 
$$t=t_{n-1}$$
 to that at 
$$t = t_n$$
 and obtain the incremental form of the equation of motion:

$$M\Delta\ddot{d}_n+C\Delta\dot{d}_n+\Delta f^{int}_n-\Delta f^{ext}_n=0$$

 where 
$$\Delta X_n = X_{n}-X_{n-1}$$
. The internal and external forces can be linearized such that
$$\Delta f^{int}_n=\frac{\partial f^{int}}{\partial d}\Delta d_n = K^{int}_n\Delta d_n$$
 and
$$\Delta f^{ext}_n =\frac{\partial f^{ext}}{\partial d}\Delta d_n = K^{ext}_n\Delta d_n$$

where 
$$K^{int}$$
 is called the consistent tangent matrix and 
$$K^{ext}$$
 is the load stiffness matrix. Both these matrices are in general dependent on the displacement
$$d$$
Note that in a non-linear analysis, the matrices 
$$M$$
,
$$C$$
$$K^{int}$$
and 
$$K^{ext}$$
are dependent on the displacement
$$d$$
. It is possible to neglect this dependency by running a Linear analysis.

The equation is then solved using the Newmark-Beta method (see https://en.wikipedia.org/wiki/Newmark-beta_method). It is possible to include numerical damping to the algorithm using the HHT method (see the Damping section). Note that when deriving the incremental form of the equation of motion, it is assumed that the mass and damping properties are constant over time, which is not correct due to their dependency on the displacement 
$$d$$
. In addition, the linearization of the internal and external forces introduces an error on the matrices 
$$K^{int}$$
and
$$K^{ext}$$
. This means that the equation of motion solved at 
$$t=t_n$$
 uses matrices of 
$$t_{n-1}$$
, which in a non-linear analysis will introduce an error 
$$r_n$$
 such as

$$M\Delta\ddot{d}_n+C\Delta\dot{d}_n+(K^{int}_n-K^{ext}_n)\Delta d_n = r_n $$


This error can be reduced by using the Newton-Raphson method (see https://en.wikipedia.org/wiki/Newton%27s_method). The Newton-Raphson iteration scheme reduces the residual error by repeatedly computing a new incremental displacement and the corresponding mass, damping and stiffness matrices. 

This iteration process is what makes the non-linear simulation be significantly more computationally intensive than the linear analysis: for each time step, several iterations of the Newton-Raphson scheme have to be performed before reaching a small enough residual error.

The Newton-Raphson iterations are stopped when the residual energy of the system is smaller than a certain criteria. From the residual error, for a given time step, the residual energy is defined as

$$R^X_E=\sqrt{\sum_{i=1}^{N_X}\left(r_n^i\mathbb{X}^i\right)}$$
 and 
$$R^\theta_E=\sqrt{\sum_{i=1}^{N_\theta}\left(r_n^i\mathbb{\Theta^i}\right)}$$
where 
$$\mathbb{X}$$
 and 
$$\mathbb{\Theta}$$
 are the translational and rotational degrees of freedom, respectively, 
$$R^X_E$$
 and 
$$R^\theta_E$$
 are the residual translational and rotational energies, respectively and 
$$N_X$$
 and 
$$N_\theta$$
 are the number of translational and rotational degrees of freedom, respectively. The superscript 
$$i$$
 indicates that the parameter applies to the 
$$i^{th}$$
degree of freedom. 
As can be seen from the previous equations, in order to consider that the equation of motion is in equilibrium, the residual error (and therefore the residual energies) must be small. If the Newton-Raphson method converges, the residual energy will decrease at each iteration. Several iterations are performed until the desired accuracy is reached. In Ashes, a treshold to consider that the equation has reached equilibrium is defined by the following criteria:

$$R^X_E+R^\theta_E<10^{\epsilon_a}$$
 

where 
$$\epsilon_a$$
 is an exponent defined as Energy tolerance (absolute) that can be modified by the user in the Convergence section of the Analysis tab in the Analysis parameters dialog. This criteria is called Absolute criteria.  

For the default templates under the default conditions, it typically takes 7-8 Newton-Raphson iterations to fulfill the default criteria. This can be checked in the Solver sensor 

If the absolute criteria is not fulfilled within a given number of iterations (defined in Maximum iterations, see Analysis), the time step is considered not converged. By default, when the time step is not converged, the simulation will stop, but it is possible to force Ashes to keep simulating even if convergence is not reached by modifying the On non-convergence parameter (see Analysis)

Note: results obtained from a simulation that does not converge can be inaccurate. Be sure to double-check the accuracy of such results

1 Time step

The time step for the FEM solver can be manually set in the Analysis tab of the Analysis parameter window, by setting the parameter Timestep scheme to User defined.
If that parameter is set to Calculated, Ashes computes a time step based on the radius 
$$r$$
of the current model. This time step is calculated so that the rotation of the rotor during one time step is not more than 2 degrees. The rotational speed of the rotor is calculated assuming a wind speed 
$$V=11.4\text{ m}\cdot{s}^{-1}$$
and a tip speed ratio 
$$TSR = 6$$

For the NREL 5-MW wind turbine for example, with an (initial, projected) rotor radius 
$$r=62.941\text{ m}$$
, the time step would be
$$dt_{unrounded} = \frac{2\pi r}{TSR\cdot V}\frac{2}{360}=0.0321\text{ s}$$

Ashes then truncates the time step to keep only one significant digit, giving 
$$dt = 0.03\text{ s}$$


Note: If the model has no rotor, the calculated time step is set to 
$$0.03\text{ s}$$