Solve wind turbine with ode solver

Hey, I'm trying to solve the following differential equation with an ode solver:

   tspan = [ta te];
      ic = 1;
      opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
      [t,u]=ode113(@(t,u) fu1_SS(t,u,u0,Pt,A,zusFG),[ta te],u0,opts);

function u=fu1_SS(t,u,u0,Pt,A,zusFG) 
A=A.';
A = interp1(Pt,A,t); % Interpolate the data set (P,A) at time t
A=A.';

u=zeros((nodes*ld+3*zusFG)*2,1);  % 1:nodes*ld Displacement, nodes*ld+1:nodes*ld*2 Velocity
du=zeros((nodes*ld+3*zusFG)*2,1); % 1:nodes*ld Velocity, nodes*ld+1:nodes*ld*2 Acceleration  
du(1:nodes*ld+3*zusFG,1)=u(nodes*ld+3*zusFG+1:(nodes*ld+3*zusFG)*2,1);   % du(1)=u(2)
u(1:nodes*ld+3*zusFG,1)=K\(A-M*du(nodes*ld+3*zusFG+1:(nodes*ld+3*zusFG)*2,1)-C*u(nodes*ld+3*zusFG+1:(nodes*ld+3*zusFG)*2,1));   %Displacement u(1)=K\(P-M*du(2)-C*u(2))

end

nodes are the Nodes of the wind turbine, ld are the degrees of freedom, zusFG are additional degrees of freedom of the soil. M, C and K are matrices which declare the system (mass, damping and stiffness). Pt is the time vector to interpolate A. A is the wind load, which consists out of a load vector for every node. The Problem is, that the displacement is becoming larger and larger instead of oscillating about one constant value...

Does anybody have an idea where the problem is?

Thanks for you help in advance!

Answer 1

[Torsten]

I don't understand this line in your code.

u(1:nodes*ld+3*zusFG,1)=K\(A-M*du(nodes*ld+3*zusFG+1:(nodes*ld+3*zusFG)*2,1)-C*u(nodes*ld+3*zusFG+1:(nodes*ld+3*zusFG)*2,1));   %Displacement u(1)=K\(P-M*du(2)-C*u(2))

u is the solution variable - you can't determine or prescribe it here. You will have to prescribe velocity in

du(1:nodes*ld+3*zusFG)

and acceleration in

du( nodes*ld+3*zusFG+1:(nodes*ld+3*zusFG)*2)

if your function has the form

function du=fu1_SS(t,u,u0,Pt,A,zusFG)

Best wishes

Torsten.

[Lisanne Meinerzhagen]

Hey Torsten, thanks again, that did already help!! Only problem is now, that even if I'm using an stiff ode solver (ode15s), matlab interrupts the calculation at t approximately equal to 0.27s, although my tspan is 45s... and now my deflections become so small, the output vector is only giving out zeros... do you have a suggestion which ode solver I could use? I'll send you the code again:

   tspan = [ta te];
      opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
      [t,u]=ode15s(@(t,u) fu1_SS(t,u,Pt,A,zusFG),[ta te],u0,opts);

function du=fu1_SS(t,u,Pt,A,zusFG) 
A=A.';
A = interp1(Pt,A,t); % Interpolate the data set (P,A) at time t
A=A.';
du(1:nodes*ld+3*zusFG,1)=u(nodes*ld+3*zusFG+1:(nodes*ld+3*zusFG)*2,1);   % du(1)=u(2)
du(nodes*ld+1+3*zusFG:(nodes*ld+3*zusFG)*2,1)=M\(A-K*u(1:nodes*ld+3*zusFG,1)-C*u(nodes*ld+3*zusFG+1:(nodes*ld+3*zusFG)*2,1));   %Acceleration %du(2)=M\(P-K*u(1)-C*u(2)) 
end

Thank you again!