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
