Numerical solution to a complicated ODE

[Satchmo]

Background: For an oscillating wing that will be used in a tidal generator I am trying to model the position (as function of time) of an underwater airfoil whose movement is limited to vertical motion. The free body diagram is very simple. When the airfoil is positioned for upward travel, the only upward force is lift, and the downwards forces are drag and the force from an electrical generator. Assuming the airfoil to be neutrally buoyant, we have

F_L - F_D - F_T = (m+vm)*$$\ddot{y}$$

(vm is virtual mass or added mass, an effect of accelerating a body through fluid)
with y(0)=0 and $$\dot{y}$$(0) = 0
However, when the airfoil is in motion its vertical velocity changes the effective angle of attack and effective velocity seen by the airfoil, so that

F_L = A * sqrt(Vs² + ($$\dot{y}$$)²) * (C_L[ $$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)] * cos(atan($$\frac{\dot{y}}{Vs}$$))

F_D = A * sqrt(Vs² + ($$\dot{y}$$)²) * (C_D[ $$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)] * sin(atan($$\frac{\dot{y}}{Vs}$$))

where A is a constant that has to do with airfoil geometry, C_L[...] and C_D[...] are the coefficient of lift and drag for the instantaneous angle of attack (angle of attack given by [$$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)])
For C_L and C_D I have in MATLAB a smoothing spline functions so that I can input an angle and have the coefficient returned.

I would prefer to be able to use that MATLAB function to return the drag and lift coefficients, but a rough approximation can be given by

C_L = 2 * pi * ($$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)
and C_D = 0.013 * ($$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)


With these approximations, the lift and drag forces become

F_L = A * sqrt(Vs² + ($$\dot{y}$$)²) * (2*pi*[ $$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)] * cos(atan($$\frac{\dot{y}}{Vs}$$))


F_D = A * sqrt(Vs² + ($$\dot{y}$$)²) * (0.013 * [ $$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)] * sin(atan($$\frac{\dot{y}}{Vs}$$))

and the entire equation is:
A * sqrt(Vs² + ($$\dot{y}$$)²) * (2*pi*[ $$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)] * cos(atan($$\frac{\dot{y}}{Vs}$$)) - A * sqrt(Vs² + ($$\dot{y}$$)²) * (0.013 * [ $$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)] * sin(atan($$\frac{\dot{y}}{Vs}$$)) - F_T = (m+vm)*$$\ddot{y}$$

This is the first time I've tried to numerically solve a DE. I've looked into Matlab's ODE solvers (ode45 and ode113), but they require having the differential equations given as
y=..., $$\dot{y}$$= ..., $$\ddot{y}$$..., which I cannot do for this problem. I've never used Maple or Mathematica before, but web searches seem to show that they require the same format for solving ODEs. Therefore, I'm stuck.

Can anyone point me in the direction of a program that could possibly solve this ODE? That or ANY other advice would be greatly appreciated!

 

[Lord Crc]

Dividing the entire equation by (m+vm) gets you $$\ddot{y} = f(y, \dot{y})$$. You know that $$dy/dt} = \dot{y}$$ and $$d\dot{y}/dt = \ddot{y} = f(y, \dot{y})$$ So now you got a system of two (coupled) ODEs that you can try to solve, one for $$y$$ and the other for $$\dot{y}$$. I've never used Matlab, but from the reference it seems you can do something like

function dy = airfoil(t,y)
  dy = zeros(2,1);    % a column vector
  dy(1) = y(2); % dy/dt = \dot{y}
  dy(2) = f(y(1), y(2)); % d\dot{y}/dt = \ddot{y} = f(y, \dot{y})

It's a bit late here so I might have missed something or messed it up :)