Transient problems in elasticity

General Discussion
#1

Dear all,

first, many thanks for your wonderful software.
I’m trying to solve transient problems in elasticity, more in particular simulation of mechanical resonators which have some damping (e.g. the time dependent response of a bell that is stroke at part of the boundary at time t=0, or simulation of ultrasonic transducers).
I read the manual and looked for examples and also in the community blog but I wasn’t able to find any example of solutions of transient problems in elasticity.

I know that usually there are two ways to perform such task:

  • direct time integration
  • modal superposition

I guess that Freefem++ is able to solve such problems. Do you have any example or documentation which could help implementing such type of problems in Freefem++?
Thanks a lot!!

Riccardo

#2

Dear all,

I tried to build a code to solve transient problems in elasticity (a simple bar with a fixed end, gravity force and a periodic force on the tip) but either it gives nonsense results or freefem crashes.
The time integration algorithm is Newmark (see e.g. Bathe, Finite Element Procedures, paragraph 9.2.3 “The Newmark method”).

I’m struggling to find any working solution.
Please could you help me?
Below the code.

Thanks a lot
Riccardo

load “lapack”

//Parameters
real dt=0.00001;
real T=0.001;

real Rho = 8000.; //Density
real E = 210.e9; //Young modulus
real Nu = 0.27; //Poisson ratio
real Gravity = -9.81; //Gravity
real ForceOnTip = 1e-10;

//Newmark constants
real delta = 0.5;
real alpha = 0.25 * (0.5+delta);
real alpha0=1/(alphadt^2);
real alpha1=delta/(alphadt);
real alpha2=1/(alphadt);
real alpha3=1/(2.alpha)-1.;
real alpha4=delta/alpha-1.;
real alpha5=dt/2.(delta/alpha-2.);
real alpha6=dt(1-delta);
real alpha7=delta*dt;

//Mesh quality
real L = 20.; //Beam length
real H = 1.; //Beam height
int Fixed = 1; //Beam fixed label
int Free = 2; //Beam free label
int Forced = 3; //Beam free label
border b1(t=0., L){x=t; y=0.; label=Free;};
border b2(t=0., H){x=L; y=t; label=Forced;};
border b3(t=0., L){x=L-t; y=H; label=Free;};
border b4(t=0., H){x=0.; y=H-t; label=Fixed;};

int nnL = 200; //max(2., nnL);
int nnH = 10; //max(2., nnH);
mesh Th = buildmesh(b1(nnL) + b2(nnH) + b3(nnL) + b4(nnH));

//Fespace
func Pk = P2;
fespace Uh(Th, [Pk, Pk]);
Uh [ux, uy]=[0,0];
real[int] UU(Uh.ndof), UUt(Uh.ndof), UUtt(Uh.ndof), temp(Uh.ndof);
real[int] ElasticityForce(Uh.ndof), UUold(Uh.ndof), UUttold(Uh.ndof);
UU=ux[]; // cause also the copy of uy
UUt=0;

//Macro
real sqrt2 = sqrt(2.);

macro Epsilon(ux, uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt2] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //

//Problem
real Mu = E/(2.(1.+Nu));
real Lambda = ENu/((1.+ Nu)*(1.-2.*Nu));
real freq = 20000;

varf vElasticity ([ux,uy], [vx, vy]) = int2d(Th)( Lambda * Divergence(vx, vy) * Divergence(ux, uy) + 2. * Mu * ( Epsilon(vx, vy)’ * Epsilon(ux, uy) ) );
varf vBC ([ux,uy], [vx, vy]) = on(Fixed, ux=0, uy=0) ;
varf vMass ([ux,uy], [vx, vy]) = int2d(Th)( Rho * (ux * vx + uy * vy)/(dt*dt) );

matrix Mass = vMass(Uh, Uh, solver=sparsesolver);
matrix Elasticity1 = vElasticity(Uh, Uh);
matrix Elasticity2 = vBC(Uh, Uh);
matrix Elasticity = Elasticity1 + Elasticity2; // I’m not sure if boundary conditions should be applied to Elasticity matrix
matrix Damping = 0Mass+0Elasticity;

// modified stiffness matrix for Newmark algorithm
matrix ElasticityStiffness = Elasticity + alpha0Mass + alpha1Damping;

// Solve static problem
set (Elasticity, solver=sparsesolver);
real[int] ElasticityBoundary = vElasticity(0, Uh);
real[int] ElasticityBC = vBC(0, Uh);
real[int] ElasticityBF = ElasticityForce; //ElasticityBoundary + ElasticityForce;
real[int] uyy = Elasticity^-1 * ElasticityBF;

// Prepare dynamic (transient problem)
// computation of acceleration (UUtt)
temp = DampingUUt;
temp += ElasticityUU;
temp = ElasticityForce - temp;
UUtt = Mass^-1 * temp;

for (real t=0;t<T;t+=dt)
{
// Newmark algorithm
varf vForce ([ux,uy], [vx, vy]) = int2d(Th)( Rho * Gravity * vy ) + on(Forced, ux=0, uy=ForceOnTipsin(23.14159freq(t+dt)));

// computation of modified elasticity (stiffness) matrix 
ElasticityForce = vForce(0, Uh); 	
temp = alpha0*UU;
temp += alpha2*UUt;
temp += alpha3*UUtt;
ElasticityForce += Mass * temp;
temp = alpha1*UU;
temp += alpha4*UUt;
temp += alpha5*UUtt;
ElasticityForce += Damping * temp;

UUold   = UU;
UUttold = UUtt;
UU   = ElasticityStiffness^-1 * ElasticityForce;

temp = UU-UUold;
UUtt = alpha0 * temp;
temp = alpha2*UUt;
temp += alpha3*UUttold;
UUtt -= temp;

UUt  += alpha6*UUttold;
UUt  += alpha7*UUtt;

}
Uh [uux, uuy] = UU;
real dxmin = uux[].min;
real dymin = uuy[].min;
cout << " - dep. max x = “<< dxmin << " y=” << dymin << endl;
cout << " dep. (20, 0) = " << uux(20, 0) << " " << uuy(20, 0) << endl;

//[ux, uy] = UU;
//Movemesh
real coef = 2e11;
mesh Th1 = movemesh(Th, [x+uuxcoef, y+uuycoef]);

[uux, uuy] = [uux, uuy];
plot(Th1, [uux, uuy]);