RK4 implementation

Wind & Water Structural Analysis
1

Dear all,

I have written MATLAB code to simulate wind turbine. I have considered only one mode (fore-aft) of tower and verified results with FAST simulation results. They are agreeing reasonably well. In this case I have used RK4 numerical integration in both MATLAB and FAST with a time step of 0.05 sec. To improve fidelity of MATLAB model, I have included second bending mode (fore-aft) of the tower also. Now, if i use the same integrator settings (RK4, stepsize 0.05) tower top displacement results are blowing up (to 10^200). In order to get meaningful results I had to reduce step size to 0.0001s but the results are meaningful but the computation time is prohibitively high (in days). But FAST is able to give meaningful results even with 0.01s timestep (RK4) even the time taken by FAST is much less (~ Hour). I am trying to understand the reason for this.
My question is, RK4 require to make 4 computations (K1, K2, K3 and K4) to move one step in time, right? Does FAST compute induction factors at all the 4 intermediate steps? Or does it do some sort of an interpolation? Because this is makes significant difference to the computational time (Approx ~0.08s for induction factor computation of at all airfoils at one time step). In my computations, I calculate induction factor at all the 4 intermediate steps. This makes my code significantly slower.

2

Dear Kumara,

While the ElastoDyn structural dynamics module of FAST can use an RK4 integrator, the AeroDyn aerodynamics module of FAST does not use an RK4 integrator. Instead, the AeroDyn aerodynamic calculations happen every DTAero seconds. In addition to RK4, ElastoDyn also has supports a 4th-order Adams-Bashforth (AB4) integrator and a 4th-order Adams-Bashforth-Moulton (ABM4) integrator. Regardless, the aerodynamic calculations are not typically the computational bottleneck of FAST.

RK4, AB4, and ABM4 are not stiff integrators, which means that the required time step (dt) is dictated by the highest natural frequency (f_max). For these integrators, my rule if thumb is dt <= 1/(10*f_max). Your time step sounds very small, which implies a very high f_max; otherwise, I would expect a bug in the implementation.

Best regards,

3

Dear Jason,

It was a bug in my code causing the trouble. Fixed it. Now I am able to get meaningful results with 0.02s (about 1/10 of the period corresponding to highest frequency) .
Sorry for the trouble.
Regards,
Kumara