How to use the Newmark algorithm to solve transient dynamic problems in FreeFEM?

General Discussion
1

Hello everyone,

I’m currently working on a transient dynamic problem using the Newmark time integration scheme in FreeFEM. While I have carefully implemented the algorithm and double-checked my code multiple times, the simulation still fails to produce correct results. The solution appears to diverge or behave unphysically after a few time steps.

If anyone has experience with transient simulations using Newmark in FreeFEM or has encountered similar issues, I would really appreciate your insights or suggestions.

The detail of the code is as following:

load "iovtk"
//properties of material
real E = 3.2e10, nu = 0.2, rho=2450;
//damping parameters
real am = 0e-5, ak = 0e-5;
//Newmark algorithm parameters
real deltat = 1e-7, delta = 1.0/4.0, gamma = 1.0/2.0;
real b0 = 1/delta/deltat^2, b1 = gamma/delta/deltat, b2 = 1/delta/deltat, b3 = 1/2/delta-1, 
     b4 = gamma/delta-1, b5 = deltat/2*(gamma/delta-2), b6 = deltat*(1-gamma), b7 = gamma*deltat;
//solver parameters
real loadinit=0e-1, loadrate = 1e3, damagethreshold=1e-2, maxIterations=1000, pressure=1e6;

//define stiffness matrix
macro Q(E,nu)  [ [E/(1-nu^2),nu*E/(1-nu^2),0],
                [nu*E/(1-nu^2), E/(1-nu^2),0],
                [0,0,0.5*E/(1+nu)] ]//
//define epsilon
macro epsilon(u,v) [dx(u),dy(v),dy(u)+dx(v)]//

//load mesh
border c1(t=0, 1){x=t; y=0.; label=1;};
border c2(t=0, 1){x=1; y=t; label=2;};
border c3(t=0, 1){x=1-t; y=1; label=3;};
border c4(t=0, 1){x=0; y=1-t; label=4;};
mesh Th = buildmesh(c1(20) + c2(20) + c3(20) + c4(20));
plot(Th,wait=true);

//define space
fespace Vh(Th,[P1,P1]), Eh(Th,P1);
Vh [u,v],[u1,v1],[vv1,vv2],[aa1,aa2],[tempu,tempv],[tempvv1,tempvv2],[tempaa1,tempaa2];
Eh uu, vv; 

//degrees of freedom
int DofP1 = Vh.ndof;
int DofP0 = Eh.ndof;

//difine varf
varf Mass([u,v],[u1,v1]) = int2d(Th)(rho*[u,v]'*[u1,v1])+on(1,u=0,v=0);
varf Stiff([u,v],[u1,v1]) = int2d(Th)(epsilon(u,v)'*Q(E,nu)*epsilon(u1,v1))+on(1,u=0,v=0);
varf Tensile([u,v],[u1,v1]) = int1d(Th,3)(pressure*v1);
//varf Tensile([u,v],[u1,v1]) = on(1,v=-displacement)+on(2,v=displacement);

//initiation
[u,v]=[0,0]; [vv1,vv2]=[0,0]; [tempu,tempv]=[0,0]; [tempvv1,tempvv2]=[0,0];
matrix Muu = Mass(Vh,Vh,solver = sparsesolver);
matrix Kuu = Stiff(Vh,Vh,solver = sparsesolver);
real[int] bu = Tensile(0,Vh);
real[int] vec0(DofP1);
real[int] tempvec = Kuu*u[];
for(int j = 0; j < DofP1; j++)
{
    vec0(j) = bu[j]-tempvec(j);
}
aa1[] = Muu^-1*vec0;
[tempaa1,tempaa2]=[aa1,aa2];
//plot(aa2,fill=1,value=1,wait=1,cmm="dis"+" aa2"+0);

//save data
ofstream Fx("output/displacement-force.txt");
Fx<<"displacement-force"<<endl;
int[int] vtuOrder=[1,1];
uu=0; vv=0;
//savevtk("output/phasefield/damage000.vtu", Th, uu, vv, dataname="u v", order=vtuOrder);

//time step solver
int index = 0;
real displacement,Force;
displacement=loadinit;
while (index<maxIterations-1)
{
    real damageError = 1;
    int iter=0;
    displacement = displacement+loadrate*deltat;
    while (damageError > damagethreshold)
    {
        //solve displacement
        Muu = Mass(Vh,Vh, solver = sparsesolver);
        Kuu = Stiff(Vh,Vh, solver = sparsesolver);
        matrix Cuu = am*Muu + ak*Kuu;
        matrix Keq = Kuu + b0*Muu + b1*Cuu;
        set(Keq, solver = sparsesolver);
        bu = Tensile(0,Vh);
        real[int] beq(DofP1);
        //equivalent load
        real[int] vec1(DofP1), vec2(DofP1);
        for(int j = 0; j < DofP1; j++)
        {
            vec1(j) = b0*tempu[](j)+b2*tempvv1[](j)+b3*tempaa1[](j);
            vec2(j) = b1*tempu[](j)+b4*tempvv1[](j)+b5*tempaa1[](j);
        }
        real[int] vec3 = Muu*vec1, vec4 = Cuu*vec2;
        for(int i = 0; i < DofP1; i++)
        {
            beq(i) = bu(i) + vec3(i) + vec4(i);
        }
        u[] = Keq^-1*beq;
        //u[] = Kuu^-1*bu;
        damageError = 0;
    }
    //update velosity and accelerate
    for(int i = 0; i < DofP1; i++)
    {
        aa1[](i) = b0*(u[](i)-tempu[](i)) - b2*tempvv1[](i) - b3*tempaa1[](i);
        vv1[](i) = tempvv1[](i) + b6*tempaa1[](i) + b7*aa1[](i);
    }
    [tempu,tempv] = [u,v]; [tempvv1,tempvv2] = [vv1,vv2]; [tempaa1, tempaa2] = [aa1,aa2];
    plot(v,fill=1,value=1,wait=1,cmm="dis"+" displacement"+displacement);
    //save data
    index = index+1;
    Fx<<displacement<<"  "<<Force<<endl;
    uu=u; vv=v;
    //savevtk("output/phasefield/damage00"+index+".vtu", Th, uu, vv, dataname="u v", order=vtuOrder);
}
2

In FreeFEM, using the central difference method seems to produce stable results. Strangely, however, the unconditionally stable Newmark method yields a divergent solution.

3

I have found the mistake, which is related to the variable b3 on line 8. When I tried to output the value of b3, I found it equals to -1, but it should be 1. So I corrected it as b3=1.0/2.0/delta-1.0, and got the right result. I think b3 follows the type of the first number behind the equal sign, which made it be ’int‘ type but not ‘real’.

4

Hello,
the variable b3 is really declared as real, but when you write b3 = 1/2/delta-1, FF computes in the order of appearence 1/2 then /delta, then -1. Finally the result is affected to b3.
When you have two numbers of the same type (real or integer), the result of their addition or multiplication is of the same type. If they are not the same type (one is real and the other is integer), the integer is converted to real before performing the computation, and the result is real.
With these rules you understand the result. It is because the first operation 1/2 is performed as integers, and the result is 0 integer.
Another obvious rule is that if you declare b3 real and you affect to it an integer value, it is of course converted to real.
Similarly if you declare int a; a=1.4; the real value 1.4 is converted to int, hence you get a=1.

In your case the correct and safe writing is as you say b3=1./2./delta-1. so that everything is real and computed as real.

By the way, do you have some reference for the proof of stability of the Newmark algorithm? Does it work for nonsymmetric matrices, or do you need symmetric matrices?