User login


You are here

XFEM: crack tip function derivatives


I am using a X-FEM Matlab code found here"]http://people.civ...

In branch.m/branch_node.m (both functions of radius r and angle theta) crack tip fields are defined...and questions arise.

% Functions

f(1) = r2 * st2 ;
f(2) = r2 * ct2;
f(3) = r2 * st2 * ct;
f(4) = r2 * ct2 * ct;

-these functions are not the same defined e.g. in (Fleming,1997) or (Sukumar,2000) etc...why?

 their derivatives are then considered:

(here for function f(3))

% third function
dPhidr  = fac * st2 * (2*st*st-ct) ;
dPhidt  = fac * ct* cos(3*theta/2.);
dfdx(3) = dPhidr * drdx + dPhidt * dtdx  ;
dfdy(3) = dPhidr * drdy + dPhidt * dtdy  ;


-these not seem correct...(e.g. dPhidt should simply give fac * st2 * ct). But if I change these values, crack intensity factors are not correctly computed in the Matlab code. Do you have any idea? (I miss something for sure, but...)

In the following I give some notations: 


r2 = sqrt(r);

fac  = 0.5/r2 ;

st2  = sin(theta/2.);
ct2  = cos(theta/2.);
st   = sin(theta);
ct   = cos(theta);

drdx = cos(alpha);
drdy = sin(alpha);
dtdx = -sin(alpha);
dtdy = cos(alpha);


thanks for helping! 

phunguyen's picture

Hi Stephano,

 Sorry for very late reply however I did not notice your post.

I agree that the third and fourth functions are different from the one of another authors. However I did not find big differences between them.

 Concerning the derivatives, I used bad names for the variables. Using the chain rule, we write, f(r,theta)

dfdx =  dfdx1 * dx1dx + dfdx2 * dx2dx

dfdy =  dfdx1 * dx1dy + dfdx2 * dx2dy

 What you see in my code, dfdr is in fact dfdx1 and calculated in advance by hand as 

dfdx1 = dfdr * drdx1 + dfdt * dtdx1 

I just posted the corrected code with renamed variables.

Hope that it was a bit of help.



damoonmoetamedi's picture

Hi nguyen,

I am working on your matlab source.

 why you use crack angle (alpha) for calculating jacobian of drdx in branch.m.

Later, you use gauss point angle (tetha) to calculate drdx in j-integ?

Another problem,is not related to your matlab source, do you know how could I calculate dynamic J-Integ, because there are terms which are related to velocity of auxiliary displacement and strain.

thanks for your time,


Subscribe to Comments for "XFEM: crack tip function derivatives"

Recent comments

More comments


Subscribe to Syndicate