Dear professionals,
I'm self-developing explicit 3-D hybrid FEM/DEM code using GPGPU for my research regarding Penetration of projectile into rock.
In the code, I'm incorporating JH-2(Johnson holmquist) constitutive model.
Although this constitutive mode has been implemented in LS-DYNA and Abacus Explicit, some implementation detail has not been geven in exsisting literature. I mean, I reviewed several literature but I only found a flowchart of implementation of JH2 from Table 1 in
http://handle.dtic.mil/100.2/ADA402347
So, based on this document, I did my verification tests described in
http://www.dynalook.com/european-conf-2003/implementation-and-validation-of-the-johnson.pdf
and successfully passed the Case A with little error. But I failed in the Case B in which evaluation of effective plastic strain increment and damage evoluation is the key. My implementation of the model gave too rapid increase of damage. This could be due to a difference between my and Johnson and Holmquist’s definition of the plastic strain increment used to drive the damage. Actually same problem has been reported in https://www.diva-portal.org/smash/get/diva2:613844/FULLTEXT01.pdf where the author used Autodyn.
Anyway l regarded all the inelastic deformations as plastic , which seems to be not working. I think I have to consider strain linked to the (change of yield stress due to) damage increase. But the question is how??? Any knd help would be sincerely appreciated.
For your reference, the fortran code I made is shown here. I identified that the problem is coming from the hilightened part by red.
Thank you in advance.
Best regards,
Daisuke Fukuda
!@J H2 treatment(suppose trial cauthy stress Tij is T(1:3,1:3))
A_JH2 = 0.93d0 !@ Intact strength constant
B_JH2 = 0.0d0 !@ Fracture strength constant
C_JH2 = 0.0d0 !@ Strain rate constant
M_JH2 = 0.0d0 !@ Fracture strength constant
N_JH2 = 0.6d0 !@ Intact strength constant
D1_JH2 = 0.005d0 !@ Damage coefficient
D2_JH2 = 1.0d0 !@ Damage exponent
HEL_JH2 = 2.79d9 !@ Hugoniot elastic limit (HEL)
P_HEL_JH2 = 1.46d9 !@ Pressure corresponding to HEL
K1_JH2 = 130.95d9 !@ bulk modulus K, <-- In JH2 model, this is termed as K1
K2_JH2 = 0.0d0 !@ EOS(Equation of state) constant
K3_JH2 = 0.0d0 !@ EOS(Equation of state) constant
beta_JH2 = 1.0d0 !@ Bulking factor
threshhold_eps_rate = 1.0d0
tensileSt = 200.0d6
lamda_lame= K1_JH2 - 2.0d0/3.0d0* mu_lame !lame's constant:λ=(K-(2/3)μ)
mu_lame = 90.16d9 !@ lame's constant:μ, a.k.a. shear modulus G
Sig_HEL_JH2 = three/two*(HEL_JH2-P_HEL_JH2) !σHEL
!@ mean stress
mean_stress=(T(1,1)+T(2,2)+T(3,3))/three
!@ deviatric stress
S(1,1)=T(1,1)-mean_stress
S(2,2)=T(2,2)-mean_stress
S(3,3)=T(3,3)-mean_stress
!@ Equivalent stress σeq
temp_val = S(1,1)*S(1,1)+S(2,2)*S(2,2)+S(3,3)*S(3,3)&
+ two*(S(1,2)*S(1,2)+S(2,3)*S(2,3)+S(3,1)*S(3,1))
sig_eq = sqrt(three/two*temp_val)
!@ Effective strain rate, effective_eps_rate , D_tensor is rate of deformation tensor
temp_val = D_tensor(1,1)*D_tensor(1,1)+D_tensor(2,2)*D_tensor(2,2)+D_tensor(3,3)*D_tensor(3,3)&
+ two*(D_tensor(1,2)*D_tensor(1,2)+D_tensor(2,3)*D_tensor(2,3)+D_tensor(3,1)*D_tensor(3,1))
effective_eps_rate = sqrt(two/three*temp_val)
!@ normalized strain ε*
normalized_eps_rate = effective_eps_rate / threshhold_eps_rate
!@ Damage value of tetra element whose element number is "itet"
Damage = Damage_tet(itet)
!@ Pressure value of tetra element whose element number is "itet"
Press = Press_tet(itet)
!@ Normalized pressure P*
P_ast = Press / P_HEL_JH2
!@ Normalized tensile strenght T*
T_ast = tensileSt / P_HEL_JH2
!@ Normalized Intact strength σi*
if(normalized_eps_rate > one)then
temp_val = one + C_JH2*log(normalized_eps_rate)
else
temp_val = one
end if
Sig_i_ast = A_JH2 * (P_ast + T_ast)**N_JH2 * temp_val
!@ Normalized Fractured strength σf*
Sig_f_ast = B_JH2 * (P_ast**M_JH2) * temp_val
!@ Normalized strength σ*
Sig_ast = Sig_i_ast - Damage * ( Sig_i_ast - Sig_f_ast )
!@ Yield function, f_yield
f_yield = sig_eq - Sig_HEL_JH2 * Sig_ast
if(f_yield > zero)then !When yielding occurs
!@ return deviatoric stresses to the Yield surface using the Radial Return Algorithm,
S = S * (Sig_HEL_JH2 * Sig_ast / sig_eq)
!@ Equivalent stress σeq(new)
temp_val = S(1,1)*S(1,1)+S(2,2)*S(2,2)+S(3,3)*S(3,3)&
+ two*(S(1,2)*S(1,2)+S(2,3)*S(2,3)+S(3,1)*S(3,1))
sig_eq = sqrt(three/two*temp_val)
if(D1_JH2>zero.and.damage<one)then !when Damage gradually changes from 0 to 1
!@ compute the equivalent plastic strain increment：Δεp
del_eps_p =( S(1,1)*D_tensor(1,1)+ S(2,2)*D_tensor(2,2)+ S(3,3)*D_tensor(3,3)&
+two*( S(1,2)*D_tensor(1,2)+ S(2,3)*D_tensor(2,3)+ S(3,1)*D_tensor(3,1)) )&
/sig_eq*dt
if(del_eps_p<zero)del_eps_p=zero
!@ compute theplastic strain to failure (εp)_f
eps_p_f = D1_JH2 * ( P_ast +T_ast )**D2_JH2
!@ update the damage variable
Damage = Damage + del_eps_p/eps_p_f
if(Damage>one)Damage =one
else
Damage =one
end if
end if
!@ Compressibility factor μ_comp
mu_comp = (voli-volc)/volc !@ ,
!@ Regardless of sign of μ_comp, P includes the term "K1×μ_comp"
Press = K1_JH2 * mu_comp
if(mu_comp > zero)then
!@ compression-->mu_comp>0
!@ P=K1×μ_comp + K2×(μ_comp*μ_comp)+ K3×(μ_comp*μ_comp*μ_comp)
Press = Press + K2_JH2 * mu_comp* mu_comp + K3_JH2 * mu_comp * mu_comp * mu_comp
end if
del_P = zero
if(Damage > zero.and.mu_comp>zero)then
!@ σ* at previous time step is saved to temp_val
temp_val = Sig_ast_tet(itet)
!@ Energy loss due to damage:ΔU =(σHEL)^2/(6G)*{(previous σ*)^2-(current σ*)^2}
del_U = Sig_HEL_JH2 * Sig_HEL_JH2 /six /mu * ( temp_val*temp_val- Sig_ast*Sig_ast )
if(del_U<zero)del_U=zero
!@ additional pressure increment , note: del_P_tet(itet) is pressure increment at previous step
temp_val = K1_JH2 * mu_comp + del_P_tet(itet)
temp_val = temp_val*temp_val + two * beta_JH2 * K1_JH2 * del_U
del_P = -K1_JH2 * mu_comp + sqrt(temp_val)
end if
Press = Press + del_P
!@ save the following variables for next timestep
Damage_tet(itet) = Damage
Sig_ast_tet(itet) = Sig_ast
Press_tet(itet) = Press
del_P_tet(itet) = del_P
!@ σij(cauthy stress) = Sij(deviatoric stress) - Pδij
T = S
T(1,1) = T(1,1) +(- Press)
T(2,2) = T(2,2) +(- Press)
T(3,3) = T(3,3) +(- Press)
!@ end of stress modification by JH2
How can I implement JH-2(Johnson holmquist) constitutive equation in explicit FEM? - ResearchGate. Available from: https://www.researchgate.net/post/How_can_I_implement_JH-2Johnson_holmquist_constitutive_equation_in_explicit_FEM [accessed Oct 11, 2015].