!======================================================================! ! ! ! Software Name : FrontCOMP_cure Ver. 3.1 ! ! ! ! Module Name : m_static_NLGEOM_elpl ! ! Category : Static Analysis ! ! ! ! Developed based on "FrontSTR" of RSS21 project ! ! ! ! Written by Tomotaka Ogasawara, 2006/06/01 ! ! Tomotaka Ogasawara, 2013/03/26 ! ! ! ! Contact address : IIS,The University of Tokyo, CISS ! ! ! ! "Composite Material Strength & Reliability Evaluation Simulator" ! ! ! !======================================================================! !C================================================================C !C-- module m_static_NLGEOM_elpl !C================================================================C module m_static_NLGEOM_elpl contains !C*** !C*** subroutines for NONLINEAR SOLID MECHANICS !C*** !C================================================================C !C-- subroutine fstrNLGEOM_output_monit !C================================================================C subroutine fstrNLGEOM_output_monit(hecMESH, fstrPARAM, fstrSOLID, & elpl_IW4, elpl_IW5, i_step, DT0, & my_rank_monit_1, my_rank_monit_2) ! use m_fstr ! implicit none !C !C-- global variable !C type ( hecmwST_local_mesh ) :: hecMESH type ( fstr_solid ) :: fstrSOLID type ( fstr_param ) :: fstrPARAM !C !C-- local variable !C integer(kind=kint) :: elpl_IW4, elpl_IW5 integer(kind=kint) :: ii,jj integer(kind=kint) :: nb_intg, nb_node, nb_elem integer(kind=kint) :: i_step integer(kind=kint) :: my_rank_monit_1, my_rank_monit_2 real(kind=kreal) :: DT0 real(kind=kreal) :: s11, s22, s33, s12, s23, s31, sbar real(kind=kreal) :: sd11, sd22, sd33, sig_y real(kind=kreal) :: e11, e22, e33, e12, e23, e31, epbar !! nb_intg = fstrSOLID%intg_monit_1 nb_node = fstrSOLID%node_monit_1 nb_elem = fstrSOLID%elem_monit_1 if( hecMESH%my_rank .eq. my_rank_monit_1) then jj=fstrPARAM%global_local_id(1,nb_node) ii=fstrPARAM%global_local_id(2,nb_node) !! !C-- displacement !! if( fstrSOLID%iout_list(1) .eq. 1 ) then write(elpl_IW4,'(i10,1pe13.4e3,i10,1p6e13.4e3)') i_step & ,DT0*i_step & ,jj ,fstrSOLID%TOTAL_DISP( hecMESH%n_dof*(ii-1)+1 : hecMESH%n_dof*ii ) end if end if ! if( hecMESH%my_rank .eq. my_rank_monit_2) then s11 = fstrSOLID%STRESS3D(1,nb_intg,nb_intg,nb_intg,nb_elem) s22 = fstrSOLID%STRESS3D(2,nb_intg,nb_intg,nb_intg,nb_elem) s33 = fstrSOLID%STRESS3D(3,nb_intg,nb_intg,nb_intg,nb_elem) s12 = fstrSOLID%STRESS3D(4,nb_intg,nb_intg,nb_intg,nb_elem) s23 = fstrSOLID%STRESS3D(5,nb_intg,nb_intg,nb_intg,nb_elem) s31 = fstrSOLID%STRESS3D(6,nb_intg,nb_intg,nb_intg,nb_elem) sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 sbar = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)) e11 = fstrSOLID%STRAIN3D_P(1,nb_intg,nb_intg,nb_intg,nb_elem) e22 = fstrSOLID%STRAIN3D_P(2,nb_intg,nb_intg,nb_intg,nb_elem) e33 = fstrSOLID%STRAIN3D_P(3,nb_intg,nb_intg,nb_intg,nb_elem) e12 = fstrSOLID%STRAIN3D_P(4,nb_intg,nb_intg,nb_intg,nb_elem) e23 = fstrSOLID%STRAIN3D_P(5,nb_intg,nb_intg,nb_intg,nb_elem) e31 = fstrSOLID%STRAIN3D_P(6,nb_intg,nb_intg,nb_intg,nb_elem) epbar = fstrSOLID%EP_BAR(nb_intg,nb_intg,nb_intg,nb_elem) !! !C-- stress , strain ,epbar !! if( fstrSOLID%iout_list(6) .eq. 1 ) then write(elpl_IW5,'(i10,1pe13.4e3,i10,1p14e13.4e3)') i_step & ,DT0*i_step ,hecMESH%global_elem_ID(nb_elem) & ! ,s11 ,sbar, epbar ,s11, s22, s33, s12, s23, s31, sbar ,e11, e22, e33, e12, e23, e31, epbar end if end if end subroutine fstrNLGEOM_output_monit !C================================================================C !C-- subroutine fstrNLGEOM_output_monit_2d !C================================================================C subroutine fstrNLGEOM_output_monit_2d(hecMESH, fstrPARAM, fstrSOLID, & elpl_IW4, elpl_IW5, i_step, DT0, & my_rank_monit_1, my_rank_monit_2) ! use m_fstr ! implicit none !C !C-- global variable !C type ( hecmwST_local_mesh ) :: hecMESH type ( fstr_solid ) :: fstrSOLID type ( fstr_param ) :: fstrPARAM !C !C-- local variable !C integer(kind=kint) :: elpl_IW4, elpl_IW5 integer(kind=kint) :: ii,jj integer(kind=kint) :: nb_intg, nb_node, nb_elem integer(kind=kint) :: i_step integer(kind=kint) :: my_rank_monit_1, my_rank_monit_2 real(kind=kreal) :: DT0 real(kind=kreal) :: s11, s22, s33, s12, s23, s31, sbar real(kind=kreal) :: sd11, sd22, sd33, sig_y real(kind=kreal) :: epbar !! nb_intg = fstrSOLID%intg_monit_1 nb_node = fstrSOLID%node_monit_1 nb_elem = fstrSOLID%elem_monit_1 if( hecMESH%my_rank .eq. my_rank_monit_1) then jj=fstrPARAM%global_local_id(1,nb_node) ii=fstrPARAM%global_local_id(2,nb_node) !! !C-- displacement !! if( fstrSOLID%iout_list(1) .eq. 1 ) then write(elpl_IW4,'(i10,1pe13.4e3,i10,1p6e13.4e3)') i_step & ,DT0*i_step & ,jj ,fstrSOLID%TOTAL_DISP( hecMESH%n_dof*(ii-1)+1 : hecMESH%n_dof*ii ) end if end if ! if( hecMESH%my_rank .eq. my_rank_monit_2) then s11 = fstrSOLID%STRESS2D(1,nb_intg,nb_intg,nb_elem) s22 = fstrSOLID%STRESS2D(2,nb_intg,nb_intg,nb_elem) s33 = fstrSOLID%STRESS2D(4,nb_intg,nb_intg,nb_elem) s12 = fstrSOLID%STRESS2D(3,nb_intg,nb_intg,nb_elem) s23 = 0.0 s31 = 0.0 sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 sbar = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)) epbar = fstrSOLID%EP_BAR(nb_intg,nb_intg,1,nb_elem) !! !C-- stress , strain ,epbar !! if( fstrSOLID%iout_list(6) .eq. 1 ) then write(elpl_IW5,'(i10,1pe13.4e3,i10,1p3e13.4e3)') i_step & ,DT0*i_step & ,hecMESH%global_elem_ID(nb_elem) ,s11 ,sbar, epbar end if end if end subroutine fstrNLGEOM_output_monit_2d !C================================================================C !C-- subroutine fstr_elpl_INTG !C !C-- yield judgement !C-- stress integration !C================================================================C subroutine fstr_elpl_INTG(DMAT,STR,STR_YLD, FLG_ELPL, h_dash, EP_BAR, EE, PP, DE & , iter,BIJ,MAT_OTHO1,MAT_ISO) ! use hecmw use m_fstr ! implicit none !C !C-- I/F variable !C INTEGER(kind=kint) :: iter INTEGER(kind=kint) :: FLG_ELPL REAL(kind=kreal) :: DMAT(6,6),STR(6),STR_YLD,h_dash,EE,PP REAL(kind=kreal) :: DE(6),DS(6) REAL(kind=kreal) :: STR_trial(6) REAL(kind=kreal) :: EP_BAR !C !C-- local variable !C INTEGER(kind=kint) :: i,j,k real(kind=kreal) :: GG, GG2 real(kind=kreal) :: s11, s22, s33, s12, s23, s31, sbar, sbar2 real(kind=kreal) :: ds11, ds22, ds33, ds12, ds23, ds31 real(kind=kreal) :: sd11, sd22, sd33, sig_y real(kind=kreal) :: dsd11, dsd22, dsd33 real(kind=kreal) :: x_ganma, x_A, x_B, x_E real(kind=kreal) :: depbar INTEGER(kind=kint) :: MAT_ISO REAL(kind=kreal) :: BIJ(3,3),MAT_OTHO1(9) !! !!-- check parameters !! if(h_dash .lt. 0.0) then if(myrank .eq. 0) then write(imsg,*) 'stop due to fstrSOLID%h_dash < 0.0' end if call hecmw_abort( hecmw_comm_get_comm()) end if !! GG = EE/(1.0+PP)/2.0 !! !!---- elastic !! if( FLG_ELPL .eq. 0) then call fstr_elpl_DMAT(DMAT,STR,FLG_ELPL, h_dash, EE, PP,BIJ,MAT_OTHO1,MAT_ISO) DS = 0.0 do i = 1,6 do j = 1,6 DS(i)=DS(i)+DMAT(i,j)*DE(j) end do end do STR_trial = STR + DS s11 = STR_trial(1) s22 = STR_trial(2) s33 = STR_trial(3) s12 = STR_trial(4) s23 = STR_trial(5) s31 = STR_trial(6) sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 sbar2 = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)) if( sbar2 .lt. STR_YLD .or. iter .gt. 1) then STR = STR_trial else if( sbar2 .ge. STR_YLD ) then FLG_ELPL = 1 s11 = STR(1) s22 = STR(2) s33 = STR(3) s12 = STR(4) s23 = STR(5) s31 = STR(6) sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 ds11 = DS(1) ds22 = DS(2) ds33 = DS(3) ds12 = DS(4) ds23 = DS(5) ds31 = DS(6) dsd11 = ds11 -(ds11+ds22+ds33)/3.0 dsd22 = ds22 -(ds11+ds22+ds33)/3.0 dsd33 = ds33 -(ds11+ds22+ds33)/3.0 x_A = 0.5*(dsd11**2 +dsd22**2 +dsd33**2 +2.0*ds12**2 +2.0*ds23**2 & +2.0*ds31**2) x_E = 0.5*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2) x_E = x_E -1.0/3.0*STR_YLD**2 x_B = 0.5*(sd11*dsd11 +sd22*dsd22 +sd33*dsd33 +2.0*s12*ds12 +2.0*s23*ds23 & +2.0*s31*ds31) x_ganma = (dsqrt(x_B**2-x_A*x_E)-x_B)/x_A STR = STR + x_ganma*DS call fstr_elpl_DMAT(DMAT,STR,FLG_ELPL, h_dash, EE, PP,BIJ,MAT_OTHO1,MAT_ISO) DS = 0.0 do i = 1,6 do j = 1,6 DS(i)=DS(i)+DMAT(i,j)*(1.0-x_ganma)*DE(j) end do end do STR_trial = STR + DS s11 = STR_trial(1) s22 = STR_trial(2) s33 = STR_trial(3) s12 = STR_trial(4) s23 = STR_trial(5) s31 = STR_trial(6) sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 sbar2 = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)) ! DMAT = 0.0 !-- elastic matrix DMAT(1,1)=EE*(1.0-PP)/(1.0-2.0*PP)/(1.0+PP) DMAT(1,2)=EE*PP/(1.0-2.0*PP)/(1.0+PP) DMAT(1,3)=DMAT(1,2) DMAT(2,1)=DMAT(1,2) DMAT(2,2)=DMAT(1,1) DMAT(2,3)=DMAT(1,2) DMAT(3,1)=DMAT(1,3) DMAT(3,2)=DMAT(2,3) DMAT(3,3)=DMAT(1,1) DMAT(4,4)=EE/(1.0+PP)*0.5 DMAT(5,5)=EE/(1.0+PP)*0.5 DMAT(6,6)=EE/(1.0+PP)*0.5 ! DS = 0.0 do i = 1,6 do j = 1,6 DS(i)=DS(i)+DMAT(i,j)*(1.0-x_ganma)*DE(j) end do end do ds11 = DS(1) ds22 = DS(2) ds33 = DS(3) ds12 = DS(4) ds23 = DS(5) ds31 = DS(6) depbar = sd11*ds11 +sd22*ds22 +sd33*ds33 +2.0*s12*ds12 +2.0*s23*ds23 & +2.0*s31*ds31 depbar = depbar/(2.0*GG+2.0/3.0*h_dash)/sbar2 EP_BAR = EP_BAR + depbar sig_y = STR_YLD + h_dash*depbar STR = sig_y/sbar2 * STR_trial STR_YLD = sig_y end if !! !!---- plastic !! else if( FLG_ELPL .eq. 1) then call fstr_elpl_DMAT(DMAT,STR,FLG_ELPL, h_dash, EE, PP,BIJ,MAT_OTHO1,MAT_ISO) DS = 0.0 do i = 1,6 do j = 1,6 DS(i)=DS(i)+DMAT(i,j)*DE(j) end do end do STR_trial = STR + DS s11 = STR_trial(1) s22 = STR_trial(2) s33 = STR_trial(3) s12 = STR_trial(4) s23 = STR_trial(5) s31 = STR_trial(6) sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 sbar2 = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)) ! if( sbar2 .ge. STR_YLD .or. iter .gt. 1) then ! DMAT = 0.0 !-- elastic matrix DMAT(1,1)=EE*(1.0-PP)/(1.0-2.0*PP)/(1.0+PP) DMAT(1,2)=EE*PP/(1.0-2.0*PP)/(1.0+PP) DMAT(1,3)=DMAT(1,2) DMAT(2,1)=DMAT(1,2) DMAT(2,2)=DMAT(1,1) DMAT(2,3)=DMAT(1,2) DMAT(3,1)=DMAT(1,3) DMAT(3,2)=DMAT(2,3) DMAT(3,3)=DMAT(1,1) DMAT(4,4)=EE/(1.0+PP)*0.5 DMAT(5,5)=EE/(1.0+PP)*0.5 DMAT(6,6)=EE/(1.0+PP)*0.5 ! DS = 0.0 do i = 1,6 do j = 1,6 DS(i)=DS(i)+DMAT(i,j)*DE(j) end do end do ds11 = DS(1) ds22 = DS(2) ds33 = DS(3) ds12 = DS(4) ds23 = DS(5) ds31 = DS(6) depbar = sd11*ds11 +sd22*ds22 +sd33*ds33 +2.0*s12*ds12 +2.0*s23*ds23 & +2.0*s31*ds31 depbar = depbar/(2.0*GG+2.0/3.0*h_dash)/sbar2 EP_BAR = EP_BAR + depbar sig_y = STR_YLD + h_dash*depbar STR = sig_y/sbar2 * STR_trial STR_YLD = sig_y else if( sbar2 .lt. STR_YLD ) then FLG_ELPL = 0 call fstr_elpl_DMAT(DMAT,STR,FLG_ELPL, h_dash, EE, PP,BIJ,MAT_OTHO1,MAT_ISO) DS = 0.0 do i = 1,6 do j = 1,6 DS(i)=DS(i)+DMAT(i,j)*DE(j) end do end do STR = STR + DS end if end if ! 100 continue ! end subroutine fstr_elpl_INTG !C================================================================C !C-- subroutine fstr_elpl_INTG_2d !C !C-- yield judgement !C-- stress integration !C================================================================C subroutine fstr_elpl_INTG_2d(DMAT,STR,STR_YLD, FLG_ELPL, h_dash, EP_BAR, EE, PP, DE & ,iter, ISET) ! use hecmw use m_fstr ! implicit none !C !C-- I/F variable !C INTEGER(kind=kint) :: iter INTEGER(kind=kint) :: FLG_ELPL, ISET REAL(kind=kreal) :: DMAT(4,4),STR(4),STR_YLD,h_dash,EE,PP REAL(kind=kreal) :: DE(4),DS(4) REAL(kind=kreal) :: STR_trial(4) REAL(kind=kreal) :: EP_BAR !C !C-- local variable !C INTEGER(kind=kint) :: i,j,k real(kind=kreal) :: GG, GG2, COEF1, COEF2 real(kind=kreal) :: s11, s22, s33, s12, s23, s31, sbar, sbar2 real(kind=kreal) :: ds11, ds22, ds33, ds12, ds23, ds31 real(kind=kreal) :: sd11, sd22, sd33, sig_y real(kind=kreal) :: dsd11, dsd22, dsd33 real(kind=kreal) :: x_ganma, x_A, x_B, x_E real(kind=kreal) :: depbar !! !!-- check parameters !! if(h_dash .lt. 0.0) then if(myrank .eq. 0) then write(imsg,*) 'stop due to fstrSOLID%h_dash < 0.0' end if call hecmw_abort( hecmw_comm_get_comm()) end if !! GG = EE/(1.0+PP)/2.0 !! !!---- elastic !! if( FLG_ELPL .eq. 0) then call fstr_elpl_DMAT_2d(DMAT,STR,FLG_ELPL, h_dash, EE, PP,ISET) DS = 0.0 do i = 1,4 do j = 1,4 DS(i)=DS(i)+DMAT(i,j)*DE(j) end do end do STR_trial = STR + DS s11 = STR_trial(1) s22 = STR_trial(2) s33 = STR_trial(4) s12 = STR_trial(3) s23 = 0.0 s31 = 0.0 sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 sbar2 = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)) if( sbar2 .lt. STR_YLD .or. iter .gt. 1) then STR = STR_trial else if( sbar2 .ge. STR_YLD ) then FLG_ELPL = 1 s11 = STR(1) s22 = STR(2) s33 = STR(4) s12 = STR(3) s23 = 0.0 s31 = 0.0 sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 ds11 = DS(1) ds22 = DS(2) ds33 = DS(4) ds12 = DS(3) ds23 = 0.0 ds31 = 0.0 dsd11 = ds11 -(ds11+ds22+ds33)/3.0 dsd22 = ds22 -(ds11+ds22+ds33)/3.0 dsd33 = ds33 -(ds11+ds22+ds33)/3.0 x_A = 0.5*(dsd11**2 +dsd22**2 +dsd33**2 +2.0*ds12**2 +2.0*ds23**2 & +2.0*ds31**2) x_E = 0.5*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2) x_E = x_E -1.0/3.0*STR_YLD**2 x_B = 0.5*(sd11*dsd11 +sd22*dsd22 +sd33*dsd33 +2.0*s12*ds12 +2.0*s23*ds23 & +2.0*s31*ds31) x_ganma = (dsqrt(x_B**2-x_A*x_E)-x_B)/x_A STR = STR + x_ganma*DS call fstr_elpl_DMAT_2d(DMAT,STR,FLG_ELPL, h_dash, EE, PP,ISET) DS = 0.0 do i = 1,4 do j = 1,4 DS(i)=DS(i)+DMAT(i,j)*(1.0-x_ganma)*DE(j) end do end do STR_trial = STR + DS s11 = STR_trial(1) s22 = STR_trial(2) s33 = STR_trial(4) s12 = STR_trial(3) s23 = 0.0 s31 = 0.0 sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 sbar2 = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)) ! DMAT = 0.0 !-- elastic matrix IF(ISET.EQ.0) THEN !*PLANE STRAIN COEF1=EE/((1.0+PP)*(1.0-2.0*PP)) COEF2=EE/(2.0*(1.0+PP)) DMAT(1,1)=COEF1*(1.0-PP) DMAT(1,2)=COEF1*PP DMAT(1,3)=0.0 DMAT(2,1)=DMAT(1,2) DMAT(2,2)=DMAT(1,1) DMAT(2,3)=0.0 DMAT(3,1)=0.0 DMAT(3,2)=0.0 DMAT(3,3)=COEF2 ELSE IF(ISET.EQ.1) THEN !*PLANE STRESS COEF1=EE/(1.0-PP*PP) COEF2=0.5*(1.0-PP) DMAT(1,1)=COEF1 DMAT(1,2)=COEF1*PP DMAT(1,3)=0.0 DMAT(2,1)=DMAT(1,2) DMAT(2,2)=DMAT(1,1) DMAT(2,3)=0.0 DMAT(3,1)=0.0 DMAT(3,2)=0.0 DMAT(3,3)=COEF1*COEF2 ELSE IF(ISET.EQ.2) THEN !*AXISYMMETRIC COEF1=EE*(1.0-PP)/((1.0+PP)*(1.0-2.0*PP)) COEF2=(1.0-2.0*PP)/(2.0*(1.0-PP)) DMAT(1,1)=COEF1 DMAT(1,2)=COEF1*PP/(1.0-PP) DMAT(1,3)=0.0 DMAT(1,4)=DMAT(1,2) DMAT(2,1)=DMAT(1,2) DMAT(2,2)=DMAT(1,1) DMAT(2,3)=0.0 DMAT(2,4)=DMAT(1,2) DMAT(3,1)=0.0 DMAT(3,2)=0.0 DMAT(3,3)=COEF1*COEF2 DMAT(3,4)=0.0 DMAT(4,1)=DMAT(1,4) DMAT(4,2)=DMAT(2,4) DMAT(4,3)=0.0 DMAT(4,4)=DMAT(1,1) END IF ! DS = 0.0 do i = 1,4 do j = 1,4 DS(i)=DS(i)+DMAT(i,j)*(1.0-x_ganma)*DE(j) end do end do ds11 = DS(1) ds22 = DS(2) ds33 = DS(4) ds12 = DS(3) ds23 = 0.0 ds31 = 0.0 depbar = sd11*ds11 +sd22*ds22 +sd33*ds33 +2.0*s12*ds12 +2.0*s23*ds23 & +2.0*s31*ds31 depbar = depbar/(2.0*GG+2.0/3.0*h_dash)/sbar2 EP_BAR = EP_BAR + depbar sig_y = STR_YLD + h_dash*depbar STR = sig_y/sbar2 * STR_trial STR_YLD = sig_y end if !! !!---- plastic !! else if( FLG_ELPL .eq. 1) then call fstr_elpl_DMAT_2d(DMAT,STR,FLG_ELPL, h_dash, EE, PP,ISET) DS = 0.0 do i = 1,4 do j = 1,4 DS(i)=DS(i)+DMAT(i,j)*DE(j) end do end do STR_trial = STR + DS s11 = STR_trial(1) s22 = STR_trial(2) s33 = STR_trial(4) s12 = STR_trial(3) s23 = 0.0 s31 = 0.0 sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 sbar2 = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)) ! if( sbar2 .ge. STR_YLD .or. iter .gt. 1) then ! DMAT = 0.0 !-- elastic matrix IF(ISET.EQ.0) THEN !*PLANE STRAIN COEF1=EE/((1.0+PP)*(1.0-2.0*PP)) COEF2=EE/(2.0*(1.0+PP)) DMAT(1,1)=COEF1*(1.0-PP) DMAT(1,2)=COEF1*PP DMAT(1,3)=0.0 DMAT(2,1)=DMAT(1,2) DMAT(2,2)=DMAT(1,1) DMAT(2,3)=0.0 DMAT(3,1)=0.0 DMAT(3,2)=0.0 DMAT(3,3)=COEF2 ELSE IF(ISET.EQ.1) THEN !*PLANE STRESS COEF1=EE/(1.0-PP*PP) COEF2=0.5*(1.0-PP) DMAT(1,1)=COEF1 DMAT(1,2)=COEF1*PP DMAT(1,3)=0.0 DMAT(2,1)=DMAT(1,2) DMAT(2,2)=DMAT(1,1) DMAT(2,3)=0.0 DMAT(3,1)=0.0 DMAT(3,2)=0.0 DMAT(3,3)=COEF1*COEF2 ELSE IF(ISET.EQ.2) THEN !*AXISYMMETRIC COEF1=EE*(1.0-PP)/((1.0+PP)*(1.0-2.0*PP)) COEF2=(1.0-2.0*PP)/(2.0*(1.0-PP)) DMAT(1,1)=COEF1 DMAT(1,2)=COEF1*PP/(1.0-PP) DMAT(1,3)=0.0 DMAT(1,4)=DMAT(1,2) DMAT(2,1)=DMAT(1,2) DMAT(2,2)=DMAT(1,1) DMAT(2,3)=0.0 DMAT(2,4)=DMAT(1,2) DMAT(3,1)=0.0 DMAT(3,2)=0.0 DMAT(3,3)=COEF1*COEF2 DMAT(3,4)=0.0 DMAT(4,1)=DMAT(1,4) DMAT(4,2)=DMAT(2,4) DMAT(4,3)=0.0 DMAT(4,4)=DMAT(1,1) END IF ! DS = 0.0 do i = 1,4 do j = 1,4 DS(i)=DS(i)+DMAT(i,j)*DE(j) end do end do ds11 = DS(1) ds22 = DS(2) ds33 = DS(4) ds12 = DS(3) ds23 = 0.0 ds31 = 0.0 depbar = sd11*ds11 +sd22*ds22 +sd33*ds33 +2.0*s12*ds12 +2.0*s23*ds23 & +2.0*s31*ds31 depbar = depbar/(2.0*GG+2.0/3.0*h_dash)/sbar2 EP_BAR = EP_BAR + depbar sig_y = STR_YLD + h_dash*depbar STR = sig_y/sbar2 * STR_trial STR_YLD = sig_y else if( sbar2 .lt. STR_YLD ) then FLG_ELPL = 0 call fstr_elpl_DMAT_2d(DMAT,STR,FLG_ELPL, h_dash, EE, PP,ISET) DS = 0.0 do i = 1,4 do j = 1,4 DS(i)=DS(i)+DMAT(i,j)*DE(j) end do end do STR = STR + DS end if end if ! 100 continue ! end subroutine fstr_elpl_INTG_2d !C================================================================C !C-- subroutine fstr_elpl_DMAT !C !C-- set elasto-plastic matrix for stress integration !C================================================================C subroutine fstr_elpl_DMAT(DMAT,STR,FLG_ELPL, h_dash,EE,PP,BIJ,MAT_OTHO1,MAT_ISO) ! use hecmw use m_fstr ! implicit none !C !C-- I/F variable !C INTEGER(kind=kint) :: FLG_ELPL REAL(kind=kreal) :: DMAT(6,6),STR(6),h_dash,EE,PP !C !C-- local variable !C INTEGER(kind=kint) :: i,j,k,l,i2,j2,k2,l2 real(kind=kreal) :: GG, GG2 real(kind=kreal) :: s11, s22, s33, s12, s23, s31, sbar real(kind=kreal) :: sd11, sd22, sd33 INTEGER(kind=kint) :: MAT_ISO REAL(kind=kreal) :: BIJ(3,3),MAT_OTHO1(9) REAL(kind=kreal) :: CIJKL(3,3,3,3),CIJKLB(3,3,3,3) real(kind=kreal) :: E1,E2,E3,G12,G23,G31,nyu12,nyu23,nyu13 real(kind=kreal) :: nyu21,nyu32,nyu31,delta1 !! DMAT = 0.0 do i = 1,3 do j = 1,3 do k = 1,3 do l = 1,3 CIJKL (i,j,k,l) = 0.0 CIJKLB(i,j,k,l) = 0.0 end do end do end do end do !! !-- elastic matrix ! if( MAT_ISO .eq. 0) then DMAT(1,1)=EE*(1.0-PP)/(1.0-2.0*PP)/(1.0+PP) DMAT(1,2)=EE*PP/(1.0-2.0*PP)/(1.0+PP) DMAT(1,3)=DMAT(1,2) DMAT(2,1)=DMAT(1,2) DMAT(2,2)=DMAT(1,1) DMAT(2,3)=DMAT(1,2) DMAT(3,1)=DMAT(1,3) DMAT(3,2)=DMAT(2,3) DMAT(3,3)=DMAT(1,1) DMAT(4,4)=EE/(1.0+PP)*0.5 DMAT(5,5)=EE/(1.0+PP)*0.5 DMAT(6,6)=EE/(1.0+PP)*0.5 ! else if( MAT_ISO .eq. 1 .or. MAT_ISO .eq. -1) then else if( MAT_ISO .eq. 1) then E1 = MAT_OTHO1(1) E2 = MAT_OTHO1(2) E3 = MAT_OTHO1(3) G12 = MAT_OTHO1(4) G23 = MAT_OTHO1(5) G31 = MAT_OTHO1(6) nyu12 = MAT_OTHO1(7) nyu23 = MAT_OTHO1(8) nyu13 = MAT_OTHO1(9) nyu21 = E2/E1*nyu12 nyu32 = E3/E2*nyu23 nyu31 = E3/E1*nyu13 delta1 = (1.0 -nyu12*nyu21 -nyu23*nyu32 -nyu31*nyu13 -2.0*nyu12*nyu23*nyu31)/(E1*E2*E3) CIJKLB(1,1,1,1) = (1.0/E3 -nyu23*nyu23/E2)/(E2*delta1) CIJKLB(2,2,2,2) = (1.0/E3 -nyu13*nyu13/E1)/(E1*delta1) CIJKLB(3,3,3,3) = (1.0/E2 -nyu12*nyu12/E1)/(E1*delta1) CIJKLB(1,1,2,2) = (nyu12/E3 +nyu13*nyu23/E2)/(E1*delta1) CIJKLB(2,2,3,3) = (nyu23/E2 +nyu12*nyu13/E1)/(E1*delta1) CIJKLB(3,3,1,1) = (nyu13 +nyu12*nyu23)/(E1*E2*delta1) CIJKLB(1,2,1,2) = G12 CIJKLB(2,3,2,3) = G23 CIJKLB(3,1,3,1) = G31 ! CIJKLB(2,2,1,1) = CIJKLB(1,1,2,2) CIJKLB(3,3,2,2) = CIJKLB(2,2,3,3) CIJKLB(1,1,3,3) = CIJKLB(3,3,1,1) CIJKLB(2,1,1,2) = CIJKLB(1,2,1,2) CIJKLB(2,1,2,1) = CIJKLB(1,2,1,2) CIJKLB(1,2,2,1) = CIJKLB(1,2,1,2) CIJKLB(3,2,2,3) = CIJKLB(2,3,2,3) CIJKLB(3,2,3,2) = CIJKLB(2,3,2,3) CIJKLB(2,3,3,2) = CIJKLB(2,3,2,3) CIJKLB(1,3,3,1) = CIJKLB(3,1,3,1) CIJKLB(1,3,1,3) = CIJKLB(3,1,3,1) CIJKLB(3,1,1,3) = CIJKLB(3,1,3,1) ! do i = 1,3 do j = 1,3 do k = 1,3 do l = 1,3 do i2 = 1,3 do j2 = 1,3 do k2 = 1,3 do l2 = 1,3 CIJKL(i,j,k,l) = CIJKL(i,j,k,l) + BIJ(i,i2)*BIJ(j,j2)*BIJ(k,k2)*BIJ(l,l2)*CIJKLB(i2,j2,k2,l2) end do end do end do end do end do end do end do end do ! DMAT(1,1)=CIJKL(1,1,1,1) DMAT(1,2)=CIJKL(1,1,2,2) DMAT(1,3)=CIJKL(1,1,3,3) DMAT(2,1)=CIJKL(2,2,1,1) DMAT(2,2)=CIJKL(2,2,2,2) DMAT(2,3)=CIJKL(2,2,3,3) DMAT(3,1)=CIJKL(3,3,1,1) DMAT(3,2)=CIJKL(3,3,2,2) DMAT(3,3)=CIJKL(3,3,3,3) DMAT(4,4)=CIJKL(1,2,1,2) DMAT(5,5)=CIJKL(2,3,2,3) DMAT(6,6)=CIJKL(3,1,3,1) else if( MAT_ISO .eq. 2) then end if !! !-- elasto-plastic matrix ! if( FLG_ELPL .eq. 1) then GG = EE/(1.0+PP)/2.0 ! s11 = STR(1) s22 = STR(2) s33 = STR(3) s12 = STR(4) s23 = STR(5) s31 = STR(6) sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 sbar = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)) GG2 = 9.0*GG**2/(3.0*GG+h_dash)/sbar**2 !! DMAT(1,1) = DMAT(1,1) - GG2*sd11*sd11 DMAT(1,2) = DMAT(1,2) - GG2*sd11*sd22 DMAT(1,3) = DMAT(1,3) - GG2*sd11*sd33 DMAT(1,4) = DMAT(1,4) - GG2*sd11*s12 DMAT(1,5) = DMAT(1,5) - GG2*sd11*s23 DMAT(1,6) = DMAT(1,6) - GG2*sd11*s31 DMAT(2,1) = DMAT(2,1) - GG2*sd22*sd11 DMAT(2,2) = DMAT(2,2) - GG2*sd22*sd22 DMAT(2,3) = DMAT(2,3) - GG2*sd22*sd33 DMAT(2,4) = DMAT(2,4) - GG2*sd22*s12 DMAT(2,5) = DMAT(2,5) - GG2*sd22*s23 DMAT(2,6) = DMAT(2,6) - GG2*sd22*s31 DMAT(3,1) = DMAT(3,1) - GG2*sd33*sd11 DMAT(3,2) = DMAT(3,2) - GG2*sd33*sd22 DMAT(3,3) = DMAT(3,3) - GG2*sd33*sd33 DMAT(3,4) = DMAT(3,4) - GG2*sd33*s12 DMAT(3,5) = DMAT(3,5) - GG2*sd33*s23 DMAT(3,6) = DMAT(3,6) - GG2*sd33*s31 DMAT(4,1) = DMAT(4,1) - GG2*s12 *sd11 DMAT(4,2) = DMAT(4,2) - GG2*s12 *sd22 DMAT(4,3) = DMAT(4,3) - GG2*s12 *sd33 DMAT(4,4) = DMAT(4,4) - GG2*s12 *s12 DMAT(4,5) = DMAT(4,5) - GG2*s12 *s23 DMAT(4,6) = DMAT(4,6) - GG2*s12 *s31 DMAT(5,1) = DMAT(5,1) - GG2*s23 *sd11 DMAT(5,2) = DMAT(5,2) - GG2*s23 *sd22 DMAT(5,3) = DMAT(5,3) - GG2*s23 *sd33 DMAT(5,4) = DMAT(5,4) - GG2*s23 *s12 DMAT(5,5) = DMAT(5,5) - GG2*s23 *s23 DMAT(5,6) = DMAT(5,6) - GG2*s23 *s31 DMAT(6,1) = DMAT(6,1) - GG2*s31 *sd11 DMAT(6,2) = DMAT(6,2) - GG2*s31 *sd22 DMAT(6,3) = DMAT(6,3) - GG2*s31 *sd33 DMAT(6,4) = DMAT(6,4) - GG2*s31 *s12 DMAT(6,5) = DMAT(6,5) - GG2*s31 *s23 DMAT(6,6) = DMAT(6,6) - GG2*s31 *s31 ! end if ! end subroutine fstr_elpl_DMAT !C================================================================C !C-- subroutine fstr_elpl_DMAT_2d !C !C-- set elasto-plastic matrix for stress integration !C================================================================C subroutine fstr_elpl_DMAT_2d(DMAT,STR,FLG_ELPL, h_dash, EE, PP, ISET) ! use hecmw ! implicit none !C !C-- I/F variable !C INTEGER(kind=kint) :: FLG_ELPL, ISET REAL(kind=kreal) :: DMAT(4,4),STR(4),h_dash,EE,PP !C !C-- local variable !C real(kind=kreal) :: GG, GG2, COEF1, COEF2 real(kind=kreal) :: s11, s22, s33, s12, s23, s31, sbar real(kind=kreal) :: sd11, sd22, sd33 real(kind=kreal) :: ss11, ss22, ss12 !! DMAT = 0.0 !! !-- elastic matrix ! IF(ISET.EQ.0) THEN !*PLANE STRAIN COEF1=EE/((1.0+PP)*(1.0-2.0*PP)) COEF2=EE/(2.0*(1.0+PP)) DMAT(1,1)=COEF1*(1.0-PP) DMAT(1,2)=COEF1*PP DMAT(1,3)=0.0 DMAT(2,1)=DMAT(1,2) DMAT(2,2)=DMAT(1,1) DMAT(2,3)=0.0 DMAT(3,1)=0.0 DMAT(3,2)=0.0 DMAT(3,3)=COEF2 ELSE IF(ISET.EQ.1) THEN !*PLANE STRESS COEF1=EE/(1.0-PP*PP) COEF2=0.5*(1.0-PP) DMAT(1,1)=COEF1 DMAT(1,2)=COEF1*PP DMAT(1,3)=0.0 DMAT(2,1)=DMAT(1,2) DMAT(2,2)=DMAT(1,1) DMAT(2,3)=0.0 DMAT(3,1)=0.0 DMAT(3,2)=0.0 DMAT(3,3)=COEF1*COEF2 ELSE IF(ISET.EQ.2) THEN !*AXISYMMETRIC COEF1=EE*(1.0-PP)/((1.0+PP)*(1.0-2.0*PP)) COEF2=(1.0-2.0*PP)/(2.0*(1.0-PP)) DMAT(1,1)=COEF1 DMAT(1,2)=COEF1*PP/(1.0-PP) DMAT(1,3)=0.0 DMAT(1,4)=DMAT(1,2) DMAT(2,1)=DMAT(1,2) DMAT(2,2)=DMAT(1,1) DMAT(2,3)=0.0 DMAT(2,4)=DMAT(1,2) DMAT(3,1)=0.0 DMAT(3,2)=0.0 DMAT(3,3)=COEF1*COEF2 DMAT(3,4)=0.0 DMAT(4,1)=DMAT(1,4) DMAT(4,2)=DMAT(2,4) DMAT(4,3)=0.0 DMAT(4,4)=DMAT(1,1) END IF !! !-- elasto-plastic matrix ! if( FLG_ELPL .eq. 1) then ! IF(ISET.EQ.0) THEN !*PLANE STRAIN GG = EE/(1.0+PP)/2.0 ! s11 = STR(1) s22 = STR(2) s33 = STR(4) s12 = STR(3) s23 = 0.0 s31 = 0.0 sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 sbar = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)) GG2 = 9.0*GG**2/(3.0*GG+h_dash)/sbar**2 !! DMAT(1,1) = DMAT(1,1) - GG2*sd11*sd11 DMAT(1,2) = DMAT(1,2) - GG2*sd11*sd22 DMAT(1,3) = DMAT(1,3) - GG2*sd11*s12 DMAT(2,1) = DMAT(2,1) - GG2*sd22*sd11 DMAT(2,2) = DMAT(2,2) - GG2*sd22*sd22 DMAT(2,3) = DMAT(2,3) - GG2*sd22*s12 DMAT(3,1) = DMAT(3,1) - GG2*s12 *sd11 DMAT(3,2) = DMAT(3,2) - GG2*s12 *sd22 DMAT(3,3) = DMAT(3,3) - GG2*s12 *s12 ELSE IF(ISET.EQ.1) THEN !*PLANE STRESS GG = EE/(1.0+PP)/2.0 ! s11 = STR(1) s22 = STR(2) s33 = STR(4) s12 = STR(3) s23 = 0.0 s31 = 0.0 sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 sbar = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)) ss11 = DMAT(1,1)*sd11 + DMAT(1,2)*sd22 + 2.0*DMAT(1,3)*s12 ss22 = DMAT(2,1)*sd11 + DMAT(2,2)*sd22 + 2.0*DMAT(2,3)*s12 ss12 = DMAT(3,1)*sd11 + DMAT(3,2)*sd22 + 2.0*DMAT(3,3)*s12 GG2 = 1.0 /((ss11*sd11 + ss22*sd22 + 2.0*ss12*s12) + 4.0/9.0*(sbar**2)*h_dash) !! DMAT(1,1) = DMAT(1,1) - GG2*ss11*ss11 DMAT(1,2) = DMAT(1,2) - GG2*ss11*ss22 DMAT(1,3) = DMAT(1,3) - GG2*ss11*ss12 DMAT(2,1) = DMAT(2,1) - GG2*ss22*ss11 DMAT(2,2) = DMAT(2,2) - GG2*ss22*ss22 DMAT(2,3) = DMAT(2,3) - GG2*ss22*ss12 DMAT(3,1) = DMAT(3,1) - GG2*ss12*ss11 DMAT(3,2) = DMAT(3,2) - GG2*ss12*ss22 DMAT(3,3) = DMAT(3,3) - GG2*ss12*ss12 ELSE IF(ISET.EQ.2) THEN !*AXISYMMETRIC GG = EE/(1.0+PP)/2.0 ! s11 = STR(1) s22 = STR(2) s33 = STR(4) s12 = STR(3) s23 = 0.0 s31 = 0.0 sd11 = s11 -(s11+s22+s33)/3.0 sd22 = s22 -(s11+s22+s33)/3.0 sd33 = s33 -(s11+s22+s33)/3.0 sbar = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)) GG2 = 9.0*GG**2/(3.0*GG+h_dash)/sbar**2 !! DMAT(1,1) = DMAT(1,1) - GG2*sd11*sd11 DMAT(1,2) = DMAT(1,2) - GG2*sd11*sd22 DMAT(1,3) = DMAT(1,3) - GG2*sd11*s12 DMAT(1,4) = DMAT(1,4) - GG2*sd11*sd33 DMAT(2,1) = DMAT(2,1) - GG2*sd22*sd11 DMAT(2,2) = DMAT(2,2) - GG2*sd22*sd22 DMAT(2,3) = DMAT(2,3) - GG2*sd22*s12 DMAT(2,4) = DMAT(2,4) - GG2*sd22*sd33 DMAT(3,1) = DMAT(3,1) - GG2*s12 *sd11 DMAT(3,2) = DMAT(3,2) - GG2*s12 *sd22 DMAT(3,3) = DMAT(3,3) - GG2*s12 *s12 DMAT(3,4) = DMAT(3,4) - GG2*s12 *sd33 DMAT(4,1) = DMAT(4,1) - GG2*sd33*sd11 DMAT(4,2) = DMAT(4,2) - GG2*sd33*sd22 DMAT(4,3) = DMAT(4,3) - GG2*sd33*s12 DMAT(4,4) = DMAT(4,4) - GG2*sd33*sd33 END IF ! end if ! end subroutine fstr_elpl_DMAT_2d !C================================================================C !C-- subroutine table_nlsta !C================================================================C subroutine table_nlsta(hecMESH, fstrSOLID, i_step, ig0, f_t, flag_u) use hecmw use m_fstr implicit none !C !C-- global variable !C type ( hecmwST_local_mesh ) :: hecMESH type ( fstr_solid ) :: fstrSOLID !C !C-- local variable !C integer(kind=kint) :: i, ig0, i_step integer(kind=kint) :: jj_n_amp, jj1, jj2 integer(kind=kint) :: s1, s2, flag, flag_u real(kind=kreal) :: t_1, t_2, t_t, f_1, f_2, f_t !C !C-- zero clear !C i = 0 s1 = 0 s2 = 0 flag = 0 t_1 = 0.0 t_2 = 0.0 t_t = 0.0 f_1 = 0.0 f_2 = 0.0 f_t = 0.0 !! ! fstrSOLID%NLSTATIC_ngrp_amp = 1 ! jj_n_amp = -1 !! jj_n_amp = fstrSOLID%NLSTATIC_ngrp_amp if( jj_n_amp .eq. 0 ) then t_t = DT*i_step f_t = t_t else jj1 = hecMESH%amp%amp_index(jj_n_amp - 1) jj2 = hecMESH%amp%amp_index(jj_n_amp) jj1 = jj1 + 2 t_t = DT*i_step ! if(jj2 .eq. 0) then ! f_t = 1.0 if(t_t .gt. hecMESH%amp%amp_table(jj2)) then f_t = hecMESH%amp%amp_val(jj2) else if(t_t .le. hecMESH%amp%amp_table(jj2)) then do i = jj1, jj2 if(t_t .le. hecMESH%amp%amp_table(i)) then s2 = i s1 = i - 1 flag = 1 end if if( flag .eq. 1 ) exit end do t_2 = hecMESH%amp%amp_table(s2) t_1 = hecMESH%amp%amp_table(s1) f_2 = hecMESH%amp%amp_val(s2) f_1 = hecMESH%amp%amp_val(s1) if( t_2-t_1 .lt. 1.0e-20) then if(myrank .eq. 0) then write(imsg,*) 'stop due to t_2-t_1 <= 0' end if call hecmw_abort( hecmw_comm_get_comm()) end if f_t = ((t_2*f_1 - t_1*f_2) + (f_2 - f_1)*t_t) / (t_2 - t_1) end if end if end subroutine table_nlsta !C================================================================C !C-- subroutine fstrNLGEOM_get_prop !C================================================================C !--- get material properties SUBROUTINE fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick) use m_fstr IMPLICIT REAL(kind=kreal) (A-H,O-Z) TYPE (hecmwST_local_mesh) :: hecMESH !*EHM June 24 04 !Get thickness n_item = hecMESH%section%sect_R_index(isect)-hecMESH%section%sect_R_index(isect-1) !Print *,'n_item for thickness: ',n_item !PAUSE ihead = hecMESH%section%sect_R_index(isect-1) !Print *,'ihead for thickness: ',ihead !PAUSE !do i = 1,n_item thick = hecMESH%section%sect_R_item(ihead+1) IF(thick.LE.0.0) STOP "Zero thickness <= 0 is illegal" !Print *,'cval:',cval !PAUSE !end do !C** material ID mid=hecMESH%section%sect_mat_ID_item(isect) !C** Number of Item n_item=hecMESH%material%mat_ITEM_index(mid)-hecMESH%material%mat_ITEM_index(mid-1) !C** Head possition ihead=hecMESH%material%mat_ITEM_index(mid-1) !C Get ITEM of Meterial (Young's Modulus & Possion's Ratio pp=0.0 if( n_item .lt. 1 ) then write(IMSG,*) 'n_item=',n_item write(IMSG,*) '###Error 1' stop endif if ( n_item.ge.1 ) then !C** Number of Sub Item n_subitem = hecMESH%material%mat_subITEM_index(ihead+1) & - hecMESH%material%mat_subITEM_index(ihead) !C** Head possition mpos=hecMESH%material%mat_subITEM_index(ihead) !C Get SUBITEM of Meterial if( n_subitem .lt. 1 ) then write(IMSG,*) '###Error 2' stop endif if ( n_subitem.ge.1 ) then ee=hecMESH%material%mat_val(mpos+1) endif if ( n_subitem.ge.2 ) then pp=hecMESH%material%mat_val(mpos+2) endif endif !C Get ITEM of Meterial (Density) rho=0.0 if ( n_item.ge.2 ) then !C** Number of Sub Item n_subitem=hecMESH%material%mat_subITEM_index(ihead+2) & -hecMESH%material%mat_subITEM_index(ihead+1) !C** Head possition mpos=hecMESH%material%mat_subITEM_index(ihead+1) !C Get SUBITEM of Meterial if( n_subitem .lt. 1 ) then write(IMSG,*) '###Error 3' stop endif if ( n_subitem.ge.1 ) then rho=hecMESH%material%mat_val(mpos+1) endif endif !C Get ITEM of Meterial (Thermal Expansion) alpha=0.0 if ( n_item.ge.3 ) then !C** Number of Sub Item n_subitem=hecMESH%material%mat_subITEM_index(ihead+3) & -hecMESH%material%mat_subITEM_index(ihead+2) !C** Head possition mpos=hecMESH%material%mat_subITEM_index(ihead+2) !C Get SUBITEM of Meterial if( n_subitem .lt. 1 ) then write(IMSG,*) '###Error 4' stop endif if ( n_subitem.ge.1 ) then alpha=hecMESH%material%mat_val(mpos+1) endif endif end subroutine fstrNLGEOM_get_prop !C================================================================C !C-- subroutine fstr_elpl_eig_jcb !C !C-- eigen solver for 3 x 3 matrix !C-- jacobi method !C================================================================C subroutine fstr_elpl_eig_jcb(ARRAY_2, ARRAY_PRIS_2,i_step,G_ELM_ID) ! use hecmw use m_fstr ! implicit none !C !C-- I/F variable !C REAL(kind=kreal) :: ARRAY_2(6),ARRAY_PRIS_2(12) !C !C-- local variable !C INTEGER(kind=kint), parameter :: n = 3 INTEGER(kind=kint), parameter :: imax = 100 real(kind=kreal), parameter :: eps_eig = 1.0e-06 INTEGER(kind=kint) :: i,j,k,l,m,i2,j2,k2,l2 INTEGER(kind=kint) :: ic, ic2, i_step, G_ELM_ID INTEGER(kind=kint) :: idx_2(3), idx_2_iv(3) real(kind=kreal) :: a(n,n), x(n,n), amax, akm, q, co, si, ww real(kind=kreal) :: amax2 !! do i = 1,12 ARRAY_PRIS_2(i) = 0.0 end do do i = 1,3 do j = 1,3 x(i,j) = 0.0 end do end do a(1,1) = ARRAY_2(1) a(1,2) = ARRAY_2(4) a(1,3) = ARRAY_2(6) a(2,1) = ARRAY_2(4) a(2,2) = ARRAY_2(2) a(2,3) = ARRAY_2(5) a(3,1) = ARRAY_2(6) a(3,2) = ARRAY_2(5) a(3,3) = ARRAY_2(3) !!! do k = 1, n do m = 1, n if(k .eq. m) then x(k,m) = 1.0 else x(k,m) = 0.0 end if end do end do ! do ic = 1, imax ic2 = 0 amax = 0.0 ! do k = 1, n-1 do m = k+1, n if(m .le. n) then akm =abs(a(k,m)) if(akm .gt. amax) then i = k j = m amax = akm end if end if end do end do ! if(amax .lt. eps_eig) then ic2 = 1 exit end if ! if(a(i,i) .eq. a(j,j)) then q = 0.5 * acos(0.0) else q = 0.5 * atan(2.0*a(j,i)/(a(i,i)-a(j,j))) end if co = cos(q) si = sin(q) ! do k = 1, n ww = a(i,k) a(i,k) = ww*co + a(j,k)*si a(j,k) =-ww*si + a(j,k)*co if((k.ne.i) .and. (k.ne.j)) then a(k,i) = a(i,k) a(k,j) = a(j,k) end if end do a(i,i) = a(i,i)*co + a(i,j)*si a(j,j) =-a(j,i)*si + a(j,j)*co a(i,j) = 0.0 a(j,i) = 0.0 ! do k = 1, n ww = x(k,i) x(k,i) = ww*co + x(k,j)*si x(k,j) =-ww*si + x(k,j)*co end do end do ! if(i_step .ge. -100) then if(ic2 .eq. 0) then write(idbg,'(''not convergence at eigen solver of principal value: '' & ,'' step ='',i10,'' global_elem_ID(ic) ='',i10,'' imax ='',i10 )') & i_step, G_ELM_ID, ic end if end if !!! idx_2(1) = 1 idx_2_iv(1) = 1 amax2 = a(1,1) do i = 2,3 if(a(i,i) .gt. amax2) then idx_2(i) = 1 idx_2_iv(1) = i amax2 = a(i,i) end if end do if(idx_2_iv(1) .eq. 1) then if(a(2,2) .gt. a(3,3)) then idx_2(2) = 2 idx_2(3) = 3 idx_2_iv(2) = 2 idx_2_iv(3) = 3 else if(a(2,2) .le. a(3,3)) then idx_2(2) = 3 idx_2(3) = 2 idx_2_iv(2) = 3 idx_2_iv(3) = 2 end if else if(idx_2_iv(1) .eq. 2) then if(a(3,3) .gt. a(1,1)) then idx_2(3) = 2 idx_2(1) = 3 idx_2_iv(2) = 3 idx_2_iv(3) = 1 else if(a(3,3) .le. a(1,1)) then idx_2(3) = 3 idx_2(1) = 2 idx_2_iv(2) = 1 idx_2_iv(3) = 3 end if else if(idx_2_iv(1) .eq. 3) then if(a(1,1) .gt. a(2,2)) then idx_2(1) = 2 idx_2(2) = 3 idx_2_iv(2) = 1 idx_2_iv(3) = 2 else if(a(1,1) .le. a(2,2)) then idx_2(1) = 3 idx_2(2) = 2 idx_2_iv(2) = 2 idx_2_iv(3) = 1 end if end if do i = 1,3 ARRAY_PRIS_2( i) = a(idx_2_iv(i),idx_2_iv(i)) end do do i = 1,3 ARRAY_PRIS_2(3+i) = x(i,idx_2_iv(1)) ARRAY_PRIS_2(6+i) = x(i,idx_2_iv(2)) ARRAY_PRIS_2(9+i) = x(i,idx_2_iv(3)) end do !! do i = 1,3 ! ARRAY_PRIS_2( i) = a(i,i) ! end do ! do i = 1,3 ! ARRAY_PRIS_2(3+i) = x(i,1) ! ARRAY_PRIS_2(6+i) = x(i,2) ! ARRAY_PRIS_2(9+i) = x(i,3) !! end do ! end subroutine fstr_elpl_eig_jcb !C================================================================C !C-- subroutine fstr_elpl_eig_qr !C !C-- eigen solver for 3 x 3 matrix !C-- jacobi method (2011.06.01) !C-- QR method (not yet included) !C================================================================C subroutine fstr_elpl_eig_qr(ARRAY_2, ARRAY_PRIS_2,i_step,G_ELM_ID) ! use hecmw use m_fstr ! implicit none !C !C-- I/F variable !C REAL(kind=kreal) :: ARRAY_2(6),ARRAY_PRIS_2(12) !C !C-- local variable !C INTEGER(kind=kint), parameter :: n = 3 INTEGER(kind=kint), parameter :: imax = 100 real(kind=kreal), parameter :: eps_eig = 1.0e-06 INTEGER(kind=kint) :: i,j,k,l,m,i2,j2,k2,l2 INTEGER(kind=kint) :: ic, ic2, i_step, G_ELM_ID INTEGER(kind=kint) :: idx_2(3), idx_2_iv(3) real(kind=kreal) :: a(n,n), x(n,n), amax, akm, q, co, si, ww real(kind=kreal) :: amax2 !! do i = 1,12 ARRAY_PRIS_2(i) = 0.0 end do do i = 1,3 do j = 1,3 x(i,j) = 0.0 end do end do a(1,1) = ARRAY_2(1) a(1,2) = ARRAY_2(4) a(1,3) = ARRAY_2(6) a(2,1) = ARRAY_2(4) a(2,2) = ARRAY_2(2) a(2,3) = ARRAY_2(5) a(3,1) = ARRAY_2(6) a(3,2) = ARRAY_2(5) a(3,3) = ARRAY_2(3) !!! do k = 1, n do m = 1, n if(k .eq. m) then x(k,m) = 1.0 else x(k,m) = 0.0 end if end do end do ! do ic = 1, imax ic2 = 0 amax = 0.0 ! do k = 1, n-1 do m = k+1, n if(m .le. n) then akm =abs(a(k,m)) if(akm .gt. amax) then i = k j = m amax = akm end if end if end do end do ! if(amax .lt. eps_eig) then ic2 = 1 exit end if ! if(a(i,i) .eq. a(j,j)) then q = 0.5 * acos(0.0) else q = 0.5 * atan(2.0*a(j,i)/(a(i,i)-a(j,j))) end if co = cos(q) si = sin(q) ! do k = 1, n ww = a(i,k) a(i,k) = ww*co + a(j,k)*si a(j,k) =-ww*si + a(j,k)*co if((k.ne.i) .and. (k.ne.j)) then a(k,i) = a(i,k) a(k,j) = a(j,k) end if end do a(i,i) = a(i,i)*co + a(i,j)*si a(j,j) =-a(j,i)*si + a(j,j)*co a(i,j) = 0.0 a(j,i) = 0.0 ! do k = 1, n ww = x(k,i) x(k,i) = ww*co + x(k,j)*si x(k,j) =-ww*si + x(k,j)*co end do end do ! if(i_step .ge. -100) then if(ic2 .eq. 0) then write(idbg,'(''not convergence at eigen solver of principal value: '' & ,'' step ='',i10,'' global_elem_ID(ic) ='',i10,'' imax ='',i10 )') & i_step, G_ELM_ID, ic end if end if !!! idx_2(1) = 1 idx_2_iv(1) = 1 amax2 = a(1,1) do i = 2,3 if(a(i,i) .gt. amax2) then idx_2(i) = 1 idx_2_iv(1) = i amax2 = a(i,i) end if end do if(idx_2_iv(1) .eq. 1) then if(a(2,2) .gt. a(3,3)) then idx_2(2) = 2 idx_2(3) = 3 idx_2_iv(2) = 2 idx_2_iv(3) = 3 else if(a(2,2) .le. a(3,3)) then idx_2(2) = 3 idx_2(3) = 2 idx_2_iv(2) = 3 idx_2_iv(3) = 2 end if else if(idx_2_iv(1) .eq. 2) then if(a(3,3) .gt. a(1,1)) then idx_2(3) = 2 idx_2(1) = 3 idx_2_iv(2) = 3 idx_2_iv(3) = 1 else if(a(3,3) .le. a(1,1)) then idx_2(3) = 3 idx_2(1) = 2 idx_2_iv(2) = 1 idx_2_iv(3) = 3 end if else if(idx_2_iv(1) .eq. 3) then if(a(1,1) .gt. a(2,2)) then idx_2(1) = 2 idx_2(2) = 3 idx_2_iv(2) = 1 idx_2_iv(3) = 2 else if(a(1,1) .le. a(2,2)) then idx_2(1) = 3 idx_2(2) = 2 idx_2_iv(2) = 2 idx_2_iv(3) = 1 end if end if do i = 1,3 ARRAY_PRIS_2( i) = a(idx_2_iv(i),idx_2_iv(i)) end do do i = 1,3 ARRAY_PRIS_2(3+i) = x(i,idx_2_iv(1)) ARRAY_PRIS_2(6+i) = x(i,idx_2_iv(2)) ARRAY_PRIS_2(9+i) = x(i,idx_2_iv(3)) end do !! do i = 1,3 ! ARRAY_PRIS_2( i) = a(i,i) ! end do ! do i = 1,3 ! ARRAY_PRIS_2(3+i) = x(i,1) ! ARRAY_PRIS_2(6+i) = x(i,2) ! ARRAY_PRIS_2(9+i) = x(i,3) !! end do ! end subroutine fstr_elpl_eig_qr !C================================================================C !C-- subroutine fstrNLGEOM_ReadTemp !C================================================================C subroutine fstrNLGEOM_ReadTemp(hecMESH,hecMAT,fstrSOLID,fstrPARAM & ,heat_sta_IW1,heat_sta_IW2 & ,i_step, DT0) use m_fstr implicit none type (hecmwST_local_mesh) :: hecMESH type (hecmwST_matrix ) :: hecMAT type (fstr_solid ) :: fstrSOLID type (fstr_param ) :: fstrPARAM ! Local variable integer(kind=kint) :: i,j,inod,ielem,g1,g2,g3 integer(kind=kint) :: heat_sta_IW1, heat_sta_IW2 integer(kind=kint) :: i_step, incr real(kind=kreal) :: bb(27) real(kind=kreal) :: DT0, CTIME character(len= 30) :: cdummy,cdummy2(6) character(len=100) :: cdummy3 character(len=4) :: cdummy4(2) character(len=8) :: cdummy5(2) character(len=1) :: cdummy6(2) ! character(len=*) :: cdummy7(6) character :: cdummy7(6) !C !C-- TEMPERATURE !C ! read (heat_sta_IW1,*) cdummy2(1) !! read (heat_sta_IW1 ) cdummy2(1) read (heat_sta_IW1 ) cdummy4(1) ! write(idbg ,*) cdummy4(1) ! read (heat_sta_IW1,*) (cdummy2(j),j=1,6) ! read (heat_sta_IW1 ) (cdummy2(j),j=1,6) read (heat_sta_IW1 ) cdummy5(1),incr,cdummy5(2),CTIME ! write(idbg ,*) cdummy5(1),incr,cdummy5(2),CTIME ! write(idbg ,*) (cdummy2(j),j=1,6) do i= 1, hecMESH%nn_internal inod=fstrPARAM%global_local_id(1,i) ! read (heat_sta_IW1,*) inod,fstrSOLID%TEMPC(i) read (heat_sta_IW1 ) inod,fstrSOLID%TEMPC(i) ! write(idbg ,*) inod,fstrSOLID%TEMPC(i) enddo !C !C-- CURE !C ! read (heat_sta_IW2,*) cdummy2(1) ! read (heat_sta_IW2 ) cdummy2(1) read (heat_sta_IW2 ) cdummy4(1) ! write(idbg ,*) cdummy4(1) ! read (heat_sta_IW2,*) (cdummy2(j),j=1,6) ! read (heat_sta_IW2 ) (cdummy2(j),j=1,6) read (heat_sta_IW2 ) cdummy5(1),incr,cdummy5(2),CTIME ! write(idbg ,*) cdummy5(1),incr,cdummy5(2),CTIME ! write(idbg ,*) (cdummy2(j),j=1,6) do i= 1, hecMESH%nn_internal inod=fstrPARAM%global_local_id(1,i) ! read (heat_sta_IW2,*) inod,fstrSOLID%CURE3D_nod(i) read (heat_sta_IW2 ) inod,fstrSOLID%CURE3D_nod(i) ! write(idbg ,*) inod,fstrSOLID%CURE3D_nod(i) enddo end subroutine fstrNLGEOM_ReadTemp end module m_static_NLGEOM_elpl