# Numerical solution to a complicated ODE

1. Dec 6, 2009

### 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

FL - FD - FT = (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

FL = A * sqrt(Vs2 + ($$\dot{y}$$)2) * (CL[ $$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)] * cos(atan($$\frac{\dot{y}}{Vs}$$))

FD = A * sqrt(Vs2 + ($$\dot{y}$$)2) * (CD[ $$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)] * sin(atan($$\frac{\dot{y}}{Vs}$$))

where A is a constant that has to do with airfoil geometry, CL[...] and CD[...] 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 CL and CD 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

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

With these approximations, the lift and drag forces become

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

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

and the entire equation is:
A * sqrt(Vs2 + ($$\dot{y}$$)2) * (2*pi*[ $$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)] * cos(atan($$\frac{\dot{y}}{Vs}$$)) - A * sqrt(Vs2 + ($$\dot{y}$$)2) * (0.013 * [ $$\alpha$$ - atan($$\frac{\dot{y}}{Vs}$$)] * sin(atan($$\frac{\dot{y}}{Vs}$$)) - FT = (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!

File size:
3.8 KB
Views:
73
2. Dec 6, 2009

### 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
Code (Text):
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 :)

Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook