Finite Element Method in Aeroelasticity (Flutter)

General Discussion
1

Dear Frederic

I’m working on a fluid structure interaction problem , its an airfoil on air flow to simulate the flutter phenomenon.

The problem consist on solve the Navier Stokes equations couple to system of nonlinear ordinary differential equations of secod order that describes de airfoil motion, the airfoil have two degrees of freedom, it can be move in the vertical direction and rotate around the elastic axis.

Request:
I´m trying to make that the airfoil can rotate around the elastic axis using the rotational matrix (Euler Angles) and it can move in vertical direction in order to couple the ODE. The idea is to use the option movemesh in Freefem.

image
image
image

image

In this equations L =Lift, M= Moment and EO is the elastic axis.

This is the code, if you can explain how to do this step, I really appreciate.

// Air - Airfoil Structure Interaction with Navier-Stokes equations
// Outer Rectangle [A,B]X[C,D]
// Obstacle Airfoil.

// Variables
real L = -0.5;
real R = 5.0;
real B = -1.0;
real T = 1.0;
int Px = 40;
real Py = 28;
real S = 20;
real m = 0.086622; // Mass kg//
real Ialpha = 0.000487291; // Moment of inertia kg.m^2 //
real Salpha = -0.000779673; // Moment of inertia kg.m //
real khh = 105.109; // Rigidez N/m //
real kalp = 3.695582; // Rigidez Nm/rad //
real l = 0.05; // Spam m //
real c = 0.3; // Chord m//
real rho =1.225; // Density Kg/m^3//
real v = 0.000015; // Viscosity m/s^2//
real XEO = 0.12; // Elastic axis 0.4c //
real xt = 0.111; // Center of gravity 0.37c//
real dhh = 0.105109; // Proportional Damping k_hhe with e=10^-3 //
real daplh = 0.003669558; // Proportional Damping k_alpe with e=10^-3 //
real hh = 1; // Step of ODE
real M = 1./(mIalpha- Salpha^2(cos(z))^2);
real A = cos(z);
real O = sin(z);

// Boundary Conditions
real hig = -0.05;
real dhig = 0;
real alp = 6;
real dalp = 0;

// ODE derivatives Angle and high.
real ddhig; // Second derivative of hight
real ddalp; // Second derivative of angle

//Airfoil
border upp(t=0, 1){x=t; y=0.17735sqrt(t) - 0.075597t - 0.212836*(t^2) + 0.17363*(t^3) - 0.06254*(t^4); label =5;}
border low(t=1, 0){x = t; y=-(0.17735sqrt(t) -0.075597t - 0.212836*(t^2) + 0.17363*(t^3) - 0.06254*(t^4));label=5;}

// Outer Rectangle
border Bottom(t= 0,1){x= L + (R-L)*t; y=B; label= 1;};
border Top(t= 1,0){x= L + (R-L)*t; y= T; label= 3;};
border Left(t= 1,0){x= L; y= B + (T-B)*t; label= 4;};
border Right(t= 0,1){x= R; y= B + (T-B)*t; label= 2;};

//Mesh Domain
plot( Bottom(Px) + Right(Py) + Top(-Px) + Left(-Py)+ upp(-S) + low(S), wait=1, ps=“geomoutline.eps” );
mesh Th= buildmesh( Bottom(Px) + Right(Py) + Top(Px) + Left(Py)+ upp(S) + low(S));
Th = adaptmesh(Th, 1./30., IsMetric=1, nbvx=10000);
plot(Th, wait=true, ps=“NACA0012.eps”, bw=true);
macro airfoilrot( ;
mesh th= movemesh(Th, [x+1, y+2]);

//Definition of Space
fespace Vh1( Th, P2);
fespace Vh2( Th, P1);
Vh1 u2, v2, u1, v1, uu1, uu2;
Vh2 p, q;

int numTSteps=80;
bool reuseMatrix=false;
real nu=1./1000.;
real dt=0.05;
real bndryVelocity;

// Functions of ODE
func real fun1(real hig, real dhig, real alp, real dalp)
{
real f1= dhig;
return f1;
}
func real fun2( real hig, real dhig, real alp, real dalp)
{
real f2 = (-M* Ialpha * dhh)dhig - (-MOIalpha - MSalphaA)dalp + (MIalphakhh)hig +(MSalphaA kalp)*alp;
return f2;
}
func real fun3(real hig, real dhig, real alp, real dalp)
{
real f3= dalp;
return f3;
}

func real fun4( real hig, real dhig, real alp, real dalp)
{
real f4 = (M* SalphaAdhh)dhig + ((MOSalpha)-(Mmdaplh))dalp + (MSalphaAkhh)hig + (Mkalpm)*alp ;
return f4;
}

// Weak Formulation
problem NS ([u1, u2, p] , [v1, v2, q] , solver=UMFPACK , init=reuseMatrix) =
int2d(Th)( 1./dt*( u1v1 + u2v2 )
+ nu * ( dx(u1)dx(v1) + dy(u1)dy(v1)
+ dx(u2)dx(v2) + dy(u2)dy(v2) )
+ pq(0.000001)
+ v1dx§ + v2dy§
+ dx(u1)*q + dy(u2)*q )

  • int2d(Th) ( 1./dt*(convect( [uu1, uu2] , -dt , uu1) * v1
    -convect( [uu1, uu2] , -dt , uu2) * v2 ))

    // b.c.: uniform velocity top, bottom, inlet (left)
    // “do nothing” on exit (right)

  • on(1, u1=bndryVelocity, u2=0)

  • on(3, u1=bndryVelocity, u2=0)

  • on(4, u1=bndryVelocity, u2=0)

  • on(5, u1=0, u2=0)
    ;

real[int] times(numTSteps);
real var(numTSteps);

// Iterative Process
for (int i=0 ; i < numTSteps ; i++) {
times[i] = i*dt;

// ramp up velocity from 0.0
bndryVelocity = i*dt;
if (bndryVelocity >= 0.1) bndryVelocity = 0.1;

uu1 = u1;
uu2 = u2;

// Solve the problem
NS;

// reuse the matrix in the rest of the iterations
reuseMatrix = true;

// plot the current solution
plot(coef=0.2, cmm= " [u1,u2] and p ", p, [u1,u2],
ArrowShape= 1, ArrowSize= -0.8 );

// System 4x4 ODE Euler Method

for (hh =0; hh < 50; hh = hh+1){
dhig = dhig + fun1(hig, dhig, alp, dalp)*hh;
dhig = dhig;
ddhig = ddhig + fun2(hig, dhig, alp, dalp)*hh;
ddhig = ddhig;
dalp = dalp + fun3(hig, dhig, alp, dalp)*hh;
dalp = dalp;
ddalp = ddalp + fun4(hig, dhig, alp, dalp)*hh;
ddalp = ddalp;

cout << "ddhig = " << ddhig << endl;
cout << "ddalp = " << ddalp << endl;

}
}