Calculate J-integral for edge cracked domain


253
views
1
7 months ago by
Hi, I want to calculate J-integral for the edge crack domain for the domain as:

where J is given as:


and


How can I solve this:

Thanks in advance
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
set_log_active(False)
#==================================
# mesh (concident nodes) for the crack
#===================================
nx = 1
ny = 1
num = 1
mesh = Mesh()
editor = MeshEditor()
editor.open(mesh, 2, 2)
editor.init_vertices(10)
editor.init_cells(8)
editor.add_vertex(0, 0, 0.0)
editor.add_vertex(1, 0.5, 0.0)
editor.add_vertex(2, 0.5, 0.5)
editor.add_vertex(3, 0.0, 0.5)
editor.add_vertex(4, -0.5, 0.5)
editor.add_vertex(5, -0.5, 0.0)
editor.add_vertex(6, -0.5, -0.5)
editor.add_vertex(7, 0.0, -0.5)
editor.add_vertex(8, 0.5, -0.5)
editor.add_vertex(9, -0.5, 0.0)
editor.add_cell(0, 0, 1, 3)
editor.add_cell(1, 1, 2, 3)
editor.add_cell(2, 0, 3, 4)
editor.add_cell(3, 0, 4, 5)
editor.add_cell(4, 0, 9, 7)
editor.add_cell(5, 6, 7, 9)
editor.add_cell(6, 0, 7, 8)
editor.add_cell(7, 0, 8, 1)
editor.close()
#--------------------------------------
# change num_refine for mesh refinement
#--------------------------------------
num_refine = 5
h=0.5**num_refine
for i in range(num_refine):
    mesh=refine(mesh)
print "number of unknown",mesh.num_vertices()
print "number of elements",mesh.num_cells()
#plot(mesh)
#interactive()
# Create function space
V = VectorFunctionSpace(mesh, "CG", 2)
# Create test and trial functions, and source term
u, w = TrialFunction(V), TestFunction(V)
E, nu = 10.0, 0.3
mu, lmbda = E/(2.0*(1.0 + nu)), E*nu/((1.0 + nu)*(1.0 -2.0*nu))

sigma = 2*mu*sym(grad(u)) +lmbda*tr(grad(u))*Identity(2)

F = inner(sigma, grad(w))*dx 

a, L = lhs(F), rhs(F)

#---------------------------
# Boundary conditions
#---------------------------
class top(SubDomain):
    def inside(self,x,on_boundary):
        tol = 1e-10
        return abs(x[1]-0.5) < tol and on_boundary

class bottom(SubDomain):
    def inside(self,x,on_boundary):
        tol = 1e-10
        return abs(x[1]+0.5) < tol and on_boundary

Top = top()
Bottom = bottom()

bc_b = DirichletBC(V, Constant((0.0,0.0)), Bottom)
bc_t = DirichletBC(V.sub(1), Constant(0.01), Top)
bc = [bc_b, bc_t]
u = Function(V)
solve(a==L,u,bc)
u1,u2 = split(u)
plot(u1)
plot(u2)
interactive()

# Now I want to solve J-integral
​
Community: FEniCS Project
Just a small offtopic question, but why did you manually constructed the mesh instead of using the UnitSquareMesh object?
written 7 months ago by Murilo Moreira  
1
I need a crack and manual construction allowed to have conincident node which represents a crack.
written 7 months ago by hirshikesh  
Thank you for the explanation!
written 7 months ago by Murilo Moreira  
Hi,

if you figure this one out on your own, I would appreciate if you shared your solution :)

thanks!
written 7 months ago by klunkean  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »