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]);