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(adsbygoogle = window.adsbygoogle || []).push({});

F_{L}- F_{D}- F_{T}= (m+vm)*[tex]\ddot{y}[/tex]

(vm is virtual mass or added mass, an effect of accelerating a body through fluid)

with y(0)=0 and [tex]\dot{y}[/tex](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^{2}+ ([tex]\dot{y}[/tex])^{2}) * (C_{L}[ [tex]\alpha[/tex] - atan([tex]\frac{\dot{y}}{Vs}[/tex])] * cos(atan([tex]\frac{\dot{y}}{Vs}[/tex]))

F_{D}= A * sqrt(Vs^{2}+ ([tex]\dot{y}[/tex])^{2}) * (C_{D}[ [tex]\alpha[/tex] - atan([tex]\frac{\dot{y}}{Vs}[/tex])] * sin(atan([tex]\frac{\dot{y}}{Vs}[/tex]))

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 [[tex]\alpha[/tex] - atan([tex]\frac{\dot{y}}{Vs}[/tex])])

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 * ([tex]\alpha[/tex] - atan([tex]\frac{\dot{y}}{Vs}[/tex])

and C_{D}= 0.013 * ([tex]\alpha[/tex] - atan([tex]\frac{\dot{y}}{Vs}[/tex])

With these approximations, the lift and drag forces become

F_{L}= A * sqrt(Vs^{2}+ ([tex]\dot{y}[/tex])^{2}) * (2*pi*[ [tex]\alpha[/tex] - atan([tex]\frac{\dot{y}}{Vs}[/tex])] * cos(atan([tex]\frac{\dot{y}}{Vs}[/tex]))

F_{D}= A * sqrt(Vs^{2}+ ([tex]\dot{y}[/tex])^{2}) * (0.013 * [ [tex]\alpha[/tex] - atan([tex]\frac{\dot{y}}{Vs}[/tex])] * sin(atan([tex]\frac{\dot{y}}{Vs}[/tex]))

and the entire equation is:

A * sqrt(Vs^{2}+ ([tex]\dot{y}[/tex])^{2}) * (2*pi*[ [tex]\alpha[/tex] - atan([tex]\frac{\dot{y}}{Vs}[/tex])] * cos(atan([tex]\frac{\dot{y}}{Vs}[/tex])) - A * sqrt(Vs^{2}+ ([tex]\dot{y}[/tex])^{2}) * (0.013 * [ [tex]\alpha[/tex] - atan([tex]\frac{\dot{y}}{Vs}[/tex])] * sin(atan([tex]\frac{\dot{y}}{Vs}[/tex])) - F_{T}= (m+vm)*[tex]\ddot{y}[/tex]

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=..., [tex]\dot{y}[/tex]= ..., [tex]\ddot{y}[/tex]..., 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!

**Physics Forums | Science Articles, Homework Help, Discussion**

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Numerical solution to a complicated ODE

**Physics Forums | Science Articles, Homework Help, Discussion**