Hello everyone,
I am working on a simulation regarding the compression of a pillar ( Rectangular form).
At high level of deformation, two elements touch each other and my calculation get stopped(see the figure attached).
Please any suggestion is very welcome.
Hello Jalal,
It is not clear what your picture means. What do the colors represent? I understand that the mesh is deformed and that your problem is its boundary.
We had a similar problem in
J. Comput. Phys., 333:387–408, 2017. DOI: 10.1016/j.jcp.2016.12.036
https://perso.math.u-pem.fr/bouchut.francois/publications/drucker-prager2D.pdf
We managed to solve it by putting some numerical surface tension, see in particular figure 4 p. 11, to prevent the folding of the free surface over itself.
This surface tension is the last term in the right-hand side of the variational formulation 2.17a
Thank you François for your help.
I found the paper very interesting and aligned perfectly with what I was searching for.
I am trying to implement the numerical tension.
The simulation is a compression test of a multilayer pillar, made of alternation of two compositions.
The colors in the picture (Blue and Orange) correspond to the yield stength values of the two materials (Compositions). I use this representation to identify the layers during the simulation animation’s.
It is a coarse mesh indeed but I have the same problem when the mesh is fine.
You make a nice computation!
In our code I remember that obtaining the local curvature of the boundary was a bit messy.
Here is the piece of code (sorry comments are in french). It is not optimal:
the construction of R is repeated on 3 pieces of the boundary
real[int] rx(Th.nbe),ry(Th.nbe),curvx(Th.nbe), curvy(Th.nbe);
int[int] ibe(Th.nbe);
real curvx0,curvy0;
int nbv=0; // compte le nombre de noeud au bord
// les noeuds du bord sont ranges dans l'ordre: label=1,2,3,4
// surface libre parcourue de droite (2) à gauche (3)
for(int n=0;n<Th.nbe;n++){
if(Th.be(n).label==2){
//-------------------- construction R(x,y) ----------------------
rx(nbv)=Th.be(n)[0].x-Th.be(n)[1].x;
ry(nbv)=Th.be(n)[0].y-Th.be(n)[1].y;
curvx(nbv)=Th.be(n)[1].x;
curvy(nbv)=Th.be(n)[1].y;
ibe(nbv)=n;
if(nbv==0){
curvx0=Th.be(n)[0].x;
curvy0=Th.be(n)[0].y;
}
nbv++;
}
}
int ncrit=nbv-1;
for(int n=0;n<Th.nbe;n++){
if(Th.be(n).label==31){
//-------------------- construction R(x,y) ----------------------
rx(nbv)=Th.be(n)[0].x-Th.be(n)[1].x;
ry(nbv)=Th.be(n)[0].y-Th.be(n)[1].y;
curvx(nbv)=Th.be(n)[1].x;
curvy(nbv)=Th.be(n)[1].y;
ibe(nbv)=n;
if(nbv==0){
curvx0=Th.be(n)[0].x;
curvy0=Th.be(n)[0].y;
}
nbv++;
}
}
for(int n=0;n<Th.nbe;n++){
if(Th.be(n).label==32){
//-------------------- construction R(x,y) ----------------------
rx(nbv)=Th.be(n)[0].x-Th.be(n)[1].x;
ry(nbv)=Th.be(n)[0].y-Th.be(n)[1].y;
curvx(nbv)=Th.be(n)[1].x;
curvy(nbv)=Th.be(n)[1].y;
ibe(nbv)=n;
if(nbv==0){
curvx0=Th.be(n)[0].x;
curvy0=Th.be(n)[0].y;
}
nbv++;
}
}
rx.resize(nbv); ry.resize(nbv);
curvx.resize(nbv); curvy.resize(nbv);
ibe.resize(nbv);
real[int] arctan(nbv+1),curv(nbv+1),num(nbv+1),denom(nbv+1),lg(nbv+1);
real normr,normm;
curv(0)=0;
for(int n=0;n<nbv-1;n++){
num(n+1)=-ry(n+1)*rx(n)+rx(n+1)*ry(n); // rt(n).r(n-1) avec rt(n)=(-ry,rx)
normr=sqrt(rx(n+1)^2+ry(n+1)^2);
normm=sqrt(rx(n)^2+ry(n)^2);
denom(n+1)=normm*normr+rx(n+1)*rx(n)+ry(n+1)*ry(n);
arctan(n+1)=2*atan(num(n+1)/denom(n+1));
lg(n+1)=0.5*(normr+normm);
curv(n+1)=arctan(n+1)/lg(n+1);
//cout << "curv("<<n<<")="<<curv(n)<<endl;
}
curv(nbv)=curv(nbv-1);
//-------------------- COURBURE ---------------------
real h=0.025, cst=0., deltan=16.;
fcurv[]=0.; fcurv=0.;
fcurv[][int(Th.be(ibe(0))[0])]=curv(0);
for(int n=0;n<nbv;n++)
fcurv[][int(Th.be(ibe(n))[1])]=h*curv(n+1)*max(0.,1.-((n-ncrit)/deltan)^2)^2;
{
ofstream fichier("curvature.dat");
for(int n=0;n<nbv;n++){
fichier << curv(n) << " " << ibe(n) << " " << curvx(n) << " " << curvy(n) << " " << num(n) << " " << denom(n)<< " " << arctan(n)<< " " << lg(n)<<"\n";
}
}
full code
drucker-prager_ale.edp (9.9 KB)