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

3

Hello,

I am interested in the above problem, have you found a solution for it?

4

Dear Markov,

yes I did.
I tried the Newmark algorithm from a well known reference in the field (the Bathe book “Finite Element Procedures”) but it did not work. However in the same chapter there was the description of an improvement made by Bathe, which he called “Bathe algorithm”.
It is slightly more complex as it involves computation at half time steps too, but it is much more accurate so according to the Author it proved to be useful not only for the original task (non linear stress, if I remember) but also for general problems.
Also my implementation of this algorithm was giving faulty results as diverging results were obtained.

However, after inspecting carefully the book, I discovered that there was a “step by step Example” that was supplied short afterwards, and checking carefully I found it had a different value for one constant.
So I guessed that the algorithm as described had an error (that’s strange, as I was looking in the latest edition of the book). But the example worked like a charm! And in fact it using that value for the constant the algorithm worked really well in any condition!

This is my implementation. It is a 3D fully working example of a tube excited sinusoidally at one end and fixed at the other.

load “msh3”
load “medit”
load “lapack”
load “iovtk”
load “UMFPACK64”
//load “PETSc”

real RhoTi=4430;
real ETi=114e9;
real NuTi=0.342;

real TuboL=1081.;
real TuboDi=49.;
real TuboDe=60.5;
real TuboRi=TuboDi/2;
real TuboRe=TuboDe/2;

real freqDelta=30; //hz
real freqForce=7.500E-03TuboLTuboL - 3.141E+01*TuboL + 4.651E+04;
real freqstart=freqForce-freqDelta/2.;
real freqend=freqForce+freqDelta/2.;
real toffset=-4e-3;

real DampingMass=0;
real DampingStiffness=25e-9;

int lbCera=110;
int lbTubo=803;
int lbTrasdEnd=106;

real MeshAspectRatio=3.;
real NodesPerMeter=450.;
real MagnificationOfDeformation=5.;
real Gravity=0.;
real ForceOnTip=50e5;
real dt=0.000001;
real T=dt*(18000);
string OutputFileName =“output_ThTrasdTube_” + “Di”+TuboDi1000.+“De”+TuboDe1000.+“L”+TuboL*1000. +
“dampM”+DampingMass+“dampK”+DampingStiffness+
“freqStart”+freqstart+“freqEnd”+freqend+“offset”+(-toffset)+".txt";

real work=0.;
real workmax=0.;
real defradmax=0.;
int NumberOfModesToCompute=5;

real Spessore=(TuboDe-TuboDi)/2.0;
int Fixed = 100; //Beam fixed label
int Free = 2; //Beam free label
int Forced = 300; //Beam free label
int[int] labs=[Fixed, 12, Forced, 13];

int NSteps;
NSteps=T/dt+1;

border B801(t=0,1){x=(0)+(TuboL)*t;y=(TuboRi)+(0)*t;label=801;};
border B802(t=0,1){x=(TuboL)+(0)*t;y=(TuboRi)+(5.75)*t;label=802;};
border B803(t=0,1){x=(TuboL)+(-TuboL)*t;y=(TuboRe)+(0)*t;label=803;};
border B804(t=0,1){x=(0)+(0)*t;y=(TuboRe)+(-5.75)*t;label=804;};

real nn=0.15;
real h= 1./nn;

mesh Th2=buildmesh(
B801(TuboLnn) + B803(TuboLnn) + B802(20nn) + B804(20nn)
);

int MaxLayersT=8;//(int(2piRR/h)/4)4;
func zminT = pi/2.;
func zmaxT = 0;
func fx= 0.001ysin(z);// / max( abs(cos(z) ), abs(sin(z)));
func fy= 0.001ycos(z);// / max( abs(cos(z) ), abs(sin(z)));
func fz= 0.001x;

int[int] r1T=[0,0], r2T=[0,0,2,2];
int[int] r4T=[0,2];
int labelUp=999, labelDown=998;
int[int] lbup = [0, labelUp, 6, labelUp, 8, labelUp, 17, labelUp, 23, labelUp, 25, labelUp, 28, labelUp, 35, labelUp];
int[int] lbdown = [0, labelDown, 6, labelDown, 8, labelDown, 17, labelDown, 23, labelDown, 25, labelDown, 28, labelDown, 35, labelDown];
mesh3 Th3T=buildlayers(Th2,coef= max(.1,y/min(max(20.,80-x),40.) ), MaxLayersT,zbound=[zminT,zmaxT],transfo=[fx,fy,fz],facemerge=0, labelup=lbup, labeldown=lbdown);

//
//Bethe constants
real delta = 0.5;
real alpha = 0.25 * (0.5+delta);
real alpha0=16./(dt^2);
real alpha1=4./(dt);
real alpha2=9./(dt^2);
real alpha3=3./dt;
real alpha4=2.*alpha1;
real alpha5=12./(dt^2);
real alpha6=-3./(dt^2);
real alpha7=-1./(dt);

ofstream Outputfile(OutputFileName);
Outputfile << “Lunghezza:\t” << TuboL << endl;
Outputfile << “DiametroInterno:\t” << TuboDi << endl;
Outputfile << “TuboDe:\t” << TuboDe << endl;

//Problem

func Pk = P2;
fespace Uh(Th3T, [Pk,Pk,Pk]);
fespace U0(Th3T, P0);
Uh [ux, uy, uz]=[0,0,0];

real sqrt2 = sqrt(2.);
macro Epsilon(u1,u2,u3) [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3))/sqrt2,(dz(u1)+dx(u3))/sqrt2,(dy(u1)+dx(u2))/sqrt2] // EOM
macro Divergence(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM

U0 fRho = (RhoTi);

U0 fLambda = (ETiNuTi/((1.+ NuTi)(1.-2.*NuTi)));

U0 fMu = (ETi/(2.*(1.+NuTi)));

varf vElasticity ([ux,uy,uz], [vx,vy,vz], tgv=1e18) = int3d(Th3T)( fLambda * Divergence(vx, vy, vz) * Divergence(ux, uy, uz) + 2. * fMu * ( Epsilon(vx, vy, vz)’ * Epsilon(ux, uy, uz) ) );
varf vBC ([ux,uy,uz], [vx,vy,vz], tgv=1e18) = on(801, uz=0) +on(labelUp, ux=0)+on(labelDown, uy=0) ;
varf vMass ([ux,uy,uz], [vx,vy,vz]) = int3d(Th3T)( fRho * (ux * vx + uy * vy + uz * vz) );

matrix Mass = vMass(Uh, Uh, solver=sparsesolver);
matrix Elasticity1 = vElasticity(Uh, Uh);
matrix Elasticity2 = vBC(Uh, Uh);
matrix Elasticity = Elasticity1 + Elasticity2;
matrix Damping = DampingMassMass+DampingStiffnessElasticity;

real[int] ElasticityForce(Uh.ndof);
real[int] UU(Uh.ndof), UUt(Uh.ndof), UUtt(Uh.ndof), temp(Uh.ndof);
real[int] ForceStat(Uh.ndof), ForceHalf(Uh.ndof), ForceFull(Uh.ndof), UUold(Uh.ndof), UUtold(Uh.ndof), UUttold(Uh.ndof);
real[int] UUhalf(Uh.ndof), UUthalf(Uh.ndof), UUtthalf(Uh.ndof);
real[int] UUfull(Uh.ndof), UUtfull(Uh.ndof), UUttfull(Uh.ndof);
UU=ux[]; // cause also the copy of uy
UU=0;
UUt=0;
UUtt=0;

varf vForce0 ([ux,uy,uz], [vx, vy, vz]) = int2d(Th3T,804)(vz * ForceOnTip);

// Solve static problem
set (Elasticity, solver=sparsesolver);
ElasticityForce = vForce0(0, Uh);
real[int] ElasticityBoundary = vElasticity(0, Uh);
real[int] ElasticityBC = vBC(0, Uh);
//real[int] ElasticityBF = ElasticityForce; //ElasticityBoundary + ElasticityForce;
real[int] ElasticityBF = ElasticityBoundary + ElasticityForce;
real[int] uu0 = Elasticity^-1 * ElasticityBF;
Uh [uux0, uuy0, uuz0] = uu0;
real coef=1.;

mesh3 Th0 = movemesh3(Th3T, transfo=[x+uux0coef, y+uuy0coef, z+uuz0*coef]);
mesh3 Th1;
[uux0, uuy0, uuz0] = [uux0, uuy0, uuz0];
//plot(Th0, [uux0, uuy0, uuz0],wait = true, cmm="coef amplification = "+coef);
real dxmin0 = uux0[].min;
real dymin0 = uuy0[].min;
real dzmin0 = uuz0[].min;
cout << " - dep. max x0 = “<< dxmin0 << " y0=” << dymin0 << endl;

// modified stiffness matrix for Newmark algorithm
matrix Stiffness1 = Elasticity + alpha0Mass + alpha1Damping;
matrix Stiffness2 = Elasticity + alpha2Mass + alpha3Damping;
set (Stiffness1, solver=sparsesolver);
set (Stiffness2, solver=sparsesolver);

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

//ofstream file(“f.txt”);

int fmmcyclen=0;
real fmmtipvalueold;
real fmmtipvalue;
real fmmfreq;
real fmmt0;
real pmmforce,pmmdefax,pmmdefrad;
real amplitude,phase;
real[int] maxdef(24+2);

NSteps=0;
amplitude=0;
phase=0;
cout << “time,radial_def,cera_force,cera_disp,cera_speed,cera_power,cera_work” << endl;
Outputfile << “time,radial_def,cera_force,cera_disp,cera_speed,cera_power,cera_work,trasdtip_disp”;
real zi;
int ii;
for (zi=0;zi<=TuboL;zi+=TuboL/24.)
{
Outputfile << “,def(z=” << zi << “)”;
maxdef(ii)=0;
ii++;
}
Outputfile << endl;

real cerawork=0;

for (real t=0;t<T;t+=dt)
{
macro fforce(ttt) ((-sin(2pifreqstart*(ttt+toffset))+sin(2pifreqend*(ttt+toffset)))/(pi*(freqend-freqstart)2(ttt+toffset+1e-9)))
// Bethe algorithm
varf vForceHalf ([ux,uy,uz], [vx, vy,vz]) = int2d(Th3T,804)( vz * ForceOnTip * fforce(t+dt/2)) ;
varf vForceFull ([ux,uy,uz], [vx, vy,vz]) = int2d(Th3T,804)( vz * ForceOnTip * fforce(t+dt)) ;
pmmforce=ForceOnTip * fforce(t+dt);
// computation of modified elasticity (stiffness) matrix
ForceHalf = vForceHalf(0, Uh);
temp = alpha0UU;
temp += alpha4UUt;
temp += 1.UUtt;
ForceHalf += Mass * temp;
temp = alpha1UU;
temp += 1.*UUt;
ForceHalf += Damping * temp;

UUold   = UU;
UUtold  = UUt;
UUttold = UUtt;
UUhalf  = Stiffness1^-1 * ForceHalf;
UUthalf = UUhalf-UU;
UUthalf *=alpha1;
UUthalf -=UUt;
UUtthalf= UUthalf-UUt;
UUtthalf*=alpha1;
UUtthalf-=UUtt;

// computation of second modified elasticity (stiffness) matrix 
ForceFull = vForceFull(0, Uh); 	
temp = alpha5*UUhalf;
temp += alpha6*UU;
temp += alpha1*UUthalf;
temp += alpha7*UUt;
ForceFull += Mass * temp;
temp = alpha1*UUhalf;
temp += alpha7*UU;
ForceFull += Damping * temp;

UUfull    = Stiffness2^-1 * ForceFull;
UUtfull   = -alpha7*UU;
UUtfull  -=alpha1*UUhalf;
UUtfull  +=alpha3*UUfull;

UUttfull  = -alpha7*UUt;
UUttfull -=alpha1*UUthalf;
UUttfull +=alpha3*UUtfull;

UU   = UUfull;
UUt  = UUtfull;
UUtt = UUttfull;

Uh [uux, uuy, uuz] = UU;
Uh [uutx, uuty, uutz] = UUt;

real trasdtipdisp=int2d(Th3T,lbTrasdEnd)((uux*x+uuy*y)/(x^2+y^2)^0.5);
real ceradisp=int2d(Th3T,lbCera)(uuz);
real ceraspeed=int2d(Th3T,lbCera)(uutz);
real cerapower=pmmforce*ceraspeed;
real tubodisp=int2d(Th3T,lbTubo)(abs(uux*x+uuy*y)/(x^2+y^2+1.e-5)^0.5);
cerawork+=cerapower*dt;

//work+=(pmmdefax-pmmdefaxold)*pmmforce;
//if (work>workmax) workmax=work;

cout << t+dt << "," << tubodisp << "," << pmmforce << "," << ceradisp << "," << ceraspeed << "," << cerapower << "," << cerawork << endl;
Outputfile << t+dt << "," << tubodisp << "," << pmmforce << "," << ceradisp << "," << ceraspeed << "," << cerapower << "," << cerawork << "," << trasdtipdisp;

ii=0;
for (zi=0;zi<=TuboL;zi+=TuboL/24.)
{
	Outputfile << "," << uux(TuboRe,0.,zi);
	if (uux(TuboRe,0.,zi)>maxdef(ii)) 
		maxdef(ii)=uux(TuboRe,0.,zi);
	if (maxdef(ii)>defradmax) defradmax=maxdef(ii);
	ii++;
}
Outputfile << endl;
if (1 && (NSteps % 20==0))
{
	//Movemesh 
	coef = 15000;
	Th1 = movemesh3(Th3T, transfo=[x+uux*coef, y+uuy*coef, z+uuz*coef]); 
	[uux, uuy, uuz] = [uux, uuy, uuz];
	plot(Th1, [uux, uuy, uuz], wait = false);
	int[int] fforder = [1];
	savevtk("Defo"+NSteps+".vtu", Th1, [uux, uuy, uuz], dataname = "UU", order = fforder, append = false);
	savemesh(Th1,"mesh"+NSteps+".mesh");
}
NSteps++;

}

So I struggled a lot but finally I found a version that worked well

5

Dear Riccardo,

Thank you for your detailed reply. I was also trying to implement the Newmark algorithm from Bathe’s book but I was having wrong results (maybe because on how I implemented it but I will double check). I am looking right now at the Bathe method in the book and the code that you sent me. I would like to ask you a question about the static calculation at the beginning, why exactly did you calculate it, is it a specific constraint for the problem?

Mark

6

Hello guys,

I also interested in Newmark-Beta method. I have completed the Newmark-Beta program by myself, and I am very happy to see that your program is also successful. I would like to ask if you know how to parallelize the Newmark-Beta program in FreeFem++?

I have encountered a big problem now, that is, when calculating 250,000 elements, the memory needs 160GB. Although I use a workstation to calculate, it is still too big and it takes a long time. Do you have a good way?

Best
Jack Huang

7

Hello Jack,
I have not been able to implement the Newmark Method yet, Riccardo was able to implement the Bathe method as you can see above. However, I will also try to parallelize the Newmark-Beta program on a workstation once I figure out the mistakes. I’ll let you know if I will be able to make it. On a side note do you have any references for the Newmark-Beta algorithm that you implemented?

Best
Mark