Mesh smoothing after elements removal. Failure zone propagation

General Discussion
1

I’m simulating the hole failure under external stress state via Mohr-Coulomb failure criterion. I remove the mesh cells which cannot withstand the failure criterion. I have written a loop to observe how the failure would propagate after iterative stress redistribution around a newly formed mesh. Here how it looks like below on the first iteration. These “pinchouts“ work as stress concentrators and create enourmously large and unrealistic failure zone. Here how it looks like after first iteration (with color corresponding to failure criterion values) and the consequent second iteration of mesh truncation.

image
image912×986 27.6 KB

What I suggest here is to smooth the failure zone after truncation with some kind of approximation. Is there a way to do it in FreeFEM++ and how? Maybe there are any examples already? The code is supplied below.

//--------------------------------------------------------------------
// Dimensions
real ySize = 2. ;  // y-size of the domain
real xSize = 2. ;  // x-size of the domain
real R = 0.1 ;      // wellbore radius

// Elastic constants
real E = 1e10 ;     // Young's modulus
real nu = 0.3 ;     // Poisson's ratio

real G = E/(2*(1+nu)) ; // shear modulus
real lambda = E*nu/((1+nu)*(1-2*nu)) ; // Lame constant

//Stresses
real Sx = -80e6;
real Sy = -40e6;
real Pwell = -25e6;
real phi = 30;
real C = 10e6;

real Pi=pi;

//--------------------------------------------------------------------
// First define boundaries 
border Right(t=-ySize/2,ySize/2){x=xSize/2;y=t;}
border Top(t=xSize/2,-xSize/2){x=t;y=ySize/2;}
border Left(t=ySize/2,-ySize/2){x=-xSize/2;y=t;}
border Bottom(t=-xSize/2,xSize/2){x=t;y=-ySize/2;} 
border Well(t=0,-2*pi){x=R*cos(t);y=R*sin(t);}

//SHOW DOMAIN
plot(Right(10)+Top(10)+Left(10)+Bottom(10) + Well(40), wait=true);


//--------------------------------------------------------------------
// Create mesh 
int n = 20; // number of mesh nodes on the outer borders 
int nwell = 200; // number of mesh nodes on wellbore 
mesh Omega = buildmesh (Right(n)+Top(n)+Left(n)+Bottom(n)+Well(nwell));

plot(Omega, wait=true);

// FE spaces 
fespace Displacement(Omega, P2); 
fespace Stress(Omega, P2); 


Displacement u1, u2, v1, v2;
Stress sigmaxx, sigmayy, sigmaxy, sigma1, sigma2, tm, sm, criterion;
Stress b=0;


//--------------------------------------------------------------------
// Macro for strain 
macro e(u1,u2)
	[
		dx(u1),
		(dy(u1)+dx(u2))/2 ,
		(dx(u2)+dy(u1))/2 , 
		dy(u2)
	]//eps_xx, eps_xy , eps_yx , eps_yy
// macro for stress 
macro sigma(u1,u2)
	[
		(lambda+2.*G)*e(u1,u2)[0]+lambda*e(u1,u2)[3],
		2.*G*e(u1,u2)[1],
		2.*G*e(u1,u2)[2],
		lambda*e(u1,u2)[0]+(lambda+2.*G)*e(u1,u2)[3]
	] //stress s_xx, s_xy, s_yx, s_yy

	
	//	Define	system	of	equations 
problem	Elasticity([u1,u2],[v1,v2]) = 
	int2d(Omega)(sigma(u1,u2)'*e(v1,v2)) 
	// Boundary conditions
	+ on(Left,u1=0)              // Dirichlet boundary conditions
	+ on(Bottom,u2=0)
	+ int1d(Omega,Right)(Sx*v1)  // Neumann boundary conditions
		//- int1d(Omega,Left)(Sx*v1)
	+ int1d(Omega,Top)(Sy*v2)
		//- int1d(Omega,Bottom)(Sy*v2)
	+ int1d(Omega,Well)(Pwell*(N.x*v1+N.y*v2))
	//+ int1d(Omega,1)(Pwell*(N.x*v1+N.y*v2))
	;

problem	Elasticityupd([u1,u2],[v1,v2]) = 
	int2d(Omega)(sigma(u1,u2)'*e(v1,v2)) 
	// Boundary conditions
	+ on(Left,u1=0)              // Dirichlet boundary conditions
	+ on(Bottom,u2=0)
	+ int1d(Omega,Right)(Sx*v1)  // Neumann boundary conditions
		//- int1d(Omega,Left)(Sx*v1)
	+ int1d(Omega,Top)(Sy*v2)
		//- int1d(Omega,Bottom)(Sy*v2)
	+ int1d(Omega,Well)(Pwell*(N.x*v1+N.y*v2))
	+ int1d(Omega,5)(Pwell*(N.x*v1+N.y*v2))
	;

	//Elasticity;
//--------------------------------------------------------------------	
//	Solve system
for (int i = 0; i < 10; i++){
	Elasticity;

	// Stresses 
	sigmaxx = sigma(u1,u2)[0];
	sigmayy = sigma(u1,u2)[3]; 
	sigmaxy = sigma(u1,u2)[1];	// we could	use	[2]	as	well

	sigma1 = (sigmaxx+sigmayy)/2+sqrt(((sigmaxx-sigmayy)/2)^2+sigmaxy^2);
	sigma2 = (sigmaxx+sigmayy)/2-sqrt(((sigmaxx-sigmayy)/2)^2+sigmaxy^2);

	tm = (sigma1-sigma2)/2;
	sm = (sigma1+sigma2)/2;

	criterion = tm-sm*sin(Pi*phi/180)-C*cos(Pi*phi/180);


	//--------------------------------------------------------------------
	// plot on the deformed surface
	mesh Th=movemesh(Omega,[x+10*u1,y+10*u2]);
	plot(Th,cmm="Deformed configuration",wait=1);


	// plot the deformation field and stress
	plot([u1,u2],coef=10,cmm="Displacement field",wait=1,value=true);
	plot(sigma1,fill=1, cmm="Stress sigma1",wait=1,value=true);
	plot(sigma2,fill=1, cmm="Stress sigma2",wait=1,value=true);
	plot(criterion,fill=1, cmm="Mohr-Coulomb",wait=1,value=true);

	Omega = adaptmesh(Omega, sigma1);
	plot(Omega,cmm="Adapted",wait=1);
	Omega = trunc(Omega, criterion<0, label=4);
	plot(Omega,cmm="Destructed",wait=1);
}