Calculation problem in J-integral

1

JMesh
JMesh1228×1216 118 KB

Problem description
Problem description1550×858 79.3 KB

I have been struggling with a problem for some time, which is to calculate the J-integral using the phase field fracture model. With the guidance of the program developed by Emilio Mart ́ınez-Pan ̃eda and his colleagues, I was able to perform the computations on the model. Through testing, it was found that consistent loading curves could be obtained. Specifically, when ud approaches 0.0059, rapid fracture behavior of the structure occurs.

It is well-known that fracture occurs when J equals the critical value of Gc for linear elastic materials. Therefore, I believed that using the standard definition of J-integral should yield a J value equivalent to G when ud approaches 0.0059 (where Gc is defined as 2.7). However, I encountered an unexpected result while using my own J-integral program for computation. Specifically, when ud=0.0059, the calculated J value was 4.3, which greatly exceeds Gc. This calculation result is clearly unreasonable, but I have yet to identify a viable solution to this issue.

The mesh partition used in my analysis is illustrated in Figure 1, where several rectangular enclosures are employed to compute the J-integral. The program utilized is provided in the attached file.

You can browse either J_Caculation.py or J_Caculation.ipynb, they are completely equivalent.

2

It seems that various types of files cannot be uploaded, so here is the relevant code.

J_Mesh.geo
//+

mesh1 = 0.05;
mesh2 = 0.008;

y = 1e-3;
l0 = 0.02;
y0 = 0.1;//积分点的选择

Point(1) = {0, 0, 0, mesh2};
Point(2) = {0.5, 0.5, 0, mesh1};
Point(3) = {-0.5, 0.5, 0, mesh1};
Point(4) = {-0.5, -0.5, 0, mesh1};
Point(5) = {0.5, -0.5, 0, mesh1};

Point(6) = {0.5,0.15,0,mesh2};
Point(7) = {0.5,-0.15,0,mesh2};
Point(8) = {-0.5,0.15,0,mesh2};
Point(9) = {-0.5,-0.15,0,mesh2};

Point(10) = {-0.4,0.8y,0,mesh2};
Point(11) = {0.4,0.12,0,mesh2};
Point(12) = {0.4,-0.12,0,mesh2};
Point(13) = {-0.4,-0.8y,0,mesh2};
Point(14) = {-0.4,0,0,mesh2};

Point(15) = {-0.2,0.4y,0,mesh2};
Point(16) = {0.2,0.06,0,mesh2};
Point(17) = {0.2,-0.06,0,mesh2};
Point(18) = {-0.2,-0.4y,0,mesh2};
Point(19) = {-0.2,-0,0,mesh2};

Point(20) = {0.25,0.04,0,mesh2};
Point(21) = {0.25,-0.04,0,mesh2};
Point(22) = {0.35,0.04,0,mesh2};
Point(23) = {0.35,-0.04,0,mesh2};
Point(24) = {0.25,-0,0,mesh2};

Point(25) = {-0.5,y,0,mesh2};
Point(26) = {-0.5,-y,0,mesh2};

Point(27) = {-0.4,-0.12,0,mesh2};
Point(28) = {-0.4,0.12,0,mesh2};

Point(29) = {-0.2,-0.06,0,mesh2};
Point(30) = {-0.2,0.06,0,mesh2};

Point(31) = {-0.025,0.05y,0,mesh2};
Point(32) = {-0.025,-0.05y,0,mesh2};
Point(33) = {0.025,0.05y,0,mesh2};
Point(34) = {0.025,-0.05y,0,mesh2};

Point(35) = {-0.025,-0.02,0,mesh2};
Point(36) = {-0.025,0.02,0,mesh2};

Point(37) = {0.025,0.02,0,mesh2};
Point(38) = {0.025,-0.02,0,mesh2};

//+
Line(1) = {2, 3};
//+
Line(2) = {3, 8};
//+
Line(3) = {8, 6};
//+
Line(4) = {6, 2};
//+
Line(5) = {7, 9};
//+
Line(6) = {9, 4};
//+
Line(7) = {4, 5};
//+
Line(8) = {5, 7};
//+
Line(9) = {8, 25};
//+
Line(10) = {25, 10};
//+
//+
Line(11) = {10, 15};
//+
Line(12) = {15, 31};
//+
Line(13) = {31, 1};
//+
Line(14) = {1, 32};
//+
Line(15) = {32, 18};
//+
Line(16) = {18, 13};
//+
Line(17) = {13, 26};
//+
Line(18) = {26, 9};
//+
Line(19) = {6, 7};
//+
Line(20) = {13, 27};
//+
Line(21) = {27, 12};
//+
Line(22) = {12, 11};
//+
Line(23) = {11, 28};
//+
Line(24) = {28, 10};
//+

//+
Line(26) = {18, 29};
//+
Line(27) = {29, 17};
//+
Line(28) = {17, 16};
//+
Line(29) = {16, 30};
//+
Line(30) = {30, 15};
//+
Line(31) = {32, 35};
//+
Line(32) = {35, 38};
//+
Line(33) = {38, 37};
//+
Line(34) = {37, 36};
//+
Line(35) = {36, 31};
//+
Line(36) = {24, 21};
//+
Line(37) = {21, 23};
//+
Line(38) = {23, 22};
//+
Line(39) = {22, 20};
//+
Line(40) = {20, 24};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {5, 6, 7, 8};
//+
Plane Surface(2) = {2};
//+
Curve Loop(3) = {3, 19, 5, -18, -17, 20, 21, 22, 23, 24, -10, -9};
//+
Plane Surface(3) = {3};
//+
Curve Loop(4) = {24, 11, -30, -29, -28, -27, -26, 16, 20, 21, 22, 23};
//+
Curve Loop(5) = {39, 40, 36, 37, 38};
//+
Plane Surface(4) = {4, 5};
//+
Curve Loop(6) = {29, 30, 12, -35, -34, -33, -32, -31, 15, 26, 27, 28};
//+
Plane Surface(6) = {6};
//+
Plane Surface(5) = {5};
//+
Curve Loop(7) = {34, 35, 13, 14, 31, 32, 33};
//+
Plane Surface(7) = {7};

Physical Surface(1) = {1};
Physical Surface(2) = {2};
Physical Surface(3) = {3};
Physical Surface(4) = {4};
Physical Surface(5) = {5};
Physical Surface(6) = {6};
Physical Surface(7) = {7};

Physical Line(11) = {20};
Physical Line(12) = {21};
Physical Line(13) = {22};
Physical Line(14) = {23};
Physical Line(15) = {24};

Physical Line(21) = {26};
Physical Line(22) = {27};
Physical Line(23) = {28};
Physical Line(24) = {29};
Physical Line(25) = {30};

Physical Line(31) = {31};
Physical Line(32) = {32};
Physical Line(33) = {33};
Physical Line(34) = {34};
Physical Line(35) = {35};

Physical Line(41) = {36};
Physical Line(42) = {37};
Physical Line(43) = {38};
Physical Line(44) = {39};
Physical Line(45) = {40};

3

And Calculation code:

import meshio
import dolfin as dl
import matplotlib.pyplot as plt
import numpy as np

Convert Mesh

mesh_from_file = meshio.read(“J_Mesh.msh”)
def create_mesh(mesh, cell_type, prune_z=False):
cells = np.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])
cell_data = np.hstack([mesh.cell_data_dict[“gmsh:physical”][key]
for key in mesh.cell_data_dict[“gmsh:physical”].keys() if key==cell_type])
# Remove z-coordinates from mesh if we have a 2D cell and all points have the same third coordinate
points= mesh.points
if prune_z:
points = points[:,:2]
mesh_new = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={“name_to_read”:[cell_data]})
return mesh_new

line_mesh = create_mesh(mesh_from_file, “line”, prune_z=True)
meshio.write(“facet_mesh.xdmf”, line_mesh)
triangle_mesh = create_mesh(mesh_from_file, “triangle”, prune_z=True)
meshio.write(“mesh.xdmf”, triangle_mesh)

mesh = dl.Mesh()
with dl.XDMFFile(“mesh.xdmf”) as infile:
infile.read(mesh)
mvc = dl.MeshValueCollection(“size_t”, mesh, 1)
with dl.XDMFFile(“facet_mesh.xdmf”) as infile:
infile.read(mvc, “name_to_read”)
mf = dl.cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_J = dl.Measure(“dS”, domain=mesh, subdomain_data=mf)#边界的识别

dl.plot(mesh,title=‘Mesh’);

lamda = 121.1538e3;
mu = 80.7692e3;
kappa = lamda+2/3*mu;
Gc = 2.7;
l0 = 0.02;
L = 1;

def epsilon(u):
return dl.sym(dl.grad(u))

def sigma(u):
eps = epsilon(u);
dim = 2;
return 2mueps + lamda*dl.tr(eps)*dl.Identity(dim)

def sigma_damage(u,phi):
eps = epsilon(u);
dim = 2;
sigma_undamage = 2mueps + lamda*dl.tr(eps)dl.Identity(dim); #无损伤应力
omega = (1-phi)**2;
sigma_damage = omegasigma_undamage;
return sigma_damage

def psi(u):
eps = epsilon(u);
eps_D = dl.dev(eps);
eps_S = eps-eps_D;
#psi+:
psi1 = 1/2kappa(1/2*(dl.tr(eps)+abs(dl.tr(eps))))*2 + mudl.inner(eps_D,eps_D);
#psi-:
psi2 = 1/2kappa(dl.tr(eps)-abs(dl.tr(eps)))**2;

return psi1 #return psi+

def H(u_old,u_new,H_old):

H_max = dl.conditional(psi(u_old) < psi(u_new),psi(u_new),H_old);
return H_max

class Top(dl.SubDomain):
def inside(self, x, on_boundary):
return dl.near(x[1],L/2) and on_boundary

class Bottom(dl.SubDomain):
def inside(self, x, on_boundary):
return dl.near(x[1],-L/2) and on_boundary

class Crack(dl.SubDomain):
def inside(self, x, on_boundary):
return dl.near(x[1],0,1e-3) and (x[0]<=0)

def Midpoint(x, on_boundary):
tol = 1e-3
return (abs(x[1] -(-L/2)) < tol) and (abs(x[0] -(0)) < tol)

V_u = dl.VectorFunctionSpace(mesh,‘CG’,1);
V_phi = dl.FunctionSpace(mesh,‘CG’,1);
V_H = dl.FunctionSpace(mesh,‘CG’,1);
V_psi = dl.FunctionSpace(mesh,‘CG’,1);

u = dl.TrialFunction(V_u);
v = dl.TestFunction(V_u);
phi = dl.TrialFunction(V_phi);
theta = dl.TestFunction(V_phi);

Mark

boundaries = dl.MeshFunction(“size_t”, mesh, 1)
boundaries.set_all(0)
Top().mark(boundaries, 1)
Bottom().mark(boundaries, 2)
Crack().mark(boundaries, 3)
ds = dl.Measure(‘ds’, domain=mesh, subdomain_data=boundaries) #手册P96

ud = dl.Expression(“ud”, ud = 0.00, degree=1)

Boundary Condition

BC_u = [dl.DirichletBC(V_u,dl.Constant((0,0)), boundaries, 2),
dl.DirichletBC(V_u.sub(1), ud, boundaries, 1)];
BC_phi = [dl.DirichletBC(V_phi, dl.Constant(1.), boundaries, 3)];

u_old,u_new = dl.Function(V_u),dl.Function(V_u);
phi_old,phi_new = dl.Function(V_phi),dl.Function(V_phi);
H_old = dl.Function(V_H);

Equations

Equation_u = dl.inner(epsilon(v),sigma_damage(u,phi_old))dl.dx;
Equation_phi = ( Gcl0dl.inner(dl.grad(phi),dl.grad(theta))
+ (Gc/l0 + 2.0H(u_old,u_new,H_old) )dl.inner(phi,theta)
- 2.0H(u_old,u_new,H_old)*theta )*dl.dx;

define Problems

Problem_u = dl.LinearVariationalProblem(dl.lhs(Equation_u),dl.rhs(Equation_u),u_new,BC_u);
Problem_phi = dl.LinearVariationalProblem(dl.lhs(Equation_phi),dl.rhs(Equation_phi),phi_new,);

Solver_u = dl.LinearVariationalSolver(Problem_u);
Solver_phi = dl.LinearVariationalSolver(Problem_phi);

Solve

tol = 1e-4
Max_iter = 200;
iter = 1;
error = 1;

#Initial
u_old.vector()[:] = 0;
u_new.vector()[:] = 0;
phi_old.vector()[:] = 0;
phi_new.vector()[:] = 0;
H_old.vector()[:] = 0;

ud.ud = 0.0059 #displacement
while (error > tol) and (iter < Max_iter+1):
print(‘Iterations(Unconverged)=’,iter);
Solver_u.solve();
Solver_phi.solve();

error_u = dl.errornorm(u_new,u_old,norm_type = 'l2',mesh = None);
error_phi = dl.errornorm(phi_new,phi_old,norm_type = 'l2',mesh = None);
error = max(error_u,error_phi)

u_old.assign(u_new);
phi_old.assign(phi_new);
psi_new = dl.project(psi(u_new),V_psi)
H_old.assign(psi_new);


iter = iter + 1;

print("Convergence!");

Compute J-Integral

dim = 2;

#normal Vector
N = np.zeros([5,2]);
N[0,0] = -1; N[0,1] = 0; #Line1
N[1,0] = 0; N[1,1] = -1; #Line2
N[2,0] = 1; N[2,1] = 0; #Line3
N[3,0] = 0; N[3,1] = 1; #Line4
N[4,0] = -1; N[4,1] = 0; #Line5

def J_Caculation(mesh,u,phi,Gc,l0,ds):
V1 = dl.FunctionSpace(mesh,‘CG’,1);
V2 = dl.FunctionSpace(mesh,‘CG’,2);
sigma_ = sigma_damage(u,phi);
eps = epsilon(u);

def Term1(u,phi):
    w_e = dl.project((1/2)*dl.inner(sigma_,eps),V2);
    gamma = (1/(2*l0))*( (phi)**2 + l0**2 *  dl.dot(dl.grad(phi),dl.grad(phi)));
    gamma = dl.project(gamma,V2);                 
    w_c = Gc*gamma;
    w = w_e + w_c;                      
    return w; 

def Term2(u,phi,n):
    term2 = 0;
    for i in range(0,dim):
        ui = u[I];
        ui_dx = dl.grad(ui)[0];
        Ti = 0;   
        for j in range(0,dim):
            Ti += sigma_[i,j]*n[j];
        
        term2 += Ti*ui_dx;
    term2 = dl.project(term2,V2);
    return term2

def Term3(u,phi,n):
    term3 = Gc*l0*phi.dx(0)*(phi.dx(0)*n[0]+phi.dx(1)*n[1]);
    term3 = dl.project(term3,V2)
    return term3

Num = 5; 
J = 0;
J1 = 0;J2 = 0;J3 = 0;
for i in range(0,Num):
    n = N[i,:];
    Term = Term1(u,phi)*n[0]-Term2(u,phi,n)-Term3(u,phi,n);
    J1 += dl.assemble( Term1(u,phi)*n[0]*ds(11+i));
    J2 += dl.assemble( -Term2(u,phi,n)*ds(11+i));
    J3 += dl.assemble(-Term3(u,phi,n)*ds(11+i));

Blockquote

J = J1+J2+J3

return J

J = J_Caculation(mesh,u_new,phi_new,Gc,l0,ds_J);
print(J);

4

To make it more likely that anyone can help you with your code, please use 3x` encapsulation, i.e.

```python
from dolfin import *

# Add code here with correct indentation
def f(x):
    return x
```

I would also strongly suggest that you state the reference to the model by either writing out the mathematics, or adding the relevant reference.

5

Thank you for your suggestion. I will spend some time optimizing my code. In a little while, I will provide a more detailed explanation.

1 Like