module ElastoPlasticityClass use ArrayClass use IOClass use MathClass use LinearSolverClass use FEMDomainClass use FEMSolverClass use iso_fortran_env implicit none abstract interface function ScalarFunction(x,params) result(ret) use iso_fortran_env complex(real64), intent(in) :: x ! variables should be complex numbers real(real64), intent(in) ::params(:) complex(real64) :: ret! variables should be complex numbers end function end interface abstract interface function E_PotentialFunction(ElasticStrain, params) result(ret) use iso_fortran_env complex(real64), intent(in) :: ElasticStrain(:, :)! variables should be complex numbers real(real64), intent(in) :: params(:) complex(real64) :: ret! variables should be complex numbers end function end interface abstract interface function P_PotentialFunction(CauchyStress, PlasticStrain, params) result(ret) use iso_fortran_env complex(real64), intent(in) :: CauchyStress(:, :), PlasticStrain(:, :) ! variables should be complex numbers real(real64), intent(in) :: params(:) complex(real64) :: ret! variables should be complex numbers end function end interface abstract interface function StressRatioFunction(CauchyStress, dCauchyStress, StrainRatio) result(ret) use iso_fortran_env real(real64), intent(in) :: CauchyStress(:, :), dCauchyStress(:, :), StrainRatio(:, :) real(real64),allocatable :: ret(:,:) end function end interface type :: EP_Model_ procedure(E_PotentialFunction), nopass, pointer :: ElasticPotential => null() procedure(P_PotentialFunction), nopass, pointer :: YieldFunction => null() procedure(P_PotentialFunction), nopass, pointer :: PlasticPotential => null() ! optional >> procedure(StressRatioFunction), nopass, pointer :: StressRatio => null() contains procedure,public :: StiffnessMatrix => StiffnessMatrix_EP_model !procedure,public :: SmallStrainTensor => SmallStrainTensor_EP_model end type interface to_EP_Model module procedure to_EP_Model_ElastoPlastClass end interface type :: EP_Domain_ procedure(E_PotentialFunction), nopass, pointer :: ElasticPotential => null() procedure(P_PotentialFunction), nopass, pointer :: YieldFunction => null() procedure(P_PotentialFunction), nopass, pointer :: PlasticPotential => null() type(FEMDomain_) :: femdomain real(real64), allocatable :: YieldFunction_params(:, :) real(real64), allocatable :: PlasticPotential_params(:, :) real(real64), allocatable :: ElasticPotential_params(:, :) ! basic fields real(real64), allocatable :: CauchyStress_field(:, :, :) ! (stressid, gpid, elemid) real(real64), allocatable :: Strain_field(:, :, :) ! (stressid, gpid, elemid) real(real64), allocatable :: PlasticStrain_field(:, :, :) ! (stressid, gpid, elemid) ! incremental form (working mem.) real(real64), allocatable :: dCauchyStress_field(:, :, :) ! (stressid, gpid, elemid) real(real64), allocatable :: dStrain_field(:, :, :) ! (stressid, gpid, elemid) real(real64), allocatable :: PlasticStrain_field_n(:, :, :) ! (stressid, gpid, elemid) real(real64), allocatable :: displacement(:) contains procedure, public :: importField => importFieldEpDomain procedure, public :: exportField => exportFieldEpDomain end type type :: ElastoPlasticity_ type(FEMSolver_) :: femsolver type(EP_Domain_), allocatable :: ep_domain(:) real(real64) :: tol = dble(1.0e-5) real(real64) :: gravity_accel(1:3) = [0.0d0, 0.0d0, -9.810d0] integer(int32) :: MAX_NEWTON_LOOP_ITR = 10000 contains ! >>>> solver >>>> procedure, public :: init => initElastoPlasticity ! >> static >> procedure, public :: solve => solveElastoPlasticity procedure, public :: solve_increment => solve_increment_ElastoPlasticity ! << static << ! <<<< solver <<<< procedure, pass :: edit_YF_PP_ElastoPlasticity generic :: edit => edit_YF_PP_ElastoPlasticity procedure, public :: export => exportElastoPlasticity procedure, public :: exportField => exportFieldElastoPlasticity procedure, public :: importField => importFieldElastoPlasticity procedure, public :: get_internal_force => get_internal_forceElastoPlasticity procedure, public :: update_stress_for_increment & => update_stress_for_incrementElastoPlasticity procedure, public :: get_delta_internal_force => get_delta_internal_forceElastoPlasticity procedure, public :: get_external_force => get_external_forceElastoPlasticity procedure, public :: fill_zero_at_DBC => fill_zero_at_DBCElastoPlasticity procedure, public :: getYieldFunctionTemplate => getYieldFunctionTemplateElastoPlasticity procedure, public :: getPlasticStrain => getPlasticStrain_ElastoPlasticity procedure, public :: getPlasticStrain_n => getPlasticStrain_n_ElastoPlasticity procedure, public :: getStrain => getStrain_ElastoPlasticity procedure, public :: getCauchyStress => getCauchyStress_ElastoPlasticity procedure, public :: getdCauchyStress => getdCauchyStress_ElastoPlasticity procedure, public :: getTractionForce => getTractionForce_ElastoPlasticity procedure, public :: setPlasticStrain => setPlasticStrain_ElastoPlasticity procedure, public :: setCauchyStress => setCauchyStress_ElastoPlasticity procedure, public :: setdCauchyStress => setdCauchyStress_ElastoPlasticity procedure, public :: addStrain => addStrain_ElastoPlasticity procedure, public :: setdStrain => setdStrain_ElastoPlasticity procedure, public :: reset => rezeroElastoPlasticity procedure, public :: rezero => rezeroElastoPlasticity procedure, public :: I1 => I1ElastoPlasticity procedure, public :: J2 => J2ElastoPlasticity procedure, public :: I1_e => I1_e_ElastoPlasticity procedure, public :: J2_e => J2_e_ElastoPlasticity end type interface d_dsigma module procedure d_dsigma_P_PotentialFunction end interface interface d_depsilon module procedure d_depsilon_E_PotentialFunction end interface interface d2_depsilon2 module procedure d2_depsilon2_E_PotentialFunction end interface ! ! ! interface d_depsilon_p ! module procedure d_depsion_p_P_PotentialFunc ! end interface interface to_I1 module procedure to_I1_real64, to_I1_complex64 end interface interface to_J1 module procedure to_J1_real64, to_J1_complex64 end interface interface to_J2 module procedure to_J2_real64, to_J2_complex64 end interface interface to_J3 module procedure to_J3_real64, to_J3_complex64 end interface interface to_LodeAngle module procedure to_LodeAngle_real64, to_LodeAngle_complex64 end interface contains ! ############################################# function to_StressTensor(YieldFunction, PlasticPotential, Strain, dStrain, CauchyStress, PlasticStrain, & YieldParams, PlasticParams, ElasticParams, pval, epsilon, Jmat) & result(tr_CauchyStress) procedure(P_PotentialFunction) :: YieldFunction procedure(P_PotentialFunction) :: PlasticPotential real(real64), intent(in) :: Strain(:, :), dStrain(:, :), ElasticParams(:), & YieldParams(:), PlasticParams(:), CauchyStress(:, :) real(real64), intent(inout) :: PlasticStrain(:, :) real(real64), optional, intent(inout) :: pval real(real64), optional, allocatable, intent(inout) :: Jmat(:, :) real(real64), intent(in) :: epsilon real(real64), allocatable :: tr_CauchyStress(:, :), dfds(:, :), E(:, :, :, :), dCauchyStress(:, :), & old_PlasticStrain(:, :), J_mat(:, :), X_vec(:), Y_vec(:), E_dfds(:, :), dX_vec(:), E_ee(:, :), & E11_epsilon(:, :), dfds_inv(:, :), dfds_forward(:, :), dfds_backward(:, :), ds(:, :), dE_e(:, :), & tr_sigma_vec(:), D_mat(:, :), tr_0_sigma(:, :), tr_f_CauchyStress(:, :), tr_b_CauchyStress(:, :), & E_ijkl(:, :, :, :), En(:, :), Ee(:, :), dgds(:, :) real(real64) :: f_val, dGamma, dGamma_lower, dGamma_upper, ddgamma, buf, G_val, K_val, eps, f_n_0, dg1, dg2 integer(int32) :: i, j, k, l, n, m character(:), allocatable :: algorithm real(real64) :: f_n, f_n1, dfdgamma, f_forward, f_backward integer(int32) :: max_itr type(IO_) :: f algorithm = "ReturnMapping" tr_CauchyStress = zeros(3, 3) ! (1) dCauchyStress = StVenant_ConstModel(ElasticStrain=dStrain, params=ElasticParams) ! (2) tr_CauchyStress = CauchyStress + dCauchyStress ! (3) f_val = real(YieldFunction(CauchyStress=dcmplx(tr_CauchyStress), & PlasticStrain=dcmplx(PlasticStrain), params=YieldParams)) dgamma = 0.0d0 f_n_0 = f_val f_n = f_val if (present(pval)) then pval = f_val end if dfds = zeros(3, 3) if (is_elastic(f_val)) then return else ! [Caution!] ! only for plastic potential with no hardening/softerning parameters associated with the plastic strain ! For such cases, Return mapping is required. if (algorithm == "ReturnMapping") then ! only for f(\sigma), not for f(\sigma, K, Hm ...etc.) max_itr = 4 !tr_0_sigma = tr_CauchyStress dgamma = 0.0d0 !dfds = d_dSigma(PlasticPotential=PlasticPotential,CauchyStress=tr_CauchyStress,& ! PlasticStrain=PlasticStrain,params=PlasticParams,epsilon=epsilon) !dfds = dfds/norm(dfds) dfdgamma = 0.0d0 f_n = f_val !print *, zfill(0,4),dgamma,dfdgamma,f_val dgds = d_dSigma(PlasticPotential=PlasticPotential, CauchyStress=tr_CauchyStress, & PlasticStrain=PlasticStrain, params=PlasticParams, epsilon=epsilon) dgds = dgds/norm(dgds) En = StVenant_ConstModel(ElasticStrain=dgds, params=ElasticParams) Ee = StVenant_ConstModel(ElasticStrain=dStrain, params=ElasticParams) dfds = d_dSigma(PlasticPotential=YieldFunction, CauchyStress=tr_CauchyStress, & PlasticStrain=PlasticStrain, params=YieldParams, epsilon=epsilon) dgamma = tensordot(dfds, Ee)/tensordot(dfds, En) dCauchyStress = StVenant_ConstModel(ElasticStrain=dStrain - dgamma*dgds, params=ElasticParams) tr_CauchyStress = CauchyStress + dCauchyStress f_n = real(YieldFunction(CauchyStress=dcmplx(tr_CauchyStress), & PlasticStrain=dcmplx(PlasticStrain), params=YieldParams)) PlasticStrain = PlasticStrain + dgamma*dgds return elseif (algorithm == "ForwardEuler") then ! Forward Euler dfds = d_dSigma(PlasticPotential=PlasticPotential, CauchyStress=tr_CauchyStress, & PlasticStrain=PlasticStrain, params=PlasticParams, epsilon=epsilon) E = StVenant_StiffnessMatrix(params=ElasticParams) dGamma_upper = 0.0d0 dGamma_lower = 0.0d0 do i = 1, 3 do j = 1, 3 do k = 1, 3 do l = 1, 3 dGamma_upper = dGamma_upper + dfds(i, j)*E(i, j, k, l)*dStrain(k, l) dGamma_lower = dGamma_lower + dfds(i, j)*E(i, j, k, l)*dfds(k, l) end do end do end do end do dGamma = dGamma_upper/dGamma_lower PlasticStrain = old_PlasticStrain + dGamma*dfds ! Forward Euler ! ds_{ij} = E_{ijkl}( d\epsilon_{kl} - d \gamma*df/ds_{kl} ) dCauchyStress = StVenant_ConstModel(ElasticStrain=dStrain - dgamma*dfds, params=ElasticParams) tr_CauchyStress = CauchyStress + dCauchyStress ! check yield function f_val = real(PlasticPotential(CauchyStress=dcmplx(tr_CauchyStress), & PlasticStrain=dcmplx(PlasticStrain), params=PlasticParams)) if (present(pval)) then pval = f_val end if end if end if end function to_StressTensor ! ############################################# ! ############################################# function to_dStressTensor(YieldFunction, PlasticPotential, dStrain, CauchyStress, PlasticStrain, & YieldParams, PlasticParams, ElasticParams, pval, epsilon, new_PlasticStrain) & result(dCauchyStress) procedure(P_PotentialFunction) :: YieldFunction procedure(P_PotentialFunction) :: PlasticPotential !! return increment of Cauchy tensor real(real64), intent(in) :: dStrain(:, :), ElasticParams(:), & YieldParams(:), PlasticParams(:), CauchyStress(:, :) real(real64), intent(in) :: PlasticStrain(:, :) real(real64), intent(in) :: epsilon real(real64), optional, intent(inout) :: pval real(real64), optional, allocatable, intent(inout) :: new_PlasticStrain(:, :) ! local variables real(real64), allocatable :: tr_CauchyStress(:, :), dfds(:, :), E(:, :, :, :), dCauchyStress(:, :), & old_PlasticStrain(:, :), J_mat(:, :), X_vec(:), Y_vec(:), E_dfds(:, :), dX_vec(:), E_ee(:, :), & E11_epsilon(:, :), dfds_inv(:, :), dfds_forward(:, :), dfds_backward(:, :), ds(:, :), dE_e(:, :), & tr_sigma_vec(:), D_mat(:, :), tr_0_sigma(:, :), tr_f_CauchyStress(:, :), tr_b_CauchyStress(:, :), & E_ijkl(:, :, :, :), En(:, :), Ee(:, :), dgds(:, :) real(real64) :: f_val, dGamma, dGamma_lower, dGamma_upper, ddgamma, buf, G_val, K_val, eps, f_n_0, dg1, dg2 integer(int32) :: i, j, k, l, n, m character(:), allocatable :: algorithm real(real64) :: f_n, f_n1, dfdgamma, f_forward, f_backward integer(int32) :: max_itr type(IO_) :: f new_PlasticStrain = PlasticStrain algorithm = "ReturnMapping" tr_CauchyStress = zeros(3, 3) ! (1) dCauchyStress = StVenant_ConstModel(ElasticStrain=dStrain, params=ElasticParams) ! (2) tr_CauchyStress = CauchyStress + dCauchyStress ! (3) f_val = real(YieldFunction(CauchyStress=dcmplx(tr_CauchyStress), & PlasticStrain=dcmplx(new_PlasticStrain), params=YieldParams)) dgamma = 0.0d0 f_n_0 = f_val f_n = f_val if (present(pval)) then pval = f_val end if dfds = zeros(3, 3) if (is_elastic(f_val)) then return else ! [Caution!] ! only for plastic potential with no hardening/softerning parameters associated with the plastic strain ! For such cases, Return mapping is required. if (algorithm == "ReturnMapping") then ! only for f(\sigma), not for f(\sigma, K, Hm ...etc.) max_itr = 4 dgamma = 0.0d0 dfdgamma = 0.0d0 f_n = f_val dgds = d_dSigma(PlasticPotential=PlasticPotential, CauchyStress=tr_CauchyStress, & PlasticStrain=new_PlasticStrain, params=PlasticParams, epsilon=epsilon) dgds = dgds/norm(dgds) En = StVenant_ConstModel(ElasticStrain=dgds, params=ElasticParams) Ee = StVenant_ConstModel(ElasticStrain=dStrain, params=ElasticParams) dfds = d_dSigma(PlasticPotential=YieldFunction, CauchyStress=tr_CauchyStress, & PlasticStrain=new_PlasticStrain, params=YieldParams, epsilon=epsilon) dgamma = tensordot(dfds, Ee)/tensordot(dfds, En) dCauchyStress = StVenant_ConstModel(ElasticStrain=dStrain - dgamma*dgds, params=ElasticParams) tr_CauchyStress = CauchyStress + dCauchyStress f_n = real(YieldFunction(CauchyStress=dcmplx(tr_CauchyStress), & PlasticStrain=dcmplx(new_PlasticStrain), params=YieldParams)) new_PlasticStrain = new_PlasticStrain + dgamma*dgds !dCauchyStress = tr_CauchyStress - CauchyStress return else print *, "ERROR :: to_dStressTensor >> no such algorithm as", algorithm stop end if end if end function to_dStressTensor ! ############################################# ! ############################################# function neoHookean(ElasticStrain, params) result(ret) ! E_PotentialFunction ! Simo and Pister, 1984 complex(real64), intent(in) :: ElasticStrain(:, :) real(real64), intent(in) :: params(:) complex(real64) :: ret real(real64) :: lambda, mu complex(real64) :: det_Ce, J complex(real64),allocatable :: C_e(:,:) integeR(int32) :: n lambda = params(1) mu = params(2) n = size(ElasticStrain,1) C_e = 2.0d0*ElasticStrain + eyes(n,n) det_Ce = det_mat(C_e,size(ElasticStrain,1)) J = sqrt(det_Ce) ret = lambda/2.0d0*(log(J))**2 - mu*log(J) + mu/2.0d0*(trace(C_e) - 3.0d0) end function ! ############################################# function neoHookean_Vladimirov(ElasticStrain, params) result(ret) ! E_PotentialFunction ! Vladimirov, I.N., Pietryga, M.P., Reese, S., 2010. Anisotropic finite elastoplasticity ! with nonlinear kinematic and isotropic hardening and application to sheet metal forming. ! Int. J. Plast. 26, 659–687. complex(real64), intent(in) :: ElasticStrain(:, :) real(real64), intent(in) :: params(:) complex(real64) :: ret real(real64) :: lambda, mu complex(real64) :: det_Ce complex(real64),allocatable :: C_e(:,:) integeR(int32) :: n lambda = params(1) mu = params(2) n = size(ElasticStrain,1) C_e = 2.0d0*ElasticStrain + eyes(n,n) det_Ce = det_mat(C_e,size(ElasticStrain,1)) ret = mu/2.0d0*(trace(C_e) - 3.0d0 ) - mu*log(sqrt(det_Ce)) & + lambda/4.0d0*( & det_Ce & - 1.0d0 & - 2.0d0*log(sqrt(det_Ce ) ) & ) end function ! ############################################# function neoHookean_Simo(ElasticStrain, params) result(ret) ! E_PotentialFunction ! Simo and Pister, 1984 complex(real64), intent(in) :: ElasticStrain(:, :) real(real64), intent(in) :: params(:) complex(real64) :: ret real(real64) :: lambda, mu complex(real64) :: det_Ce, J complex(real64),allocatable :: C_e(:,:) integeR(int32) :: n lambda = params(1) mu = params(2) n = size(ElasticStrain,1) C_e = 2.0d0*ElasticStrain + eyes(n,n) det_Ce = det_mat(C_e,size(ElasticStrain,1)) J = sqrt(det_Ce) ret = lambda/4.0d0*(det_Ce-1.0d0) - (lambda/2.0d0+mu)*log(J) + mu/2.0d0*(trace(C_e) - 3.0d0) end function ! ############################################# function StVenant(ElasticStrain, params) result(ret) ! E_PotentialFunction complex(real64), intent(in) :: ElasticStrain(:, :) real(real64), intent(in) :: params(:) complex(real64) :: ret real(real64) :: lambda, mu lambda = params(1) mu = params(2) ret = lambda/2.0d0*(trace(ElasticStrain))**2 + mu*trace(matmul(ElasticStrain,ElasticStrain)) end function ! ############################################# function StVenant_ConstModel(ElasticStrain, params) result(CauchyStress) real(real64), intent(in) :: ElasticStrain(:, :), params(:) real(real64), allocatable :: CauchyStress(:, :) real(real64) :: lambda, mu lambda = params(1) mu = params(2) CauchyStress = lambda*eyes(3, 3)*trace(ElasticStrain) & + 2.0d0*mu*ElasticStrain end function ! ############################################# ! ############################################# function StVenant_StiffnessMatrix(params) result(E) real(real64), intent(in) :: params(:) real(real64), allocatable :: E(:, :, :, :), g(:, :) real(real64) :: lambda, mu integer(int32) :: i, j, k, l allocate (E(3, 3, 3, 3)) g = eyes(3, 3) lambda = params(1) mu = params(2) !https://ss1.xrea.com/penguinitis.g1.xrea.com/study/note/elastic_coefficient.pdf do i = 1, 3 do j = 1, 3 do k = 1, 3 do l = 1, 3 E(i, j, k, l) = lambda*lambda*g(i, j)*g(k, l) & + mu*(g(i, k)*g(j, l) + g(i, l)*g(j, k)) end do end do end do end do end function ! ############################################# function StVenant_StiffnessMatrix_2D(params) result(E) real(real64), intent(in) :: params(:) real(real64), allocatable :: E(:, :) real(real64) :: lambda, mu integer(int32) :: i, j, k, l allocate (E(6, 6)) lambda = params(1) mu = params(2) E(1, 1) = 2.0d0*mu + lambda E(1, 2) = lambda E(1, 3) = lambda E(2, 1) = lambda E(2, 2) = 2.0d0*mu + lambda E(2, 3) = lambda E(3, 1) = lambda E(3, 2) = lambda E(3, 3) = 2.0d0*mu + lambda E(4, 4) = mu E(5, 5) = mu E(6, 6) = mu end function ! ############################################# function MohrCoulomb(CauchyStress, PlasticStrain, params) result(ret) ! ! P_PotentialFunction complex(real64), intent(in) :: CauchyStress(:, :), PlasticStrain(:, :) real(real64), intent(in) :: params(:) complex(real64) :: ret complex(real64) :: c, phi, I_1, J_2, J_3, theta ! https://static.rocscience.cloud/assets/verification-and-theory/RSData/mohr-coulomb-model.pdf c = params(1) phi = params(2) I_1 = to_I1(CauchyStress) J_2 = to_J2(CauchyStress) J_3 = to_J3(CauchyStress) theta = to_LodeAngle(CauchyStress) ret = I_1/ 3.0d0*sin(phi) & + sqrt(J_2)*(cos(theta) - 1.0d0/sqrt(3.0d0)*sin(theta)*sin(phi)) & - c*cos(phi) end function ! ################################################## function DruckerPrager(CauchyStress, PlasticStrain, params) result(ret) ! ! P_PotentialFunction complex(real64), intent(in) :: CauchyStress(:, :), PlasticStrain(:, :) real(real64), intent(in) :: params(:) complex(real64) :: ret complex(real64) :: c, phi, I_1, J_2, J_3, theta, A, B ! https://static.rocscience.cloud/assets/verification-and-theory/RSData/mohr-coulomb-model.pdf c = params(1) phi = params(2) I_1 = to_I1(CauchyStress) J_2 = to_J2(CauchyStress) J_3 = to_J3(CauchyStress) theta = to_LodeAngle(CauchyStress) !https://en.wikipedia.org/wiki/Drucker%E2%80%93Prager_yield_criterion ! middle circumscribes A = 6.0d0*c*cos(phi)/sqrt(3.0d0)/(3.0d0 + sin(phi)) B = 2.0d0*sin(phi)/sqrt(3.0d0)/(3.0d0 + sin(phi)) !ret = sqrt(J_2) - A - B*I_1 ret = sqrt(J_2) - abs(A) - abs(B)*I_1 end function ! ################################################## function VonMises(CauchyStress, PlasticStrain, params) result(ret) ! ! P_PotentialFunction complex(real64), intent(in) :: CauchyStress(:, :), PlasticStrain(:, :) real(real64), intent(in) :: params(:) complex(real64) :: ret complex(real64) :: k, J_2 !https://www.engineersedge.com/material_science/von_mises.htm k = params(1) J_2 = to_J2(CauchyStress) ret = sqrt(J_2) - k end function ! ################################################## function Tresca(CauchyStress, PlasticStrain, params) result(ret) ! ! P_PotentialFunction complex(real64), intent(in) :: CauchyStress(:, :), PlasticStrain(:, :) real(real64), intent(in) :: params(:) complex(real64) :: ret complex(real64) :: c, phi, J_2, theta,f_t !https://www.engineersedge.com/material_science/von_mises.htm c = params(1) phi = params(2) J_2 = to_J2(CauchyStress) theta = to_LodeAngle(CauchyStress) f_t = 2.0d0*c*cos(phi)/(1.0d0+sin(phi)) ret = sqrt(J_2) - abs(f_t/(2.0d0*cos(theta))) !ret = sqrt(J_2) - sqrt(abs(f_t/(2.0d0*cos(theta)))) !ret = 4.0d0*(J_2**3) & ! - 27.0d0*(J_3**2) & ! - 36.0d0*(k**2)*(J_2**2) & ! + 96.0d0*(k**4)*(J_2) & ! - 64.0d0*(k**6) end function ! ################################################## function CamClay(CauchyStress, PlasticStrain, params) result(ret) ! ! P_PotentialFunction complex(real64), intent(in) :: CauchyStress(:, :), PlasticStrain(:, :) real(real64), intent(in) :: params(:) complex(real64) :: ret complex(real64) :: c, phi, J_2, J_3, p, q, M, theta, pc !https://www.engineersedge.com/material_science/von_mises.htm !http://manual.midasuser.com/JP_Common/FEANX/110/FEA_NX/%E3%83%A1%E3%83%83%E3%82%B7%E3%83%A5/%E7%89%B9%E6%80%A7_%E5%BA%A7%E6%A8%99%E7%B3%BB_%E9%96%A2%E6%95%B0/%E6%9D%90%E6%96%99/%E6%9D%90%E6%96%99%E4%B8%80%E8%88%AC/Modified_Cam_Clay.htm stop "ERROR >> CamClay is not correctly implemented." c = params(1) phi = params(2) pc = params(3) J_2 = to_J2(CauchyStress) J_3 = to_J3(CauchyStress) theta = to_LodeAngle(CauchyStress) p = -3.0d0*to_I1(CauchyStress) q = 3.0d0*J_2 M = 6.0d0*sin(phi)/(3.0d0 - sin(phi)) ! 厳密には正しくない ret = 3.0d0*J_2 + M*p*log(p/pc) end function ! ################################################## function ModifiedCamClay(CauchyStress, PlasticStrain, params) result(ret) complex(real64), intent(in) :: CauchyStress(:, :), PlasticStrain(:, :) real(real64), intent(in) :: params(:) complex(real64) :: ret complex(real64) :: c, phi, J_2, J_3, p, q, M, theta, pc !https://www.engineersedge.com/material_science/von_mises.htm !http://manual.midasuser.com/JP_Common/FEANX/110/FEA_NX/%E3%83%A1%E3%83%83%E3%82%B7%E3%83%A5/%E7%89%B9%E6%80%A7_%E5%BA%A7%E6%A8%99%E7%B3%BB_%E9%96%A2%E6%95%B0/%E6%9D%90%E6%96%99/%E6%9D%90%E6%96%99%E4%B8%80%E8%88%AC/Modified_Cam_Clay.htm !http://docs.itascacg.com/3dec700/common/models/camclay/doc/modelcamclay.html c = params(1) phi = params(2) pc = params(3) J_2 = to_J2(CauchyStress) J_3 = to_J3(CauchyStress) theta = to_LodeAngle(CauchyStress) p = -3.0d0*to_I1(CauchyStress) !q = sqrt(3.0d0*J_2) M = 6.0d0*sin(phi)/(3.0d0 - sin(phi)) ! 厳密には正しくない ret = 3.0d0*J_2 + M*p*(p - pc) end function ! ################################################## ! ################################################## function to_I1_complex64(CauchyStress) result(ret) complex(real64), intent(in) :: CauchyStress(:, :) complex(real64) :: ret ret = trace(CauchyStress)/3.0d0 end function ! ################################################## ! ################################################## function to_I1_real64(CauchyStress) result(ret) real(real64), intent(in) :: CauchyStress(:, :) real(real64) :: ret ret = trace(CauchyStress)/3.0d0 end function ! ################################################## ! ################################################## function tensordot(a_ij, b_ij) result(ret) real(real64), intent(in) :: a_ij(:, :), b_ij(:, :) real(real64) :: ret integer(int32) :: i, j ret = 0.0d0 do i = 1, size(a_ij, 1) do j = 1, size(a_ij, 1) ret = ret + a_ij(i, j)*b_ij(i, j) end do end do end function ! ################################################## ! ################################################## function tensorSelfDot(a_ij) result(ret) real(real64), intent(in) :: a_ij(:, :) real(real64) :: ret integer(int32) :: i, j ret = 0.0d0 do i = 1, size(a_ij, 1) do j = 1, size(a_ij, 1) ret = ret + a_ij(i, j)*a_ij(i, j) end do end do end function ! ################################################## ! ################################################## function to_J1_complex64(CauchyStress) result(ret) complex(real64), intent(in) :: CauchyStress(:, :) complex(real64) :: ret ret = 0.0d0 end function ! ################################################## ! ################################################## function to_J1_real64(CauchyStress) result(ret) real(real64), intent(in) :: CauchyStress(:, :) real(real64) :: ret ret = 0.0d0 end function ! ################################################## ! ################################################## function to_J2_complex64(CauchyStress) result(ret) complex(real64), intent(in) :: CauchyStress(:, :) complex(real64) :: ret ! https://en.wikipedia.org/wiki/Cauchy_stress_tensor !ret = tensorSelfDot(to_DeviatricStress(CauchyStress))*0.50d0 ret = 0.50d0*( & trace(matmul(CauchyStress, CauchyStress)) & - trace(CauchyStress)*trace(CauchyStress)/3.0d0 & ) end function ! ################################################## ! ################################################## function to_J2_real64(CauchyStress) result(ret) real(real64), intent(in) :: CauchyStress(:, :) real(real64) :: ret ! https://en.wikipedia.org/wiki/Cauchy_stress_tensor !ret = tensorSelfDot(to_DeviatricStress(CauchyStress))*0.50d0 ret = 0.50d0*( & trace(matmul(CauchyStress, CauchyStress)) & - trace(CauchyStress)*trace(CauchyStress)/3.0d0 & ) end function ! ################################################## ! ################################################## function to_J3_complex64(CauchyStress) result(ret) complex(real64), intent(in) :: CauchyStress(:, :) complex(real64) :: ret ! https://en.wikipedia.org/wiki/Cauchy_stress_tensor ret = ( & trace(matmul(CauchyStress, matmul(CauchyStress, CauchyStress))) & - trace(matmul(CauchyStress, CauchyStress))*trace(CauchyStress) & + 2.0d0/9.0d0*trace(CauchyStress)*trace(CauchyStress)*trace(CauchyStress) & )/3.0d0 end function ! ################################################## ! ################################################## function to_J3_real64(CauchyStress) result(ret) real(real64), intent(in) :: CauchyStress(:, :) real(real64) :: ret ! https://en.wikipedia.org/wiki/Cauchy_stress_tensor ret = ( & trace(matmul(CauchyStress, matmul(CauchyStress, CauchyStress))) & - trace(matmul(CauchyStress, CauchyStress))*trace(CauchyStress) & + 2.0d0/9.0d0*trace(CauchyStress)*trace(CauchyStress)*trace(CauchyStress) & )/3.0d0 end function ! ################################################## ! ################################################## function to_LodeAngle_complex64(CauchyStress) result(ret) complex(real64), intent(in) :: CauchyStress(:, :) complex(real64) :: ret, J_2, J_3 J_2 = to_J2(CauchyStress) J_3 = to_J3(CauchyStress) if (J_2 == 0.0d0) then ret = 0.0d0 else if (abs(3.0d0*sqrt(3.0d0)/2.0d0*J_3/(J_2)**(1.50d0)) > 1.0d0) then if (abs(3.0d0*sqrt(3.0d0)/2.0d0*J_3/(J_2)**(1.50d0)) > 0.0d0) then ret = 1.0d0/3.0d0*asin(1.0d0) else ret = 1.0d0/3.0d0*asin(-1.0d0) end if else ret = 1.0d0/3.0d0*asin(-3.0d0*sqrt(3.0d0)/2.0d0*J_3/(J_2)**(1.50d0)) end if end if end function ! ################################################## ! ################################################## function to_LodeAngle_real64(CauchyStress) result(ret) real(real64), intent(in) :: CauchyStress(:, :) real(real64) :: ret, J_2, J_3 J_2 = to_J2(CauchyStress) J_3 = to_J3(CauchyStress) if (J_2 == 0.0d0) then ret = 0.0d0 else if (abs(3.0d0*sqrt(3.0d0)/2.0d0*J_3/(J_2)**(1.50d0)) > 1.0d0) then if (3.0d0*sqrt(3.0d0)/2.0d0*J_3/(J_2)**(1.50d0) > 0.0d0) then ret = 1.0d0/3.0d0*asin(1.0d0) else ret = 1.0d0/3.0d0*asin(-1.0d0) end if else ret = 1.0d0/3.0d0*asin(-3.0d0*sqrt(3.0d0)/2.0d0*J_3/(J_2)**(1.50d0)) end if end if end function ! ################################################## ! ################################################## function to_DeviatricStress(CauchyStress) result(ret) real(real64), intent(in) :: CauchyStress(:, :) real(real64), allocatable :: ret(:, :) ret = CauchyStress - trace(CauchyStress)/3.0d0*eyes(size(CauchyStress, 1), size(CauchyStress, 1)) end function ! ################################################## function d_dsigma_P_PotentialFunction(PlasticPotential, CauchyStress, PlasticStrain, params, epsilon) result(ret) procedure(P_PotentialFunction) :: PlasticPotential real(real64), intent(in) :: CauchyStress(:, :), PlasticStrain(:, :), params(:) real(real64), optional, intent(in) :: epsilon real(real64), allocatable :: dsigma_tensor(:, :), ret(:, :) real(real64) :: dsigma, df integer(int32) :: i, j type(Math_) :: math ret = zeros(size(CauchyStress, 1), size(CauchyStress, 2)) dsigma_tensor = zeros(size(CauchyStress, 1), size(CauchyStress, 2)) if (present(epsilon)) then dsigma = epsilon/2.0d0 else dsigma = dble(1.0e-16)/2.0d0 end if do i = 1, size(CauchyStress, 1) do j = i, size(CauchyStress, 2) dsigma_tensor(:, :) = 0.0d0 dsigma_tensor(i, j) = dsigma df = aimag(& PlasticPotential( dcmplx(CauchyStress) + dcmplx(dsigma_tensor)*math%i, dcmplx(PlasticStrain), params) & - PlasticPotential( dcmplx(CauchyStress) - dcmplx(dsigma_tensor)*math%i, dcmplx(PlasticStrain), params)& ) df = df/(dsigma*2.0d0) ret(i, j) = df ret(j, i) = df end do end do end function ! ################################################## ! ################################################## function d_depsilon_E_PotentialFunction(ElasticPotential, ElasticStrain, params, epsilon) result(ret) procedure(E_PotentialFunction) :: ElasticPotential real(real64), intent(in) :: ElasticStrain(:, :), params(:) real(real64), optional, intent(in) :: epsilon real(real64), allocatable :: dsigma_tensor(:, :), ret(:, :) real(real64) :: dsigma real(real64) :: df integer(int32) :: i, j type(Math_) :: math ret = zeros(size(ElasticStrain, 1), size(ElasticStrain, 2)) dsigma_tensor = zeros(size(ElasticStrain, 1), size(ElasticStrain, 2)) if (present(epsilon)) then dsigma = epsilon/2.0d0 else dsigma = dble(1.0e-16)/2.0d0 end if do i = 1, size(ElasticStrain, 1) do j = i, size(ElasticStrain, 2) dsigma_tensor(:, :) = 0.0d0 dsigma_tensor(i, j) = dsigma df = aimag(& ElasticPotential( dcmplx(ElasticStrain) + dcmplx(dsigma_tensor)*math%i, params) & - ElasticPotential( dcmplx(ElasticStrain) - dcmplx(dsigma_tensor)*math%i, params)& ) df = df/(dsigma*2.0d0) ret(i, j) = df ret(j, i) = df end do end do end function ! ################################################## ! ################################################## function d2_depsilon2_E_PotentialFunction(ElasticPotential, ElasticStrain, params, epsilon) result(ret) ! stiffness tensor procedure(E_PotentialFunction) :: ElasticPotential real(real64), intent(in) :: ElasticStrain(:, :), params(:) real(real64), optional, intent(in) :: epsilon real(real64), allocatable :: ret(:, :, :, :),dElasticStrain(:, :), dS(:,:) real(real64) :: depsilon integer(int32) :: i, j, k, l type(Math_) :: math complex(real64) :: df type(IO_) :: debug ret = zeros(size(ElasticStrain, 1), size(ElasticStrain, 2),& size(ElasticStrain, 1), size(ElasticStrain, 2)) dElasticStrain = zeros(size(ElasticStrain, 1), size(ElasticStrain, 2)) if (present(epsilon)) then depsilon = epsilon/2.0d0 else depsilon = dble(1.0e-8)/2.0d0 end if ! 中心差分 do k = 1, size(ElasticStrain, 1) do l = 1, size(ElasticStrain, 2) dElasticStrain(:,:) = 0.0d0 dElasticStrain(k, l) = depsilon dElasticStrain(l, k) = depsilon dS = d_depsilon(ElasticPotential,ElasticStrain+dElasticStrain, params, epsilon) & - d_depsilon(ElasticPotential,ElasticStrain-dElasticStrain, params, epsilon) do i=1,size(ElasticStrain, 1) do j=1,size(ElasticStrain, 2) ret(i,j,k,l) = ret(i,j,k,l) + dS(i,j)/(2.0d0*depsilon) enddo enddo enddo enddo ! ! (i,j) == (k,l) ! do i = 1, size(ElasticStrain, 1) ! do j = i, size(ElasticStrain, 2) ! dElasticStrain(:,:) = 0.0d0 ! dElasticStrain(i, j) = dsigma ! ! if(i/=j)then ! dElasticStrain(j, i) = dsigma ! endif ! ! ! df = real(ElasticPotential( dcmplx(ElasticStrain), params)) & ! - real(ElasticPotential( dcmplx(ElasticStrain+math%i*dElasticStrain), params) ) ! ! df = 2.0d0*df/(dsigma)/(dsigma) ! if(i/=j)then ! df = df/4.0d0 ! endif ! ! ret(i, j, i, j) = real(df) ! ret(i, j, j, i) = real(df) ! ret(j, i, i, j) = real(df) ! ret(j, i, j, i) = real(df) ! ! end do ! end do ! return ! ! (i,j) /= (k,l) ! do i = 1, size(ElasticStrain, 1) ! do j = i, size(ElasticStrain, 2) ! do k = 1, size(ElasticStrain, 1) ! do l = k, size(ElasticStrain, 2) ! if(i==k .and. j==l) cycle ! if(i==l .and. j==k) cycle ! ! dsigma_tensor(:, :) = 0.0d0 ! dsigma_tensor(i, j) = dsigma ! dsigma_tensor(j, i) = dsigma ! ! dElasticStrain(:, :) = 0.0d0 ! dElasticStrain(k, l) = dsigma ! dElasticStrain(l, k) = dsigma ! ! df = 1.0d0/dsigma/dsigma*ElasticPotential( dcmplx(ElasticStrain), params) & ! - 1.0d0/dsigma/dsigma*& ! real(ElasticPotential( & ! dcmplx(ElasticStrain+math%i*dsigma_tensor+math%i*dElasticStrain), params)) & ! - 0.50d0*ret(i, j, i, j) - 0.50d0*ret(k, l, k, l) ! ! ! ! df = V = dW/dE|_(E) ! ! ddf = dV/dE = im(dW/dE|_(E+ih) - dW/dE|_(E+ih))/2h ! ret(i, j, k, l) = abs(real(df)) ! !ret(j, i, k, l) = abs(real(df)) ! !ret(k, l, i, j) = abs(real(df)) ! !ret(l, k, i, j) = abs(real(df)) ! end do ! end do ! end do ! end do ! end function ! ################################################## ! ################################################## subroutine getYieldFunctionTemplateElastoPlasticity(this, name) class(ElastoPlasticity_), intent(in) :: this character(*), intent(in) :: name type(IO_) :: f if (index(name, ".f90") == 0) then call f%open(name + ".f90", "w") else call f%open(name, "w") end if call f%write("function MyYieldFunction(CauchyStress,PlasticStrain,params) result(ret)") call f%write(" use iso_fortran_env") call f%write(" real(real64),intent(in) :: CauchyStress(:,:),PlasticStrain(:,:),params(:)") call f%write(" real(real64) :: ret,k,J_2") call f%write(" k = params(1)") call f%write(" J_2 = to_J2(CauchyStress)") call f%write(" ret = J_2 - k**2") call f%write("end function") call f%close() end subroutine function get_Return_mapping_tangent_matrix(PlasticPotential, CauchyStress, PlasticStrain, Strain, & ElasticParams, gamma, PlasticParams, epsilon) result(ret) procedure(P_PotentialFunction) :: PlasticPotential real(real64), intent(in) :: CauchyStress(:, :), PlasticStrain(:, :), ElasticParams(:), PlasticParams(:), & gamma, Strain(:, :) integer(int32), allocatable :: indx(:, :) real(real64), optional, intent(in) :: epsilon real(real64), allocatable :: dsigma_tensor(:, :), ret(:, :) real(real64) :: dsigma, df integer(int32) :: i, j complex(real64) :: forward, backward real(real64), allocatable :: ds(:, :), dfds_forward(:, :), dfds_backward(:, :), dfds_vec(:), & dfds_backward_vec(:), dfds_forward_vec(:), E_dfds_backward(:), E_dfds_forward(:), & E_dfds_vec(:), dfds(:, :), X_vec(:), Y_vec(:), J_mat(:, :), tr_CauchyStress(:, :) type(Math_) :: math indx = zeros(6, 2) indx(1, 1:2) = [1, 1] indx(2, 1:2) = [2, 2] indx(3, 1:2) = [3, 3] indx(4, 1:2) = [1, 2] indx(5, 1:2) = [2, 3] indx(6, 1:2) = [1, 3] ret = zeros(6 + 1, 6 + 1) ds = zeros(3, 3) ret(1, 1) = 1.0d0 ret(2, 2) = 1.0d0 ret(3, 3) = 1.0d0 ret(4, 4) = 1.0d0 ret(5, 5) = 1.0d0 ret(6, 6) = 1.0d0 dfds_forward_vec = zeros(6) dfds_backward_vec = zeros(6) dfds_vec = zeros(6) do J = 1, 6 ! forward ds(:, :) = 0.0d0 ds(indx(J, 1), indx(J, 2)) = epsilon dfds_forward = d_dSigma(PlasticPotential=PlasticPotential, CauchyStress=CauchyStress + ds, & PlasticStrain=PlasticStrain, params=PlasticParams, epsilon=epsilon) dfds_backward = d_dSigma(PlasticPotential=PlasticPotential, CauchyStress=CauchyStress - ds, & PlasticStrain=PlasticStrain, params=PlasticParams, epsilon=epsilon) do I = 1, 6 dfds_forward_vec = dfds_forward(indx(I, 1), indx(I, 2)) dfds_backward_vec = dfds_backward(indx(I, 1), indx(I, 2)) end do do I = 1, 6 E_dfds_forward = gamma*matmul(StVenant_StiffnessMatrix_2D(params=ElasticParams), dfds_forward_vec) E_dfds_backward = gamma*matmul(StVenant_StiffnessMatrix_2D(params=ElasticParams), dfds_backward_vec) forward = E_dfds_forward(I) backward = E_dfds_backward(I) ret(I, J) = ret(I, J) + real(forward - backward)/(2.0d0*epsilon) end do I = 7 forward = PlasticPotential(CauchyStress=dcmplx(CauchyStress + ds*math%i), & PlasticStrain=dcmplx(PlasticStrain), & params=PlasticParams) backward = PlasticPotential(CauchyStress=dcmplx(CauchyStress - ds*math%i), & PlasticStrain=dcmplx(PlasticStrain), & params=PlasticParams) ret(I, J) = ret(I, J) + aimag((forward - backward)/(2.0d0*epsilon)) end do J = 7 dfds = d_dSigma(PlasticPotential=PlasticPotential, CauchyStress=CauchyStress, & PlasticStrain=PlasticStrain, params=PlasticParams, epsilon=epsilon) do I = 1, 6 dfds_vec = dfds_backward(indx(I, 1), indx(I, 2)) end do E_dfds_vec = matmul(StVenant_StiffnessMatrix_2D(params=ElasticParams), dfds_vec) ret(1:6, J) = E_dfds_vec(:) ret(7, 7) = 0.0d0 end function ! ################################################### subroutine initElastoPlasticity(this, femdomains, & default_YieldFunction, default_YieldFunction_params, & default_PlasticPotential, default_PlasticPotential_params) class(ElastoPlasticity_), intent(inout) :: this type(FEMDomain_), target, intent(in) :: femdomains(:) real(real64), intent(in) :: default_YieldFunction_params(:) real(real64), intent(in) :: default_PlasticPotential_params(:) procedure(P_PotentialFunction) :: default_YieldFunction procedure(P_PotentialFunction) :: default_PlasticPotential integer(int32) :: i, j, n, ne, ngp if (allocated(this%ep_domain)) deallocate (this%ep_domain) n = size(femdomains) allocate (this%ep_domain(n)) do i = 1, n this%ep_domain(i)%femdomain = femdomains(i) this%ep_domain(i)%YieldFunction => default_YieldFunction this%ep_domain(i)%PlasticPotential => default_PlasticPotential this%ep_domain(i)%YieldFunction_params & = zeros(femdomains(i)%ne(), size(default_YieldFunction_params)) do j = 1, femdomains(i)%ne() this%ep_domain(i)%YieldFunction_params(j, :) = default_YieldFunction_params(:) end do this%ep_domain(i)%PlasticPotential_params & = zeros(femdomains(i)%ne(), size(default_PlasticPotential_params)) do j = 1, femdomains(i)%ne() this%ep_domain(i)%PlasticPotential_params(j, :) = default_PlasticPotential_params(:) end do ngp = this%ep_domain(i)%femdomain%ngp() ne = this%ep_domain(i)%femdomain%ne() this%ep_domain(i)%ElasticPotential_params = zeros(ne, 2) ! E, v this%ep_domain(i)%CauchyStress_field = zeros(6, ngp, ne) this%ep_domain(i)%Strain_field = zeros(6, ngp, ne) this%ep_domain(i)%PlasticStrain_field = zeros(6, ngp, ne) this%ep_domain(i)%displacement = & zeros(this%ep_domain(i)%femdomain%nn()*this%ep_domain(i)%femdomain%nd()) end do end subroutine ! ################################################### ! ################################################### subroutine edit_YF_PP_ElastoPlasticity(this, DomainID, YieldFunction, PlasticPotential) class(ElastoPlasticity_), intent(inout) :: this integer(int32), intent(in) :: DomainID procedure(P_PotentialFunction) :: YieldFunction procedure(P_PotentialFunction) :: PlasticPotential if (associated(this%ep_domain(DomainID)%YieldFunction)) nullify (this%ep_domain(DomainID)%YieldFunction) this%ep_domain(DomainID)%YieldFunction => YieldFunction if (associated(this%ep_domain(DomainID)%PlasticPotential)) nullify (this%ep_domain(DomainID)%PlasticPotential) this%ep_domain(DomainID)%PlasticPotential => PlasticPotential end subroutine ! ################################################### subroutine solveElastoPlasticity(this, & YoungModulus, PoissonRatio, Density, & fix_node_list_x, & fix_node_list_y, & fix_node_list_z, & fix_value_list_x, & fix_value_list_y, & fix_value_list_z) class(ElastoPlasticity_), intent(inout) :: this real(real64), intent(in) :: YoungModulus(:) real(real64), intent(in) :: PoissonRatio(:) real(real64), intent(in) :: Density(:) integer(int32), optional, intent(in) :: fix_node_list_x(:) integer(int32), optional, intent(in) :: fix_node_list_y(:) integer(int32), optional, intent(in) :: fix_node_list_z(:) real(real64), optional, intent(in) :: fix_value_list_x(:) real(real64), optional, intent(in) :: fix_value_list_y(:) real(real64), optional, intent(in) :: fix_value_list_z(:) real(real64), allocatable :: disp_tr(:), d_disp(:), & F_int(:), & F_ext(:), & Residual_vector(:) integer(int32) :: newton_loop_itr integer(int32) :: ElementID type(CRS_) :: K_matrix if (size(this%ep_domain) == 1) then this%ep_domain(1)%ElasticPotential_params(:, 1) = YoungModulus(:)*PoissonRatio(:)/(1.0d0 + PoissonRatio(:)) & /(1.0d0 - 2.0d0*PoissonRatio(:)) this%ep_domain(1)%ElasticPotential_params(:, 2) = YoungModulus(:)/2.0d0/(1.0d0 + PoissonRatio(:)) if (.not. this%femsolver%initialized) then call this%femsolver%init(NumDomain=1) call this%femsolver%setDomain(FEMDomain=this%ep_domain(1)%femdomain, DomainID=1) call this%femsolver%setCRS(DOF=3) else call this%femsolver%zeros() end if !$OMP parallel !$OMP do do ElementID = 1, this%ep_domain(1)%femdomain%ne() call this%femsolver%setMatrix(DomainID=1, ElementID=ElementID, DOF=3, & Matrix=this%ep_domain(1)%femdomain%StiffnessMatrix(ElementID=ElementID, & E=YoungModulus(ElementID), & v=PoissonRatio(ElementID))) call this%femsolver%setVector(DomainID=1, ElementID=ElementID, DOF=3, & Vector=this%ep_domain(1)%femdomain%MassVector( & ElementID=ElementID, & DOF=this%ep_domain(1)%femdomain%nd(), & Density=Density(ElementID), & Accel=this%gravity_accel & ) & ) end do !$OMP end do !$OMP end parallel call this%femsolver%fix(DomainID=1, IDs=fix_node_list_x*3 - 2, FixValues=fix_value_list_x) call this%femsolver%fix(DomainID=1, IDs=fix_node_list_y*3 - 1, FixValues=fix_value_list_y) call this%femsolver%fix(DomainID=1, IDs=fix_node_list_z*3 - 0, FixValues=fix_value_list_z) disp_tr = this%femsolver%solve() K_matrix = this%femsolver%getCRS() d_disp = disp_tr print *, "[ok] trial disp done!" ! perform modified Newton-Raphson method do newton_loop_itr = 1, this%MAX_NEWTON_LOOP_ITR !F_int = this%get_internal_force(dU=d_disp) ! \int_{\Omega} B \sigma d \Omega !dU -> d_eps -> F_int = this%get_internal_force(dU=d_disp) ! \int_{\Omega} B \sigma d \Omega F_ext = this%get_external_force() ! Traction force F_ext = 0.0d0 Residual_vector = F_ext - F_int if (norm(Residual_vector) == 0.0d0) then exit end if call this%fill_zero_at_DBC(values=Residual_vector, & idx=(fix_node_list_x*3 - 2) & //(fix_node_list_y*3 - 1)//(fix_node_list_z*3 - 0)) ! solve [K]{du} = {R} call K_matrix%BiCGSTAB(x=d_disp, b=Residual_vector) disp_tr = disp_tr - d_disp print *, newton_loop_itr, norm(d_disp), norm(Residual_vector), norm(F_int) if (norm(d_disp) < this%TOL) then exit end if end do this%ep_domain(1)%displacement = disp_tr else print *, "[ERROR] solveElastoPlasticity >> only for single-domain problem" stop end if end subroutine ! ################################################### ! ################################################### subroutine solve_increment_ElastoPlasticity(this, & YoungModulus, PoissonRatio, delta_Density, & fix_node_list_x, & fix_node_list_y, & fix_node_list_z, & fix_value_list_x, & fix_value_list_y, & fix_value_list_z) class(ElastoPlasticity_), intent(inout) :: this real(real64), intent(in) :: YoungModulus(:) real(real64), intent(in) :: PoissonRatio(:) real(real64), intent(in) :: delta_Density(:) integer(int32), optional, intent(in) :: fix_node_list_x(:) integer(int32), optional, intent(in) :: fix_node_list_y(:) integer(int32), optional, intent(in) :: fix_node_list_z(:) real(real64), optional, intent(in) :: fix_value_list_x(:) real(real64), optional, intent(in) :: fix_value_list_y(:) real(real64), optional, intent(in) :: fix_value_list_z(:) real(real64), allocatable :: d_disp_tr(:), d_d_disp(:), & dF_int(:), & dF_ext(:), & Residual_vector(:) integer(int32) :: newton_loop_itr integer(int32) :: ElementID, DomainID type(CRS_) :: K_matrix DomainID = 1 ! only for single-domain if (.not. allocated(this%ep_domain(domainID)%PlasticStrain_field_n)) then this%ep_domain(domainID)%PlasticStrain_field_n = this%ep_domain(domainID)%PlasticStrain_field end if if (.not. allocated(this%ep_domain(DomainID)%dCauchyStress_field)) then this%ep_domain(DomainID)%dCauchyStress_field = this%ep_domain(DomainID)%CauchyStress_field end if if (.not. allocated(this%ep_domain(DomainID)%dStrain_field)) then this%ep_domain(DomainID)%dStrain_field = this%ep_domain(DomainID)%Strain_field end if this%ep_domain(domainID)%PlasticStrain_field = this%ep_domain(domainID)%PlasticStrain_field_n this%ep_domain(DomainID)%dCauchyStress_field(:, :, :) = 0.0d0 this%ep_domain(DomainID)%dStrain_field(:, :, :) = 0.0d0 if (size(this%ep_domain) == 1) then this%ep_domain(1)%ElasticPotential_params(:, 1) = YoungModulus(:)*PoissonRatio(:)/(1.0d0 + PoissonRatio(:)) & /(1.0d0 - 2.0d0*PoissonRatio(:)) this%ep_domain(1)%ElasticPotential_params(:, 2) = YoungModulus(:)/2.0d0/(1.0d0 + PoissonRatio(:)) if (.not. this%femsolver%initialized) then call this%femsolver%init(NumDomain=1) call this%femsolver%setDomain(FEMDomain=this%ep_domain(1)%femdomain, DomainID=1) call this%femsolver%setCRS(DOF=3) else call this%femsolver%zeros() end if do ElementID = 1, this%ep_domain(1)%femdomain%ne() call this%femsolver%setMatrix(DomainID=1, ElementID=ElementID, DOF=3, & Matrix=this%ep_domain(1)%femdomain%StiffnessMatrix(ElementID=ElementID, & E=YoungModulus(ElementID), & v=PoissonRatio(ElementID))) call this%femsolver%setVector(DomainID=1, ElementID=ElementID, DOF=3, & Vector=this%ep_domain(1)%femdomain%MassVector( & ElementID=ElementID, & DOF=this%ep_domain(1)%femdomain%nd(), & Density=delta_Density(ElementID), & Accel=this%gravity_accel & ) & ) end do call this%femsolver%fix(DomainID=1, IDs=fix_node_list_x*3 - 2, FixValues=fix_value_list_x) call this%femsolver%fix(DomainID=1, IDs=fix_node_list_y*3 - 1, FixValues=fix_value_list_y) call this%femsolver%fix(DomainID=1, IDs=fix_node_list_z*3 - 0, FixValues=fix_value_list_z) d_disp_tr = this%femsolver%solve() K_matrix = this%femsolver%getCRS() d_d_disp = d_disp_tr dF_int = zeros(size(d_disp_tr)) print *, "[ok] trial delta-disp done!" ! perform modified Newton-Raphson method do newton_loop_itr = 1, this%MAX_NEWTON_LOOP_ITR !dF_int = this%get_internal_force(dU=d_d_disp) ! \int_{\Omega} B \sigma d \Omega !dU -> d_eps -> call this%update_stress_for_increment(dU=d_d_disp) ! \int_{\Omega} B \sigma d \Omega dF_int = this%get_delta_internal_force() dF_ext = this%get_external_force() ! Traction force dF_ext = 0.0d0 call this%fill_zero_at_DBC(values=dF_ext, & idx=(fix_node_list_x*3 - 2) & //(fix_node_list_y*3 - 1)//(fix_node_list_z*3 - 0)) !print *, "norm(dF_ext)",norm(dF_ext) Residual_vector = dF_ext - dF_int if (norm(Residual_vector) == 0.0d0) then exit end if call this%fill_zero_at_DBC(values=Residual_vector, & idx=(fix_node_list_x*3 - 2) & //(fix_node_list_y*3 - 1)//(fix_node_list_z*3 - 0)) ! solve [K]{du} = {R} call K_matrix%BiCGSTAB(x=d_d_disp, b=Residual_vector) d_disp_tr = d_disp_tr - d_d_disp print *, newton_loop_itr, norm(d_d_disp), norm(Residual_vector), norm(dF_int) if (norm(d_d_disp) < this%TOL) then exit end if end do this%ep_domain(1)%displacement = this%ep_domain(1)%displacement + d_disp_tr else print *, "[ERROR] solveElastoPlasticity >> only for single-domain problem" stop end if this%ep_domain(DomainID)%Strain_field = this%ep_domain(DomainID)%Strain_field & + this%ep_domain(DomainID)%dStrain_field this%ep_domain(domainID)%PlasticStrain_field_n = & this%ep_domain(domainID)%PlasticStrain_field this%ep_domain(DomainID)%CauchyStress_field = & this%ep_domain(DomainID)%CauchyStress_field & + this%ep_domain(DomainID)%dCauchyStress_field this%ep_domain(DomainID)%dCauchyStress_field(:, :, :) = 0.0d0 this%ep_domain(DomainID)%dStrain_field(:, :, :) = 0.0d0 end subroutine ! ################################################### function get_internal_forceElastoPlasticity(this, dU) result(ret) class(ElastoPlasticity_), intent(inout) :: this real(real64), intent(in) :: dU(:) real(real64), allocatable :: ret(:), PlasticStrain(:, :), dStrain(:, :), Te(:), & StressVector(:), Bmat(:, :), CauchyStress(:, :), dU_mat(:, :) type(ShapeFunction_) :: sf integer(int32) :: ElementID, GaussPointID, num_node, i, j num_node = size(dU)/3 dU_mat = zeros(num_node, 3) do i = 1, num_node do j = 1, 3 dU_mat(i, j) = dU((i - 1)*3 + j) end do end do if (size(this%ep_domain) == 1) then if (norm(dU) /= 0.0d0) then ! Update stress and plastic strain !$OMP parallel default(shared) private(GaussPointID,PlasticStrain,dStrain,CauchyStress) !$OMP do do ElementID = 1, this%ep_domain(1)%femdomain%ne() do GaussPointID = 1, this%ep_domain(1)%femdomain%ngp() PlasticStrain = this%getPlasticStrain(ElementID=ElementID, GaussPointID=GaussPointID) dStrain = this%ep_domain(1)%femdomain%getStrainTensor( & ElementID=ElementID, GaussPointID=GaussPointID, & displacement=dU_mat) !print *, minval(dU_mat),maxval(dU_mat) CauchyStress = to_StressTensor( & YieldFunction=this%ep_domain(1)%YieldFunction, & PlasticPotential=this%ep_domain(1)%PlasticPotential, & Strain=this%getStrain(ElementID=ElementID, GaussPointID=GaussPointID), & dStrain=dStrain, & CauchyStress=this%getCauchyStress(ElementID=ElementID, GaussPointID=GaussPointID), & PlasticStrain=PlasticStrain, & YieldParams=this%ep_domain(1)%YieldFunction_params(ElementID, :), & PlasticParams=this%ep_domain(1)%PlasticPotential_params(ElementID, :), & ElasticParams=this%ep_domain(1)%ElasticPotential_params(ElementID, :), & epsilon=dble(1.0e-4) & ) call this%setPlasticStrain( & ElementID=ElementID, & GaussPointID=GaussPointID, & PlasticStrain=PlasticStrain) call this%setCauchyStress( & ElementID=ElementID, & GaussPointID=GaussPointID, & CauchyStress=CauchyStress) call this%addStrain( & ElementID=ElementID, & GaussPointID=GaussPointID, & dStrain=dStrain) end do end do !$OMP end do !$OMP end parallel end if ret = zeros(this%ep_domain(1)%femdomain%nn()*this%ep_domain(1)%femdomain%nd()) !$OMP parallel default(shared) private(GaussPointID,sf,Bmat,StressVector,Te) !$OMP do reduction(+:ret) ! Perform gauss integral to compute ret(:) do ElementID = 1, this%ep_domain(1)%femdomain%ne() do GaussPointID = 1, this%ep_domain(1)%femdomain%ngp() ! get shape function sf = this%ep_domain(1)%FEMdomain%getShapeFunction( & ElementID=ElementID, GaussPointID=GaussPointID) ! get B-matrix Bmat = this%ep_domain(1)%FEMdomain%BMatrix( & shapefunction=sf, ElementID=ElementID) ! get Stress vector StressVector = this%ep_domain(1)%CauchyStress_field(:, GaussPointID, ElementID) Te = matmul(transpose(Bmat), StressVector)*sf%detJ ret = ret + & this%ep_domain(1)%FEMdomain%asGlobalVector(LocalVector=Te, ElementID=ElementID, & DOF=this%ep_domain(1)%FEMdomain%nd()) end do end do !$OMP end do !$OMP end parallel else print *, "[ERROR] get_internal_forceElastoPlasticity >> only for single-domain problem" stop end if end function ! ################################################### ! ################################################### ! ################################################### subroutine update_stress_for_incrementElastoPlasticity(this, dU) class(ElastoPlasticity_), intent(inout) :: this real(real64), intent(in) :: dU(:) real(real64), allocatable :: ret(:), PlasticStrain(:, :), PlasticStrain_n(:, :), dStrain(:, :), dTe(:), & dStressVector(:), Bmat(:, :), dCauchyStress(:, :), CauchyStress(:, :), dU_mat(:, :), new_PlasticStrain(:, :) type(ShapeFunction_) :: sf integer(int32) :: ElementID, GaussPointID, num_node, i, j num_node = size(dU)/3 dU_mat = zeros(num_node, 3) do i = 1, num_node do j = 1, 3 dU_mat(i, j) = dU((i - 1)*3 + j) end do end do if (size(this%ep_domain) == 1) then if (norm(dU) /= 0.0d0) then ! Update stress and plastic strain !$OMP parallel default(shared) private(GaussPointID,PlasticStrain_n,new_PlasticStrain,dStrain,dCauchyStress) !$OMP do do ElementID = 1, this%ep_domain(1)%femdomain%ne() do GaussPointID = 1, this%ep_domain(1)%femdomain%ngp() PlasticStrain_n = this%getPlasticStrain_n(ElementID=ElementID, GaussPointID=GaussPointID) dStrain = this%ep_domain(1)%femdomain%getStrainTensor( & ElementID=ElementID, GaussPointID=GaussPointID, & displacement=dU_mat) !print *, minval(dU_mat),maxval(dU_mat) dCauchyStress = to_dStressTensor( & YieldFunction=this%ep_domain(1)%YieldFunction, & PlasticPotential=this%ep_domain(1)%PlasticPotential, & dStrain=dStrain, & CauchyStress=this%getCauchyStress(ElementID=ElementID, GaussPointID=GaussPointID) & + this%getdCauchyStress(ElementID=ElementID, GaussPointID=GaussPointID), & PlasticStrain=PlasticStrain_n, & YieldParams=this%ep_domain(1)%YieldFunction_params(ElementID, :), & PlasticParams=this%ep_domain(1)%PlasticPotential_params(ElementID, :), & ElasticParams=this%ep_domain(1)%ElasticPotential_params(ElementID, :), & epsilon=dble(1.0e-4), & new_PlasticStrain=new_PlasticStrain & ) call this%setPlasticStrain( & ElementID=ElementID, & GaussPointID=GaussPointID, & PlasticStrain=new_PlasticStrain) !call this%setdCauchyStress(& ! ElementID=ElementID,& ! GaussPointID=GaussPointID,& ! dCauchyStress=& ! this%getdCauchyStress(ElementID=ElementID,GaussPointID=GaussPointID) & ! + dCauchyStress) call this%setdCauchyStress( & ElementID=ElementID, & GaussPointID=GaussPointID, & dCauchyStress= & this%getdCauchyStress(ElementID=ElementID, GaussPointID=GaussPointID) & + dCauchyStress) call this%setdStrain( & ElementID=ElementID, & GaussPointID=GaussPointID, & dStrain=dStrain) end do end do !$OMP end do !$OMP end parallel end if else print *, "[ERROR] get_internal_forceElastoPlasticity >> only for single-domain problem" stop end if end subroutine ! ################################################### function get_delta_internal_forceElastoPlasticity(this) result(ret) class(ElastoPlasticity_), intent(inout) :: this real(real64), allocatable :: ret(:), PlasticStrain(:, :), PlasticStrain_n(:, :), dStrain(:, :), dTe(:), & dStressVector(:), Bmat(:, :), dCauchyStress(:, :), CauchyStress(:, :), dU_mat(:, :), new_PlasticStrain(:, :) type(ShapeFunction_) :: sf integer(int32) :: ElementID, GaussPointID, num_node, i, j if (size(this%ep_domain) == 1) then ret = zeros(this%ep_domain(1)%femdomain%nn()*this%ep_domain(1)%femdomain%nd()) !$OMP parallel default(shared) private(GaussPointID,sf,Bmat,dStressVector,dTe) !$OMP do reduction(+:ret) ! Perform gauss integral to compute ret(:) do ElementID = 1, this%ep_domain(1)%femdomain%ne() do GaussPointID = 1, this%ep_domain(1)%femdomain%ngp() ! get shape function sf = this%ep_domain(1)%FEMdomain%getShapeFunction( & ElementID=ElementID, GaussPointID=GaussPointID) ! get B-matrix Bmat = this%ep_domain(1)%FEMdomain%BMatrix( & shapefunction=sf, ElementID=ElementID) ! get Stress vector dStressVector = this%ep_domain(1)%dCauchyStress_field(:, GaussPointID, ElementID) dTe = matmul(transpose(Bmat), dStressVector)*sf%detJ ret = ret + & this%ep_domain(1)%FEMdomain%asGlobalVector(LocalVector=dTe, ElementID=ElementID, & DOF=this%ep_domain(1)%FEMdomain%nd()) end do end do !$OMP end do !$OMP end parallel else print *, "[ERROR] get_internal_forceElastoPlasticity >> only for single-domain problem" stop end if end function ! ################################################### ! ################################################## function getPlasticStrain_ElastoPlasticity(this, GaussPointID, ElementID) result(ret) class(ElastoPlasticity_), intent(in) :: this integer(int32), intent(in) :: GaussPointID, ElementID real(real64), allocatable :: ret(:, :), ret_vec(:) ret_vec = this%ep_domain(1)%PlasticStrain_field(:, GaussPointID, ElementID) ret = zeros(3, 3) ret(1, 1) = ret_vec(1) ret(2, 2) = ret_vec(2) ret(3, 3) = ret_vec(3) ret(1, 2) = ret_vec(4) ret(2, 3) = ret_vec(5) ret(1, 3) = ret_vec(6) ret(2, 1) = ret_vec(4) ret(3, 2) = ret_vec(5) ret(3, 1) = ret_vec(6) end function ! ################################################## ! ################################################## function getPlasticStrain_n_ElastoPlasticity(this, GaussPointID, ElementID) result(ret) class(ElastoPlasticity_), intent(in) :: this integer(int32), intent(in) :: GaussPointID, ElementID real(real64), allocatable :: ret(:, :), ret_vec(:) ret_vec = this%ep_domain(1)%PlasticStrain_field_n(:, GaussPointID, ElementID) ret = zeros(3, 3) ret(1, 1) = ret_vec(1) ret(2, 2) = ret_vec(2) ret(3, 3) = ret_vec(3) ret(1, 2) = ret_vec(4) ret(2, 3) = ret_vec(5) ret(1, 3) = ret_vec(6) ret(2, 1) = ret_vec(4) ret(3, 2) = ret_vec(5) ret(3, 1) = ret_vec(6) end function ! ################################################## ! ################################################## function getStrain_ElastoPlasticity(this, GaussPointID, ElementID) result(ret) class(ElastoPlasticity_), intent(in) :: this integer(int32), intent(in) :: GaussPointID, ElementID real(real64), allocatable :: ret(:, :), ret_vec(:) ret_vec = this%ep_domain(1)%Strain_field(:, GaussPointID, ElementID) ret = zeros(3, 3) ret(1, 1) = ret_vec(1) ret(2, 2) = ret_vec(2) ret(3, 3) = ret_vec(3) ret(1, 2) = ret_vec(4) ret(2, 3) = ret_vec(5) ret(1, 3) = ret_vec(6) ret(2, 1) = ret_vec(4) ret(3, 2) = ret_vec(5) ret(3, 1) = ret_vec(6) end function ! ################################################## ! ################################################## function getCauchyStress_ElastoPlasticity(this, GaussPointID, ElementID) result(ret) class(ElastoPlasticity_), intent(in) :: this integer(int32), intent(in) :: GaussPointID, ElementID real(real64), allocatable :: ret(:, :), ret_vec(:) ret_vec = this%ep_domain(1)%CauchyStress_field(:, GaussPointID, ElementID) ret = zeros(3, 3) ret(1, 1) = ret_vec(1) ret(2, 2) = ret_vec(2) ret(3, 3) = ret_vec(3) ret(1, 2) = ret_vec(4) ret(2, 3) = ret_vec(5) ret(1, 3) = ret_vec(6) ret(2, 1) = ret_vec(4) ret(3, 2) = ret_vec(5) ret(3, 1) = ret_vec(6) end function ! ################################################## ! ################################################## function getdCauchyStress_ElastoPlasticity(this, GaussPointID, ElementID) result(ret) class(ElastoPlasticity_), intent(in) :: this integer(int32), intent(in) :: GaussPointID, ElementID real(real64), allocatable :: ret(:, :), ret_vec(:) ret_vec = this%ep_domain(1)%dCauchyStress_field(:, GaussPointID, ElementID) ret = zeros(3, 3) ret(1, 1) = ret_vec(1) ret(2, 2) = ret_vec(2) ret(3, 3) = ret_vec(3) ret(1, 2) = ret_vec(4) ret(2, 3) = ret_vec(5) ret(1, 3) = ret_vec(6) ret(2, 1) = ret_vec(4) ret(3, 2) = ret_vec(5) ret(3, 1) = ret_vec(6) end function ! ################################################## subroutine setPlasticStrain_ElastoPlasticity(this, ElementID, GaussPointID, PlasticStrain) class(ElastoPlasticity_), intent(inout) :: this integer(int32), intent(in) :: ElementID, GaussPointID real(real64), intent(in) :: PlasticStrain(:, :) this%ep_domain(1)%PlasticStrain_field(1, GaussPointID, ElementID) = PlasticStrain(1, 1) this%ep_domain(1)%PlasticStrain_field(2, GaussPointID, ElementID) = PlasticStrain(2, 2) this%ep_domain(1)%PlasticStrain_field(3, GaussPointID, ElementID) = PlasticStrain(3, 3) this%ep_domain(1)%PlasticStrain_field(4, GaussPointID, ElementID) = PlasticStrain(1, 2) this%ep_domain(1)%PlasticStrain_field(5, GaussPointID, ElementID) = PlasticStrain(2, 3) this%ep_domain(1)%PlasticStrain_field(6, GaussPointID, ElementID) = PlasticStrain(1, 3) end subroutine subroutine setCauchyStress_ElastoPlasticity(this, ElementID, GaussPointID, CauchyStress) class(ElastoPlasticity_), intent(inout) :: this integer(int32), intent(in) :: ElementID, GaussPointID real(real64), intent(in) :: CauchyStress(:, :) this%ep_domain(1)%CauchyStress_field(1, GaussPointID, ElementID) = CauchyStress(1, 1) this%ep_domain(1)%CauchyStress_field(2, GaussPointID, ElementID) = CauchyStress(2, 2) this%ep_domain(1)%CauchyStress_field(3, GaussPointID, ElementID) = CauchyStress(3, 3) this%ep_domain(1)%CauchyStress_field(4, GaussPointID, ElementID) = CauchyStress(1, 2) this%ep_domain(1)%CauchyStress_field(5, GaussPointID, ElementID) = CauchyStress(2, 3) this%ep_domain(1)%CauchyStress_field(6, GaussPointID, ElementID) = CauchyStress(1, 3) end subroutine subroutine setdCauchyStress_ElastoPlasticity(this, ElementID, GaussPointID, dCauchyStress) class(ElastoPlasticity_), intent(inout) :: this integer(int32), intent(in) :: ElementID, GaussPointID real(real64), intent(in) :: dCauchyStress(:, :) this%ep_domain(1)%dCauchyStress_field(1, GaussPointID, ElementID) = dCauchyStress(1, 1) this%ep_domain(1)%dCauchyStress_field(2, GaussPointID, ElementID) = dCauchyStress(2, 2) this%ep_domain(1)%dCauchyStress_field(3, GaussPointID, ElementID) = dCauchyStress(3, 3) this%ep_domain(1)%dCauchyStress_field(4, GaussPointID, ElementID) = dCauchyStress(1, 2) this%ep_domain(1)%dCauchyStress_field(5, GaussPointID, ElementID) = dCauchyStress(2, 3) this%ep_domain(1)%dCauchyStress_field(6, GaussPointID, ElementID) = dCauchyStress(1, 3) end subroutine subroutine addStrain_ElastoPlasticity(this, ElementID, GaussPointID, dStrain) class(ElastoPlasticity_), intent(inout) :: this integer(int32), intent(in) :: ElementID, GaussPointID real(real64), intent(in) :: dStrain(:, :) this%ep_domain(1)%Strain_field(1, GaussPointID, ElementID) = this%ep_domain(1)%Strain_field(1, GaussPointID, ElementID) & + dStrain(1, 1) this%ep_domain(1)%Strain_field(2, GaussPointID, ElementID) = this%ep_domain(1)%Strain_field(2, GaussPointID, ElementID) & + dStrain(2, 2) this%ep_domain(1)%Strain_field(3, GaussPointID, ElementID) = this%ep_domain(1)%Strain_field(3, GaussPointID, ElementID) & + dStrain(3, 3) this%ep_domain(1)%Strain_field(4, GaussPointID, ElementID) = this%ep_domain(1)%Strain_field(4, GaussPointID, ElementID) & + dStrain(1, 2) this%ep_domain(1)%Strain_field(5, GaussPointID, ElementID) = this%ep_domain(1)%Strain_field(5, GaussPointID, ElementID) & + dStrain(2, 3) this%ep_domain(1)%Strain_field(6, GaussPointID, ElementID) = this%ep_domain(1)%Strain_field(6, GaussPointID, ElementID) & + dStrain(1, 3) end subroutine subroutine setdStrain_ElastoPlasticity(this, ElementID, GaussPointID, dStrain) class(ElastoPlasticity_), intent(inout) :: this integer(int32), intent(in) :: ElementID, GaussPointID real(real64), intent(in) :: dStrain(:, :) this%ep_domain(1)%dStrain_field(1, GaussPointID, ElementID) = dStrain(1, 1) this%ep_domain(1)%dStrain_field(2, GaussPointID, ElementID) = dStrain(2, 2) this%ep_domain(1)%dStrain_field(3, GaussPointID, ElementID) = dStrain(3, 3) this%ep_domain(1)%dStrain_field(4, GaussPointID, ElementID) = dStrain(1, 2) this%ep_domain(1)%dStrain_field(5, GaussPointID, ElementID) = dStrain(2, 3) this%ep_domain(1)%dStrain_field(6, GaussPointID, ElementID) = dStrain(1, 3) end subroutine ! ################################################### function get_external_forceElastoPlasticity(this) result(ret) class(ElastoPlasticity_), intent(inout) :: this real(real64), allocatable :: ret(:) ret = this%femsolver%CRS_RHS(:) end function ! ################################################### ! ################################################### subroutine fill_zero_at_DBCElastoPlasticity(this, values, idx) class(ElastoPlasticity_), intent(inout) :: this real(real64), intent(inout) :: values(:) integer(int32), intent(in) :: idx(:) values(idx(:)) = 0.0d0 end subroutine ! ################################################### ! ################################################### subroutine exportElastoPlasticity(this, name, step, amp) class(ElastoPlasticity_), intent(inout) :: this character(*), intent(in) :: name real(real64), optional, intent(in) :: amp real(real64) :: mag integer(int32), optional, intent(in) :: step integer(int32) :: i, n mag = input(default=1.0d0, option=amp) do i = 1, size(this%ep_domain) call this%ep_domain(i)%femdomain%deform(disp=mag*this%ep_domain(i)%displacement) end do if (present(step)) then do i = 1, size(this%ep_domain) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_step_"//zfill(step, 5)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_s11"//"_step_"//zfill(step, 5), & scalar=this%ep_domain(1)%CauchyStress_field(1, 1, :)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_s22"//"_step_"//zfill(step, 5), & scalar=this%ep_domain(1)%CauchyStress_field(2, 1, :)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_s33"//"_step_"//zfill(step, 5), & scalar=this%ep_domain(1)%CauchyStress_field(3, 1, :)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_s12"//"_step_"//zfill(step, 5), & scalar=this%ep_domain(1)%CauchyStress_field(4, 1, :)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_s23"//"_step_"//zfill(step, 5), & scalar=this%ep_domain(1)%CauchyStress_field(5, 1, :)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_s13"//"_step_"//zfill(step, 5), & scalar=this%ep_domain(1)%CauchyStress_field(6, 1, :)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_I1"//"_step_"//zfill(step, 5), & scalar=this%I1()) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_J2"//"_step_"//zfill(step, 5), & scalar=this%J2()) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_eI1"//"_step_"//zfill(step, 5), & scalar=this%I1_e()) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_eJ2"//"_step_"//zfill(step, 5), & scalar=this%J2_e()) end do else do i = 1, size(this%ep_domain) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & ) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_s11", scalar=this%ep_domain(1)%CauchyStress_field(1, 1, :)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_s22", scalar=this%ep_domain(1)%CauchyStress_field(2, 1, :)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_s33", scalar=this%ep_domain(1)%CauchyStress_field(3, 1, :)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_s12", scalar=this%ep_domain(1)%CauchyStress_field(4, 1, :)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_s23", scalar=this%ep_domain(1)%CauchyStress_field(5, 1, :)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_s13", scalar=this%ep_domain(1)%CauchyStress_field(6, 1, :)) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_I1", scalar=this%I1()) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_J2", scalar=this%J2()) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_eI1", scalar=this%I1_e()) call this%ep_domain(i)%femdomain%vtk(name//"_domain_"//zfill(i, 5) & //"_eJ2", scalar=this%J2_e()) end do end if call this%exportField(name=name) do i = 1, size(this%ep_domain) call this%ep_domain(i)%femdomain%deform(disp=-mag*this%ep_domain(i)%displacement) end do end subroutine ! ################################################### ! ################################################### function I1ElastoPlasticity(this) result(I1) class(ElastoPlasticity_), intent(inout) :: this real(real64), allocatable :: I1(:), stress(:, :) integer(int32) :: i, j Stress = zeros(3, 3) I1 = zeros(this%ep_domain(1)%femdomain%ne()) do i = 1, this%ep_domain(1)%femdomain%ne() stress = 0.0d0 do j = 1, this%ep_domain(1)%femdomain%ngp() stress(1, 1) = stress(1, 1) + this%ep_domain(1)%CauchyStress_field(1, j, i) stress(2, 2) = stress(2, 2) + this%ep_domain(1)%CauchyStress_field(2, j, i) stress(3, 3) = stress(3, 3) + this%ep_domain(1)%CauchyStress_field(3, j, i) stress(1, 2) = stress(1, 2) + this%ep_domain(1)%CauchyStress_field(4, j, i) stress(2, 3) = stress(2, 3) + this%ep_domain(1)%CauchyStress_field(5, j, i) stress(1, 3) = stress(1, 3) + this%ep_domain(1)%CauchyStress_field(6, j, i) end do I1(i) = to_I1(stress) end do end function ! ################################################### function J2ElastoPlasticity(this) result(J2) class(ElastoPlasticity_), intent(inout) :: this real(real64), allocatable :: J2(:), stress(:, :) integer(int32) :: i, j Stress = zeros(3, 3) J2 = zeros(this%ep_domain(1)%femdomain%ne()) do i = 1, this%ep_domain(1)%femdomain%ne() stress = 0.0d0 do j = 1, this%ep_domain(1)%femdomain%ngp() stress(1, 1) = stress(1, 1) + this%ep_domain(1)%CauchyStress_field(1, j, i) stress(2, 2) = stress(2, 2) + this%ep_domain(1)%CauchyStress_field(2, j, i) stress(3, 3) = stress(3, 3) + this%ep_domain(1)%CauchyStress_field(3, j, i) stress(1, 2) = stress(1, 2) + this%ep_domain(1)%CauchyStress_field(4, j, i) stress(2, 3) = stress(2, 3) + this%ep_domain(1)%CauchyStress_field(5, j, i) stress(1, 3) = stress(1, 3) + this%ep_domain(1)%CauchyStress_field(6, j, i) end do J2(i) = to_J2(stress) end do end function ! ################################################### function I1_e_ElastoPlasticity(this) result(I1) class(ElastoPlasticity_), intent(inout) :: this real(real64), allocatable :: I1(:), Strain(:, :) integer(int32) :: i, j Strain = zeros(3, 3) I1 = zeros(this%ep_domain(1)%femdomain%ne()) do i = 1, this%ep_domain(1)%femdomain%ne() Strain = 0.0d0 do j = 1, this%ep_domain(1)%femdomain%ngp() Strain(1, 1) = Strain(1, 1) + this%ep_domain(1)%Strain_field(1, j, i) Strain(2, 2) = Strain(2, 2) + this%ep_domain(1)%Strain_field(2, j, i) Strain(3, 3) = Strain(3, 3) + this%ep_domain(1)%Strain_field(3, j, i) Strain(1, 2) = Strain(1, 2) + this%ep_domain(1)%Strain_field(4, j, i) Strain(2, 3) = Strain(2, 3) + this%ep_domain(1)%Strain_field(5, j, i) Strain(1, 3) = Strain(1, 3) + this%ep_domain(1)%Strain_field(6, j, i) end do I1(i) = to_I1(Strain) end do end function ! ################################################### function J2_e_ElastoPlasticity(this) result(J2) class(ElastoPlasticity_), intent(inout) :: this real(real64), allocatable :: J2(:), Strain(:, :) integer(int32) :: i, j Strain = zeros(3, 3) J2 = zeros(this%ep_domain(1)%femdomain%ne()) do i = 1, this%ep_domain(1)%femdomain%ne() Strain = 0.0d0 do j = 1, this%ep_domain(1)%femdomain%ngp() Strain(1, 1) = Strain(1, 1) + this%ep_domain(1)%Strain_field(1, j, i) Strain(2, 2) = Strain(2, 2) + this%ep_domain(1)%Strain_field(2, j, i) Strain(3, 3) = Strain(3, 3) + this%ep_domain(1)%Strain_field(3, j, i) Strain(1, 2) = Strain(1, 2) + this%ep_domain(1)%Strain_field(4, j, i) Strain(2, 3) = Strain(2, 3) + this%ep_domain(1)%Strain_field(5, j, i) Strain(1, 3) = Strain(1, 3) + this%ep_domain(1)%Strain_field(6, j, i) end do J2(i) = to_J2(Strain) end do end function ! ################################################## function getTractionForce_ElastoPlasticity(this, NodeList) result(ret) class(ElastoPlasticity_), intent(inout) :: this integer(int32), intent(in) :: NodeList(:) real(real64), allocatable :: t(:) real(real64) :: ret(1:3) ! kN integer(int32) :: nn, i nn = 0 do i = 1, size(this%ep_domain) nn = nn + this%ep_domain(i)%femdomain%nn() end do t = this%get_internal_force(dU=zeros(nn*3)) ret(1) = sum(t(3*NodeList(:) - 2)) ret(2) = sum(t(3*NodeList(:) - 1)) ret(3) = sum(t(3*NodeList(:) - 0)) end function ! ###################################################### ! ###################################################### subroutine rezeroElastoPlasticity(this) class(ElastoPlasticity_), intent(inout) :: this integer(int32) :: i do i = 1, size(this%ep_domain) this%ep_domain(i)%CauchyStress_field = 0.0d0 this%ep_domain(i)%Strain_field = 0.0d0 this%ep_domain(i)%PlasticStrain_field = 0.0d0 this%ep_domain(i)%displacement = 0.0d0 end do end subroutine ! ###################################################### subroutine exportFieldElastoPlasticity(this, name) class(ElastoPlasticity_), intent(in) :: this character(*), intent(in) :: name integeR(int32):: DomainID if (.not. allocated(this%ep_domain)) then return else do DomainID = 1, size(this%ep_domain) call this%ep_domain(DomainID)%exportField(name + "_"+str(DomainID) + "_") end do end if end subroutine ! ###################################################### subroutine importFieldElastoPlasticity(this, name, num_domain) class(ElastoPlasticity_), intent(inout) :: this character(*), intent(in) :: name integer(int32), intent(in) :: num_domain integeR(int32):: DomainID if (.not. allocated(this%ep_domain)) then allocate (this%ep_domain(DomainID)) end if do DomainID = 1, size(this%ep_domain) call this%ep_domain(DomainID)%importField(name + "_"+str(DomainID) + "_") end do end subroutine ! ###################################################### subroutine exportFieldEpDomain(this, name) class(EP_Domain_), intent(in) :: this character(*), intent(in) :: name type(IO_) :: f call to_csv(name + "_CauchyStress_field", this%CauchyStress_field) call to_csv(name + "_displacement", this%displacement) end subroutine ! ###################################################### ! ###################################################### subroutine importFieldEpDomain(this, name) class(EP_Domain_), intent(inout) :: this character(*), intent(in) :: name type(IO_) :: f integer(int32) :: n1, n2, n3 n1 = 6 n2 = this%femdomain%ngp() n3 = this%femdomain%ne() this%CauchyStress_field = from_csv(name + "_CauchyStress_field", n1, n2, n3) n1 = this%femdomain%nd()*this%femdomain%nn() this%displacement = from_csv(name + "_displacement.csv", n1) end subroutine ! ####################################################### ! ####################################################### function to_EP_Model_ElastoPlastClass(ElasticPotential,YieldFunction,PlasticPotential,StressRatio) result(ret) procedure(E_PotentialFunction) :: ElasticPotential procedure(P_PotentialFunction) :: YieldFunction procedure(P_PotentialFunction) :: PlasticPotential procedure(StressRatioFunction),optional :: StressRatio type(EP_Model_) :: ret ret%ElasticPotential => ElasticPotential ret%YieldFunction => YieldFunction ret%PlasticPotential => PlasticPotential if (present(StressRatio) )then ! if StressRatio is presented, ! then, construct matrices based on ! the hypo-elasto-plasticity ret%StressRatio => StressRatio endif end function ! ####################################################### ! ####################################################### function JaumannStressRatio(sigma, d_sigma, l) result(ret) real(real64), intent(in) :: sigma(:, :), d_sigma(:, :), l(:, :) real(real64), allocatable :: w(:,:) real(real64),allocatable :: ret(:,:) w = 0.50d0*(l - transpose(l)) ret = d_sigma + matmul(sigma,w) - matmul(w,sigma) end function ! ####################################################### ! ####################################################### function TruesdellStressRatio(sigma, d_sigma, l) result(ret) real(real64), intent(in) :: sigma(:, :), d_sigma(:, :), l(:, :) real(real64),allocatable :: ret(:,:) ret = d_sigma - matmul(l,sigma) - matmul(sigma,l) + trace(l)*sigma end function ! ####################################################### ! ####################################################### function OldroydStressRatio(sigma, d_sigma, l) result(ret) real(real64), intent(in) :: sigma(:, :), d_sigma(:, :), l(:, :) real(real64),allocatable :: ret(:,:) ret = d_sigma + matmul(l,sigma) + matmul(sigma,transpose(l)) end function ! ####################################################### ! ####################################################### ! [(potential function),(Stress Ratio)] -> (element-wise Coefficient Matrix) !-------------------------------------------------------- function StiffnessMatrix_EP_model(EP_Model,ElasticParams,PlasticParams,ElasticStrain,nDim) result(ret) class(EP_Model_),intent(in) :: EP_Model real(real64),intent(in) :: ElasticParams(:), PlasticParams(:),ElasticStrain(:,:) real(real64),allocatable :: ret(:,:),stiffness_tensor(:,:,:,:) integer(int32),intent(in) :: nDim integer(int32),allocatable :: stress_matrix_order(:,:) integer(int32) :: i,j,k,l,idx,n,m if(nDim == 1 )then allocate(stress_matrix_order(1,2)) stress_matrix_order(1,1:2) = 1 else allocate(stress_matrix_order(nDim*(nDim-1)/2 + nDim,2)) idx = 0 do i = 1, nDim idx = idx + 1 stress_matrix_order(idx,1:2) = i enddo do i = 1, nDim do j=i+1, nDim idx = idx + 1 stress_matrix_order(idx,1) = i stress_matrix_order(idx,2) = j enddo enddo if(nDim==3)then stress_matrix_order(5,1:2) = [2,3] stress_matrix_order(6,1:2) = [3,1] endif endif ! とりあえず，超弾性のみ実装 ! C_{ijkl} stiffness_tensor = d2_depsilon2(EP_Model%ElasticPotential,ElasticStrain,ElasticParams) ! ret ret = zeros(size(stress_matrix_order,1),size(stress_matrix_order,1)) do m=1,size(stress_matrix_order,1) do n=1,size(stress_matrix_order,1) i = stress_matrix_order(m,1) j = stress_matrix_order(m,2) k = stress_matrix_order(n,1) l = stress_matrix_order(n,2) ret(m,n) = stiffness_tensor(i,j,k,l) enddo enddo end function ! ####################################################### !function SmallStrainTensor_EP_model( ! ####################################################### function imaginaryTimestepDerivative_c(this_func,x,params) result(ret) procedure(ScalarFunction) :: this_func real(real64),intent(in) :: x,params(:) type(Math_) :: math real(real64) :: ret,epsilon epsilon = params(1) ret = aimag( (this_func((x+math%i*epsilon),params)))/epsilon end function ! ####################################################### ! ####################################################### function imaginaryTimestep2ndDerivative(this_func,x,params) result(ret) procedure(ScalarFunction) :: this_func real(real64),intent(in) :: x real(real64),intent(in) :: params(:) type(Math_) :: math real(real64) :: ret,epsilon epsilon = params(1) ret = 2.0d0*real(this_func(dcmplx(x),params))/epsilon/epsilon - & 2.0d0*real( this_func((x+math%i*epsilon),params)/epsilon/epsilon ) end function ! ####################################################### ! ####################################################### function imaginaryTimestepDerivative(this_func,x,params) result(ret) procedure(ScalarFunction) :: this_func real(real64),intent(in) :: x,params(:) type(Math_) :: math real(real64) :: ret,epsilon epsilon = params(1) ret = aimag( this_func((x+math%i*epsilon),params)/epsilon ) end function ! ####################################################### ! ####################################################### function ForwardDifferenceDerivative(this_func,x,params) result(ret) procedure(ScalarFunction) :: this_func real(real64),intent(in) :: x,params(:) complex(real64) :: y1,y2 type(Math_) :: math real(real64) :: ret,epsilon epsilon = params(1) y1 = x+epsilon y2 = x ret = real( (this_func(y1,params)-this_func(y2,params))/epsilon ) end function ! ####################################################### ! ####################################################### function CentralDifferenceDerivative(this_func,x,params) result(ret) procedure(ScalarFunction) :: this_func real(real64),intent(in) :: x,params(:) complex(real64) :: y1,y2 type(Math_) :: math real(real64) :: ret,epsilon epsilon = params(1) y1 = x+epsilon/2.0d0 y2 = x-epsilon/2.0d0 ret = real( (this_func(y1,params)-this_func(y2,params))/epsilon ) end function ! ####################################################### ! ################################################## function is_elastic(val) result(ret) real(real64), intent(in) :: val logical :: ret if (val < 0.0d0) then ret = .true. else ret = .false. end if end function ! ################################################## end module ElastoPlasticityClass