!======================================================================! ! ! ! Software Name : FrontCOMP_cure Ver. 3.1 ! ! ! ! Module Name : m_heat_LIB_DFLUX ! ! Category : Library ! ! ! ! Developed based on "FrontSTR" of RSS21 project ! ! ! ! Written by Toshio Nagashima, 2006/06/01 ! ! Yasuji Fukahori, 2006/06/01 ! ! Tomotaka Ogasawara, 2013/03/26 ! ! ! ! Contact address : IIS,The University of Tokyo, CISS ! ! ! ! "Composite Material Strength & Reliability Evaluation Simulator" ! ! ! !======================================================================! module m_heat_LIB_DFLUX contains !C************************************************************************ !C* DFLUX_111(NN,XX,YY,ZZ,ASECT,LTYPE,VAL,VECT) !C* DFLUX_231(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT) !C* DFLUX_232(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT) !C* DFLUX_241(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT) !C* DFLUX_242(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT) !C* DFLUX_341(NN,XX,YY,ZZ,LTYPE,VAL,VECT) !C* heat_DFLUX_342(NN,XX,YY,ZZ,LTYPE,VAL,VECT & !C* ,TT,IMAT,DTIME,k_elm & !C* ,ntab2,temp2,funcA2,funcB2 & !C* ,idx_cure & !C* ,ccc_rl & !C* ,CURE3D_nod,idbg) !C* DFLUX_351(NN,XX,YY,ZZ,LTYPE,VAL,VECT) !C* DFLUX_352(NN,XX,YY,ZZ,LTYPE,VAL,VECT) !C* DFLUX_361(NN,XX,YY,ZZ,LTYPE,VAL,VECT) !C* DFLUX_362(NN,XX,YY,ZZ,LTYPE,VAL,VECT) !C* DFLUX_731(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT) !C* DFLUX_741(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT) !C* NRCV_342CURE(ILOG,XX,YY,ZZ,CURE,ARRAY) !C* cure_DFLUX_342(NN,XX,YY,ZZ,LTYPE,VAL,VECT & !C* ,TT,IMAT,DTIME,k_elm & !C* ,ntab2,temp2,funcA2,funcB2 & !C* ,idx_cure & !C* ,ccc_rl & !C* ,CURE3D_nod,idbg) !C*--------------------------------------------------------------------* SUBROUTINE heat_DFLUX_111(NN,XX,YY,ZZ,ASECT,LTYPE,VAL,VECT) !C*--------------------------------------------------------------------* !C*** !C*** SET 111 DFLUX !C*** !C* LTYPE=0 : BODY FLUX use hecmw IMPLICIT REAL(kind=kreal)(A-H,O-Z) !C* I/F VARIABLES INTEGER(kind=kint) NN, LTYPE REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),VAL,VECT(NN) !C* !C* IF( LTYPE.EQ.0 ) THEN ASECT = 1.0d0 * ASECT ELSE ASECT = 0.0d0 ENDIF DO I=1,NN VECT(I)=0.0 ENDDO !C* DX = XX(2) - XX(1) DY = YY(2) - YY(1) DZ = ZZ(2) - ZZ(1) AL = DSQRT( DX*DX + DY*DY + DZ*DZ ) VV = ASECT * AL ! VECT(1)= -VAL*VV/2.0d0 VECT(2)= -VAL*VV/2.0d0 ! RETURN end subroutine heat_DFLUX_111 !C*--------------------------------------------------------------------* SUBROUTINE heat_DFLUX_231(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT) !C*--------------------------------------------------------------------* !C*** !C*** SET 231 DFLUX !C*** !C* LTYPE=0 : BODY FLUX use hecmw IMPLICIT REAL(kind=kreal)(A-H,O-Z) !C* I/F VARIABLES INTEGER(kind=kint) NN, LTYPE REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),THICK,VAL,VECT(NN) REAL(kind=kreal) XG(2),WGT(2),H(4),HR(4),HS(4) INTEGER(kind=kint) NOD(2) !************************* ! GAUSS INTEGRATION POINT !************************* DATA XG/-0.5773502691896,0.5773502691896/ DATA WGT/1.0,1.0/ !C* IVOL=0 ISUF=0 IF( LTYPE.EQ.0 ) THEN IVOL=1 ELSE IF( LTYPE.GE.1 ) THEN ISUF=1 IF(LTYPE.EQ.1) THEN NOD(1)=1 NOD(2)=2 ELSE IF(LTYPE.EQ.2) THEN NOD(1)=2 NOD(2)=3 ELSE IF(LTYPE.EQ.3) THEN NOD(1)=3 NOD(2)=1 ENDIF ENDIF !** SURFACE LOAD IF( ISUF.EQ.1 ) THEN DO I=1,NN VECT(I)=0.0 ENDDO !** INTEGRATION OVER SURFACE DO LX=1,2 RI=XG(LX) H(1)=0.5*(1.0-RI) H(2)=0.5*(1.0+RI) HR(1)=-0.5 HR(2)= 0.5 GX=0.0 GY=0.0 DO I=1,2 GX=GX+HR(I)*XX(NOD(I)) GY=GY+HR(I)*YY(NOD(I)) ENDDO XSUM=DSQRT(GX*GX+GY*GY)*THICK DO I=1,2 VECT(NOD(I))=VECT(NOD(I))-XSUM*WGT(LX)*H(I)*VAL ENDDO ENDDO ENDIF !** VOLUME LOAD IF( IVOL.EQ.1 ) THEN DO I=1,NN VECT(I)=0.0 ENDDO V1X=XX(2)-XX(1) V1Y=YY(2)-YY(1) V1Z=ZZ(2)-ZZ(1) V2X=XX(3)-XX(1) V2Y=YY(3)-YY(1) V2Z=ZZ(3)-ZZ(1) V3X= V1Y*V2Z-V1Z*V2Y V3Y= V1Z*V2X-V1X*V2Z V3Z= V1X*V2Y-V1Y*V2X AA=0.5*DSQRT( V3X*v3X + V3Y*V3Y + V3Z*V3Z ) VV=AA*THICK VECT(1)= -VAL*VV/3.0 VECT(2)= -VAL*VV/3.0 VECT(3)= -VAL*VV/3.0 ENDIF RETURN end subroutine heat_DFLUX_231 !C*--------------------------------------------------------------------* SUBROUTINE heat_DFLUX_232(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT) !C*--------------------------------------------------------------------* !C*** !C*** SET 232 DFLUX !C*** !C* LTYPE=0 : BODY FLUX use hecmw IMPLICIT REAL(kind=kreal)(A-H,O-Z) !C* I/F VARIABLES INTEGER(kind=kint) NN, LTYPE REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),THICK,VAL,VECT(NN) REAL(kind=kreal) XG(3),WGT(3),H(6),HR(3) REAL(kind=kreal) HL1(6),HL2(6),HL3(6) INTEGER(kind=kint) NOD(3) !************************* ! GAUSS INTEGRATION POINT !************************* XG(1) = -0.7745966692 XG(2) = 0.0 XG(3) = 0.7745966692 WGT(1) = 0.5555555555 WGT(2) = 0.8888888888 WGT(3) = 0.5555555555 !C* IVOL=0 ISUF=0 IF( LTYPE.EQ.0 ) THEN IVOL=1 ELSE IF( LTYPE.GE.1 ) THEN ISUF=1 IF(LTYPE.EQ.1) THEN NOD(1)=1 NOD(2)=2 NOD(3)=4 ELSE IF(LTYPE.EQ.2) THEN NOD(1)=2 NOD(2)=3 NOD(3)=5 ELSE IF(LTYPE.EQ.3) THEN NOD(1)=3 NOD(2)=1 NOD(3)=6 ENDIF ENDIF DO I=1,NN VECT(I)=0.0 ENDDO !** SURFACE LOAD IF( ISUF.EQ.1 ) THEN !** INTEGRATION OVER SURFACE DO LX=1,3 RI=XG(LX) H(1)=-0.5*(1.0-RI)*RI H(2)= 0.5*(1.0+RI)*RI H(3)=1.0-RI*RI HR(1)= RI-0.5 HR(2)= RI+0.5 HR(3)=-2.0*RI GX=0.0 GY=0.0 DO I=1,3 GX=GX+HR(I)*XX(NOD(I)) GY=GY+HR(I)*YY(NOD(I)) ENDDO XSUM=DSQRT(GX*GX+GY*GY) WG=WGT(LX)*VAL*THICK VECT(NOD(1))=VECT(NOD(1))-H(1)*XSUM*WG VECT(NOD(2))=VECT(NOD(2))-H(2)*XSUM*WG VECT(NOD(3))=VECT(NOD(3))-H(3)*XSUM*WG ENDDO ENDIF !** VOLUME LOAD IF( IVOL.EQ.1 ) THEN !* LOOP OVER ALL INTEGRATION POINTS DO L2=1,3 XL2=XG(L2) X2 =(XL2+1.0)*0.5 DO L1=1,3 XL1=XG(L1) X1=0.5*(1.0-X2)*(XL1+1.0) ! INTERPOLATION FUNCTION X3=1.0-X1-X2 H(1)= X1*(2.0*X1-1.) H(2)= X2*(2.0*X2-1.) H(3)= X3*(2.0*X3-1.) H(4)= 4.0*X1*X2 H(5)= 4.0*X2*X3 H(6)= 4.0*X1*X3 ! DERIVATIVE OF INTERPOLATION FUNCTION ! FOR L1-COORDINATE HL1(1)=4.0*X1-1.0 HL1(2)= 0.0 HL1(3)= 0.0 HL1(4)= 4.0*X2 HL1(5)= 0.0 HL1(6)= 4.0*X3 ! FOR L2-COORDINATE HL2(1)= 0.0 HL2(2)= 4.0*X2-1.0 HL2(3)= 0.0 HL2(4)= 4.0*X1 HL2(5)= 4.0*X3 HL2(6)= 0.0 ! FOR L3-COORDINATE HL3(1)= 0.0 HL3(2)= 0.0 HL3(3)= 4.0*X3-1.0 HL3(4)= 0.0 HL3(5)= 4.0*X2 HL3(6)= 4.0*X1 ! JACOBI MATRIX XJ11=0.0 XJ21=0.0 XJ12=0.0 XJ22=0.0 DO I=1,NN XJ11=XJ11+(HL1(I)-HL3(I))*XX(I) XJ21=XJ21+(HL2(I)-HL3(I))*XX(I) XJ12=XJ12+(HL1(I)-HL3(I))*YY(I) XJ22=XJ22+(HL2(I)-HL3(I))*YY(I) ENDDO DET=XJ11*XJ22-XJ21*XJ12 WG=WGT(L1)*WGT(L2)*DET*THICK*(1.0-X2)*0.25 DO I = 1, NN VECT(I) = VECT(I) - H(I)*WG*VAL ENDDO ENDDO ENDDO ENDIF RETURN end subroutine heat_DFLUX_232 !C*--------------------------------------------------------------------* SUBROUTINE heat_DFLUX_241(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT) !C*--------------------------------------------------------------------* !C*** !C*** SET 241 DFLUX !C*** !C* LTYPE=0 : BODY FLUX use hecmw IMPLICIT REAL(kind=kreal)(A-H,O-Z) !C* I/F VARIABLES INTEGER(kind=kint) NN, LTYPE REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),THICK,VAL,VECT(NN) REAL(kind=kreal) XG(2),WGT(2),H(4),HR(4),HS(4) INTEGER(kind=kint) NOD(2) !************************* ! GAUSS INTEGRATION POINT !************************* DATA XG/-0.5773502691896,0.5773502691896/ DATA WGT/1.0,1.0/ !C* IVOL=0 ISUF=0 IF( LTYPE.EQ.0 ) THEN IVOL=1 ELSE IF( LTYPE.GE.1 ) THEN ISUF=1 IF(LTYPE.EQ.1) THEN NOD(1)=1 NOD(2)=2 ELSE IF(LTYPE.EQ.2) THEN NOD(1)=2 NOD(2)=3 ELSE IF(LTYPE.EQ.3) THEN NOD(1)=3 NOD(2)=4 ELSE IF(LTYPE.EQ.4) THEN NOD(1)=4 NOD(2)=1 ENDIF ENDIF DO I=1,NN VECT(I)=0.0 ENDDO !** SURFACE LOAD IF( ISUF.EQ.1 ) THEN !** INTEGRATION OVER SURFACE DO LX=1,2 RI=XG(LX) H(1)=0.5*(1.0-RI) H(2)=0.5*(1.0+RI) HR(1)=-0.5 HR(2)= 0.5 GX=0.0 GY=0.0 DO I=1,2 GX=GX+HR(I)*XX(NOD(I)) GY=GY+HR(I)*YY(NOD(I)) ENDDO XSUM=DSQRT(GX*GX+GY*GY)*THICK DO I=1,2 VECT(NOD(I))=VECT(NOD(I))-XSUM*WGT(LX)*H(I)*VAL ENDDO ENDDO ENDIF !** VOLUME LOAD IF( IVOL.EQ.1 ) THEN DO LX=1,2 RI=XG(LX) DO LY=1,2 SI=XG(LY) RP=1.0+RI SP=1.0+SI RM=1.0-RI SM=1.0-SI H(1)=0.25*RM*SM H(2)=0.25*RP*SM H(3)=0.25*RP*SP H(4)=0.25*RM*SP HR(1)=-.25*SM HR(2)= .25*SM HR(3)= .25*SP HR(4)=-.25*SP HS(1)=-.25*RM HS(2)=-.25*RP HS(3)= .25*RP HS(4)= .25*RM XJ11=0.0 XJ21=0.0 XJ12=0.0 XJ22=0.0 DO I=1,NN XJ11=XJ11+HR(I)*XX(I) XJ21=XJ21+HS(I)*XX(I) XJ12=XJ12+HR(I)*YY(I) XJ22=XJ22+HS(I)*YY(I) ENDDO DET=XJ11*XJ22-XJ21*XJ12 WG=WGT(LX)*WGT(LY)*DET*THICK DO I =1,NN VECT(I)=VECT(I)-WG*H(I)*VAL ENDDO ENDDO ENDDO ENDIF RETURN end subroutine heat_DFLUX_241 !C*--------------------------------------------------------------------* SUBROUTINE heat_DFLUX_242(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT) !C*--------------------------------------------------------------------* !C*** !C*** SET 242 DFLUX !C*** !C* LTYPE=0 : BODY FLUX use hecmw IMPLICIT REAL(kind=kreal)(A-H,O-Z) !C* I/F VARIABLES INTEGER(kind=kint) NN, LTYPE REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),THICK,VAL,VECT(NN) REAL(kind=kreal) XG(3),WGT(3),H(8),HR(8),HS(8) INTEGER(kind=kint) NOD(3) !************************* ! GAUSS INTEGRATION POINT !************************* XG(1) = -0.7745966692 XG(2) = 0.0 XG(3) = 0.7745966692 WGT(1) = 0.5555555555 WGT(2) = 0.8888888888 WGT(3) = 0.5555555555 !C* IVOL=0 ISUF=0 IF( LTYPE.EQ.0 ) THEN IVOL=1 ELSE IF( LTYPE.GE.1 ) THEN ISUF=1 IF(LTYPE.EQ.1) THEN NOD(1)=1 NOD(2)=2 NOD(3)=5 ELSE IF(LTYPE.EQ.2) THEN NOD(1)=2 NOD(2)=3 NOD(3)=6 ELSE IF(LTYPE.EQ.3) THEN NOD(1)=3 NOD(2)=4 NOD(3)=7 ELSE IF(LTYPE.EQ.4) THEN NOD(1)=4 NOD(2)=1 NOD(3)=8 ENDIF ENDIF DO I=1,NN VECT(I)=0.0 ENDDO !** SURFACE LOAD IF( ISUF.EQ.1 ) THEN !** INTEGRATION OVER SURFACE DO LX=1,3 RI=XG(LX) H(1)=-0.5*(1.0-RI)*RI H(2)= 0.5*(1.0+RI)*RI H(3)=1.0-RI*RI HR(1)= RI-0.5 HR(2)= RI+0.5 HR(3)=-2.0*RI GX=0.0 GY=0.0 DO I=1,3 GX=GX+HR(I)*XX(NOD(I)) GY=GY+HR(I)*YY(NOD(I)) ENDDO XSUM=DSQRT(GX*GX+GY*GY)*THICK VECT(NOD(1))=VECT(NOD(1))-H(1)*WGT(LX)*VAL*XSUM VECT(NOD(2))=VECT(NOD(2))-H(2)*WGT(LX)*VAL*XSUM VECT(NOD(3))=VECT(NOD(3))-H(3)*WGT(LX)*VAL*XSUM ENDDO ENDIF !** VOLUME LOAD IF( IVOL.EQ.1 ) THEN DO LX=1,3 RI=XG(LX) DO LY=1,3 SI=XG(LY) RP=1.0+RI SP=1.0+SI RM=1.0-RI SM=1.0-SI H(1)=0.25*RM*SM*(-1.0-RI-SI) H(2)=0.25*RP*SM*(-1.0+RI-SI) H(3)=0.25*RP*SP*(-1.0+RI+SI) H(4)=0.25*RM*SP*(-1.0-RI+SI) H(5)=0.5*(1.0-RI*RI)*(1.0-SI) H(6)=0.5*(1.0-SI*SI)*(1.0+RI) H(7)=0.5*(1.0-RI*RI)*(1.0+SI) H(8)=0.5*(1.0-SI*SI)*(1.0-RI) HR(1)=-.25*SM*(-1.0-RI-SI)-0.25*RM*SM HR(2)= .25*SM*(-1.0+RI-SI)+0.25*RP*SM HR(3)= .25*SP*(-1.0+RI+SI)+0.25*RP*SP HR(4)=-.25*SP*(-1.0-RI+SI)-0.25*RM*SP HR(5)=-RI*(1.0-SI) HR(6)= 0.5*(1.0-SI*SI) HR(7)=-RI*(1.0+SI) HR(8)=-0.5*(1.0-SI*SI) HS(1)=-.25*RM*(-1.0-RI-SI)-0.25*RM*SM HS(2)=-.25*RP*(-1.0+RI-SI)-0.25*RP*SM HS(3)= .25*RP*(-1.0+RI+SI)+0.25*RP*SP HS(4)= .25*RM*(-1.0-RI+SI)+0.25*RM*SP HS(5)=-0.5*(1.0-RI*RI) HS(6)=-SI *(1.0+RI) HS(7)= 0.5*(1.0-RI*RI) HS(8)= -SI*(1.0-RI) XJ11=0.0 XJ21=0.0 XJ12=0.0 XJ22=0.0 DO I=1,NN XJ11=XJ11+HR(I)*XX(I) XJ21=XJ21+HS(I)*XX(I) XJ12=XJ12+HR(I)*YY(I) XJ22=XJ22+HS(I)*YY(I) ENDDO DET=XJ11*XJ22-XJ21*XJ12 WG=WGT(LX)*WGT(LY)*DET*THICK DO I =1,NN VECT(I)=VECT(I)-WG*H(I)*VAL ENDDO ENDDO ENDDO ENDIF RETURN end subroutine heat_DFLUX_242 !----------------------------------------------------------------------* SUBROUTINE heat_DFLUX_341(NN,XX,YY,ZZ,LTYPE,VAL,VECT) !----------------------------------------------------------------------* !** !** SET 341 DFLUX !** ! BF LTYPE=0 :BODY FLUX ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4 use hecmw IMPLICIT NONE ! I/F VARIABLES INTEGER(kind=kint) NN,LTYPE REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),VECT(NN) ! LOCAL VARIABLES REAL(kind=kreal) H(4),HR(4),HS(4),HT(4) REAL(kind=kreal) XG(2),WGT(2) REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33,DET,WG INTEGER(kind=kint) IVOL,ISUF INTEGER(kind=kint) NOD(4) INTEGER(kind=kint) L1,L2,L3,I REAL(kind=kreal) VAL,AA REAL(kind=kreal) V1X,V1Y,V1Z REAL(kind=kreal) V2X,V2Y,V2Z REAL(kind=kreal) V3X,V3Y,V3Z REAL(kind=kreal) XL1,XL2,XL3 REAL(kind=kreal) X1,X2,X3 !************************* ! GAUSS INTEGRATION POINT !************************* DATA XG/-0.5773502691896,0.5773502691896/ DATA WGT/1.0,1.0/ ! IVOL=0 ISUF=0 IF( LTYPE.EQ.0 ) THEN IVOL=1 ELSE ISUF=1 IF(LTYPE.EQ.1) THEN NOD(1)=1 NOD(2)=2 NOD(3)=3 ELSE IF(LTYPE.EQ.2) THEN NOD(1)=4 NOD(2)=2 NOD(3)=1 ELSE IF(LTYPE.EQ.3) THEN NOD(1)=4 NOD(2)=3 NOD(3)=2 ELSE IF(LTYPE.EQ.4) THEN NOD(1)=4 NOD(2)=1 NOD(3)=3 ENDIF ENDIF ! CLEAR VECT DO I=1,NN VECT(I)=0.0 ENDDO !** SURFACE LOAD IF( ISUF.EQ.1 ) THEN V1X=XX(NOD(2))-XX(NOD(1)) V1Y=YY(NOD(2))-YY(NOD(1)) V1Z=ZZ(NOD(2))-ZZ(NOD(1)) V2X=XX(NOD(3))-XX(NOD(1)) V2Y=YY(NOD(3))-YY(NOD(1)) V2Z=ZZ(NOD(3))-ZZ(NOD(1)) V3X= V1Y*V2Z-V1Z*V2Y V3Y= V1Z*V2X-V1X*V2Z V3Z= V1X*V2Y-V1Y*V2X AA=0.5*DSQRT( V3X*v3X + V3Y*V3Y + V3Z*V3Z ) VECT(NOD(1))= -VAL*AA/3.0 VECT(NOD(2))= -VAL*AA/3.0 VECT(NOD(3))= -VAL*AA/3.0 ENDIF !** VOLUME LOAD IF( IVOL.EQ.1 ) THEN ! LOOP FOR INTEGRATION POINTS DO L3=1,2 XL3=XG(L3) X3 =(XL3+1.0)*0.5 DO L2=1,2 XL2=XG(L2) X2 =(1.0-X3)*(XL2+1.0)*0.5 DO L1=1,2 XL1=XG(L1) X1=(1.0-X2-X3)*(XL1+1.0)*0.5 ! INTERPOLATION FUNCTION H(1)=X1 H(2)=X2 H(3)=X3 H(4)=1.0-X1-X2-X3 ! DERIVATIVE OF INTERPOLATION FUNCTION ! FOR L1-COORDINATE HR(1)= 1.0 HR(2)= 0.0 HR(3)= 0.0 HR(4)=-1.0 ! FOR L2-COORDINATE HS(1)= 0.0 HS(2)= 1.0 HS(3)= 0.0 HS(4)=-1.0 ! FOR ZETA-COORDINATE HT(1)= 0.0 HT(2)= 0.0 HT(3)= 1.0 HT(4)=-1.0 ! JACOBI MATRIX XJ11=HR(1)*XX(1)+HR(2)*XX(2)+HR(3)*XX(3)+HR(4)*XX(4) XJ21=HS(1)*XX(1)+HS(2)*XX(2)+HS(3)*XX(3)+HS(4)*XX(4) XJ31=HT(1)*XX(1)+HT(2)*XX(2)+HT(3)*XX(3)+HT(4)*XX(4) XJ12=HR(1)*YY(1)+HR(2)*YY(2)+HR(3)*YY(3)+HR(4)*YY(4) XJ22=HS(1)*YY(1)+HS(2)*YY(2)+HS(3)*YY(3)+HS(4)*YY(4) XJ32=HT(1)*YY(1)+HT(2)*YY(2)+HT(3)*YY(3)+HT(4)*YY(4) XJ13=HR(1)*ZZ(1)+HR(2)*ZZ(2)+HR(3)*ZZ(3)+HR(4)*ZZ(4) XJ23=HS(1)*ZZ(1)+HS(2)*ZZ(2)+HS(3)*ZZ(3)+HS(4)*ZZ(4) XJ33=HT(1)*ZZ(1)+HT(2)*ZZ(2)+HT(3)*ZZ(3)+HT(4)*ZZ(4) ! DETERMINANT OF JACOBIAN DET=XJ11*XJ22*XJ33 & +XJ12*XJ23*XJ31 & +XJ13*XJ21*XJ32 & -XJ13*XJ22*XJ31 & -XJ12*XJ21*XJ33 & -XJ11*XJ23*XJ32 ! DO I = 1, NN VECT(I) = VECT(I) + VAL*H(I)*DET*(1.0-X3)*(1.0-X2-X3)*0.125 ENDDO enddo enddo enddo ENDIF ! RETURN end subroutine heat_DFLUX_341 !----------------------------------------------------------------------* SUBROUTINE heat_DFLUX_342(NN,XX,YY,ZZ,LTYPE,VAL,VECT & ,TT,IMAT,DTIME,k_elm,jj & ,ntab2,temp2,funcA2,funcB2 & ,idx_cure & ,ccc_rl & ,CURE3D_nod,idbg,incr,iterALL,D_TT) !----------------------------------------------------------------------* !** !** SET 342 DFLUX !** ! BF LTYPE=0 :BODY FLUX ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4 use hecmw use m_heat_LIB_CAPACITY IMPLICIT NONE ! I/F VARIABLES INTEGER(kind=kint) NN,LTYPE, IMAT,ntab2,idx_cure,k_elm,jj,idbg, incr,iterALL REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),VECT(NN),TT(NN),D_TT(NN) REAL(kind=kreal) temp2(ntab2),funcA2(ntab2+1),funcB2(ntab2+1) REAL(kind=kreal) DTIME, ccc_rl(*), Hr0,A1,A2,E1,E2,R0,m0,n0 REAL(kind=kreal) CURE3D_nod(NN),CURE3D_nod_I(NN),CURE3D_itg ! LOCAL VARIABLES REAL(kind=kreal) H(10) REAL(kind=kreal) HL1(10),HL2(10),HL3(10),HL4(10) REAL(kind=kreal) XG(3),WGT(3) REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33 REAL(kind=kreal) DET,WG INTEGER(kind=kint) IVOL,ISUF INTEGER(kind=kint) NOD(10) INTEGER(kind=kint) LX,LY,LZ,I,IG1,IG2,L1,L2,L3 REAL(kind=kreal) VAL,XSUM REAL(kind=kreal) G1X,G1Y,G1Z REAL(kind=kreal) G2X,G2Y,G2Z REAL(kind=kreal) G3X,G3Y,G3Z REAL(kind=kreal) XL1,XL2,XL3 REAL(kind=kreal) X1,X2,X3,X4 ! REAL(kind=kreal) DEN(3),CTEMP,D_CTEMP,dadt,K01,K02,K11,K22 REAL(kind=kreal) C3,C4,C40,D3,D4,D40,dtemp_dt,dtemp_dt2,LIM_d1 REAL(kind=kreal) time_reac1, REAC2(10) !!!! cure !!!!!!!!!!! Hr0 = ccc_rl( 1) A1 = ccc_rl( 2) A2 = ccc_rl( 3) E1 = ccc_rl( 4) E2 = ccc_rl( 5) R0 = ccc_rl( 6) m0 = ccc_rl( 7) n0 = ccc_rl( 8) !!!! cure !!!!!!!!!!! !************************* ! GAUSS INTEGRATION POINT !************************* XG(1) = -0.7745966692 XG(2) = 0.0 XG(3) = 0.7745966692 WGT(1) = 0.5555555555 WGT(2) = 0.8888888888 WGT(3) = 0.5555555555 ! IVOL=0 ISUF=0 IF( LTYPE.EQ.0 ) THEN IVOL=1 ELSE ISUF=1 IF(LTYPE.EQ.1) THEN NOD(1)=1 NOD(2)=2 NOD(3)=3 NOD(4)=5 NOD(5)=6 NOD(6)=7 ELSE IF(LTYPE.EQ.2) THEN NOD(1)=4 NOD(2)=2 NOD(3)=1 NOD(4)=9 NOD(5)=5 NOD(6)=8 ELSE IF(LTYPE.EQ.3) THEN NOD(1)=4 NOD(2)=3 NOD(3)=2 NOD(4)=10 NOD(5)=6 NOD(6)=9 ELSE IF(LTYPE.EQ.4) THEN NOD(1)=4 NOD(2)=1 NOD(3)=3 NOD(4)=8 NOD(5)=7 NOD(6)=10 ENDIF ENDIF ! CLEAR VECT DO I=1,10 VECT(I)=0.0 ENDDO !** SURFACE LOAD IF( ISUF.EQ.1 ) THEN ! INTEGRATION OVER SURFACE DO L2=1,3 XL2=XG(L2) X2 =(XL2+1.0)*0.5 DO L1=1,3 XL1=XG(L1) X1=0.5*(1.0-X2)*(XL1+1.0) ! INTERPOLATION FUNCTION X3=1.0-X1-X2 H(1)= X1*(2.0*X1-1.) H(2)= X2*(2.0*X2-1.) H(3)= X3*(2.0*X3-1.) H(4)= 4.0*X1*X2 H(5)= 4.0*X2*X3 H(6)= 4.0*X1*X3 ! DERIVATIVE OF INTERPOLATION FUNCTION ! FOR L1-COORDINATE HL1(1)=4.0*X1-1.0 HL1(2)= 0.0 HL1(3)= 0.0 HL1(4)= 4.0*X2 HL1(5)= 0.0 HL1(6)= 4.0*X3 ! FOR L2-COORDINATE HL2(1)= 0.0 HL2(2)= 4.0*X2-1.0 HL2(3)= 0.0 HL2(4)= 4.0*X1 HL2(5)= 4.0*X3 HL2(6)= 0.0 ! FOR L3-COORDINATE HL3(1)= 0.0 HL3(2)= 0.0 HL3(3)= 4.0*X3-1.0 HL3(4)= 0.0 HL3(5)= 4.0*X2 HL3(6)= 4.0*X1 ! JACOBI MATRIX G1X=0.0 G1Y=0.0 G1Z=0.0 G2X=0.0 G2Y=0.0 G2Z=0.0 DO I=1,6 G1X=G1X+(HL1(I)-HL3(I))*XX(NOD(I)) G1Y=G1Y+(HL1(I)-HL3(I))*YY(NOD(I)) G1Z=G1Z+(HL1(I)-HL3(I))*ZZ(NOD(I)) G2X=G2X+(HL2(I)-HL3(I))*XX(NOD(I)) G2Y=G2Y+(HL2(I)-HL3(I))*YY(NOD(I)) G2Z=G2Z+(HL2(I)-HL3(I))*ZZ(NOD(I)) ENDDO G3X=G1Y*G2Z-G1Z*G2Y G3Y=G1Z*G2X-G1X*G2Z G3Z=G1X*G2Y-G1Y*G2X XSUM=DSQRT(G3X**2+G3Y**2+G3Z**2) G3X=G3X/XSUM G3Y=G3Y/XSUM G3Z=G3Z/XSUM ! JACOBI MATRIX XJ11=G1X XJ12=G1Y XJ13=G1Z XJ21=G2X XJ22=G2Y XJ23=G2Z XJ31=G3X XJ32=G3Y XJ33=G3Z !DETERMINANT OF JACOBIAN DET=XJ11*XJ22*XJ33 & +XJ12*XJ23*XJ31 & +XJ13*XJ21*XJ32 & -XJ13*XJ22*XJ31 & -XJ12*XJ21*XJ33 & -XJ11*XJ23*XJ32 WG=WGT(L1)*WGT(L2)*DET*(1.0-X2)*0.25 DO I=1,6 VECT(NOD(I))=VECT(NOD(I))-VAL*WG*H(I) ENDDO ! ENDDO ENDDO ENDIF !** VOLUME LOAD IF( IVOL.EQ.1 ) THEN ! LOOP FOR INTEGRATION POINTS !DETERMINANT OF JACOBIAN DO L3=1,3 XL3=XG(L3) X3 =(XL3+1.0)*0.5 DO L2=1,3 XL2=XG(L2) X2 =(1.0-X3)*(XL2+1.0)*0.5 DO L1=1,3 XL1=XG(L1) X1=(1.0-X2-X3)*(XL1+1.0)*0.5 !C* INTERPOLATION FUNCTION X4=1.0-X1-X2-X3 H(1)=X1*(2.0*X1-1.0) H(2)=X2*(2.0*X2-1.0) H(3)=X3*(2.0*X3-1.0) H(4)=X4*(2.0*X4-1.0) H(5)=4.0*X1*X2 H(6)=4.0*X2*X3 H(7)=4.0*X1*X3 H(8)=4.0*X1*X4 H(9)=4.0*X2*X4 H(10)=4.0*X3*X4 !C* DERIVATIVE OF INTERPOLATION FUNCTION !C* FOR L1-COORDINATE HL1(1)= 4.0*X1-1.0 HL1(2)= 0.0 HL1(3)= 0.0 HL1(4)= 0.0 HL1(5)= 4.0*X2 HL1(6)= 0.0 HL1(7)= 4.0*X3 HL1(8)= 4.0*X4 HL1(9)= 0.0 HL1(10)= 0.0 !C* FOR L2-COORDINATE HL2(1)= 0.0 HL2(2)= 4.0*X2-1.0 HL2(3)= 0.0 HL2(4)= 0.0 HL2(5)= 4.0*X1 HL2(6)= 4.0*X3 HL2(7)= 0.0 HL2(8)= 0.0 HL2(9)= 4.0*X4 HL2(10)= 0.0 !C* FOR L3-COORDINATE HL3(1)= 0.0 HL3(2)= 0.0 HL3(3)= 4.0*X3-1.0 HL3(4)= 0.0 HL3(5)= 0.0 HL3(6)= 4.0*X2 HL3(7)= 4.0*X1 HL3(8)= 0.0 HL3(9)= 0.0 HL3(10)= 4.0*X4 !C* FOR L4-COORDINATE HL4(1)= 0.0 HL4(2)= 0.0 HL4(3)= 0.0 HL4(4)= 4.0*X4-1.0 HL4(5)= 0.0 HL4(6)= 0.0 HL4(7)= 0.0 HL4(8)= 4.0*X1 HL4(9)= 4.0*X2 HL4(10)= 4.0*X3 !C* JACOBI MATRIX XJ11=0.0 XJ21=0.0 XJ31=0.0 XJ12=0.0 XJ22=0.0 XJ32=0.0 XJ13=0.0 XJ23=0.0 XJ33=0.0 DO I=1,NN XJ11=XJ11+(HL1(I)-HL4(I))*XX(I) XJ21=XJ21+(HL2(I)-HL4(I))*XX(I) XJ31=XJ31+(HL3(I)-HL4(I))*XX(I) XJ12=XJ12+(HL1(I)-HL4(I))*YY(I) XJ22=XJ22+(HL2(I)-HL4(I))*YY(I) XJ32=XJ32+(HL3(I)-HL4(I))*YY(I) XJ13=XJ13+(HL1(I)-HL4(I))*ZZ(I) XJ23=XJ23+(HL2(I)-HL4(I))*ZZ(I) XJ33=XJ33+(HL3(I)-HL4(I))*ZZ(I) ENDDO DET=XJ11*XJ22*XJ33 & +XJ12*XJ23*XJ31 & +XJ13*XJ21*XJ32 & -XJ13*XJ22*XJ31 & -XJ12*XJ21*XJ33 & -XJ11*XJ23*XJ32 DET=-DET WG=DET*WGT(L1)*WGT(L2)*WGT(L3)*(1.-X3)*(1.0-X2-X3)*0.125 !!----- oga2009 ------------------------------------------------- if(idx_cure .eq. 0) then DO I = 1, NN VECT(I) = VECT(I) - VAL*WG*H(I) ENDDO else if(idx_cure .eq. 1) then if(IMAT .ne. 2) then VAL = 0.0D+00 DO I = 1, NN VECT(I) = VECT(I) - VAL*WG*H(I) ENDDO else if(IMAT .eq. 2) then ! CTEMP =0.0D+00 D_CTEMP =0.0D+00 CURE3D_itg=0.0D+00 DO I=1,NN CTEMP = CTEMP + TT(I)*H(I) D_CTEMP =D_CTEMP +D_TT(I)*H(I) CURE3D_itg=CURE3D_itg +CURE3D_nod(I)*H(I) END DO !!! if(CURE3D_itg .lt. 0.0D+00) then if(k_elm .eq. -999) then write(idbg,*) 'stop : CURE3D_itg .lt. 0.0 ' write(idbg,*) 'elem = ',k_elm,' CURE3D_itg= ',CURE3D_itg call flush(idbg) ! call hecmw_abort( hecmw_comm_get_comm()) end if CURE3D_itg = 0.0D+00 else if(CURE3D_itg .gt. 1.0D+00) then if(k_elm .eq. -999) then write(idbg,*) 'stop : CURE3D_itg .gt. 1.0 ' write(idbg,*) 'elem = ',k_elm,' CURE3D_itg= ',CURE3D_itg call flush(idbg) ! call hecmw_abort( hecmw_comm_get_comm()) end if CURE3D_itg = 1.0D+00 end if ! CALL heat_GET_DENSITY( CTEMP,IMAT,DEN,ntab2,temp2,funcA2,funcB2 ) K01 = A1*dexp(-E1/R0/(273.15D+00+CTEMP)) K02 = A2*dexp(-E2/R0/(273.15D+00+CTEMP)) dadt = (K01+K02*CURE3D_itg**m0)*(1.0D+00-CURE3D_itg)**n0 VAL = DEN(1)*Hr0*dadt DO I = 1, NN VECT(I) = VECT(I) - VAL*WG*H(I) ENDDO ! end if end if !!----- oga2009 ------------------------------------------------- ENDDO ENDDO ENDDO ENDIF ! RETURN end subroutine heat_DFLUX_342 !----------------------------------------------------------------------* SUBROUTINE heat_DFLUX_351(NN,XX,YY,ZZ,LTYPE,VAL,VECT) !----------------------------------------------------------------------* !** !** SET 351 DFLUX !** ! BF LTYPE=0 :BODY FLUX ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4 ! S5 LTYPE=5 :FLUX IN NORMAL-DIRECTION FOR FACE-5 use hecmw IMPLICIT REAL(kind=kreal)(A-H,O-Z) ! I/F VARIABLES INTEGER(kind=kint) NN, LTYPE REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),VAL,VECT(NN) ! LOCAL VARIABLES REAL(kind=kreal) H(6),HR(6),HS(6),HT(6),PL(6) REAL(kind=kreal) XG(2),WGT(2),XG1(3),XG2(3),WGT1(3) INTEGER(kind=kint) NOD(4) !************************* ! GAUSS INTEGRATION POINT !************************* DATA XG / -0.5773502691896,0.5773502691896 / DATA WGT / 1.0, 1.0 / DATA XG1 / 0.6666666667, 0.16666666667, 0.16666666667 / DATA XG2 / 0.1666666667, 0.66666666667, 0.16666666667 / DATA WGT1/ 0.1666666667, 0.16666666667, 0.16666666667 / ! ! SELCTION OF LOAD TYPE IVOL=0 ISUF=0 IF( LTYPE.EQ.0 ) THEN IVOL=1 ELSE ISUF=1 IF( LTYPE.EQ.1 ) THEN NOD(1) = 1 NOD(2) = 2 NOD(3) = 3 ISUF=2 ELSE IF( LTYPE.EQ.2 ) THEN NOD(1) = 6 NOD(2) = 5 NOD(3) = 4 ISUF=2 ELSE IF( LTYPE.EQ.3 ) THEN NOD(1) = 4 NOD(2) = 5 NOD(3) = 2 NOD(4) = 1 ELSE IF( LTYPE.EQ.4 ) THEN NOD(1) = 5 NOD(2) = 6 NOD(3) = 3 NOD(4) = 2 ELSE IF( LTYPE.EQ.5 ) THEN NOD(1) = 6 NOD(2) = 4 NOD(3) = 1 NOD(4) = 3 ENDIF ENDIF DO I=1,NN VECT(I)=0.0 ENDDO ! !** SURFACE LOAD ! IF( ISUF.EQ.1 ) THEN DO IG2=1,2 SI=XG(IG2) DO IG1=1,2 RI=XG(IG1) ! H(1)=0.25*(1.0-RI)*(1.0-SI) H(2)=0.25*(1.0+RI)*(1.0-SI) H(3)=0.25*(1.0+RI)*(1.0+SI) H(4)=0.25*(1.0-RI)*(1.0+SI) HR(1)=-.25*(1.0-SI) HR(2)= .25*(1.0-SI) HR(3)= .25*(1.0+SI) HR(4)=-.25*(1.0+SI) HS(1)=-.25*(1.0-RI) HS(2)=-.25*(1.0+RI) HS(3)= .25*(1.0+RI) HS(4)= .25*(1.0-RI) ! G1X=0.0 G1Y=0.0 G1Z=0.0 G2X=0.0 G2Y=0.0 G2Z=0.0 ! DO I=1,4 G1X=G1X+HR(I)*XX(NOD(I)) G1Y=G1Y+HR(I)*YY(NOD(I)) G1Z=G1Z+HR(I)*ZZ(NOD(I)) G2X=G2X+HS(I)*XX(NOD(I)) G2Y=G2Y+HS(I)*YY(NOD(I)) G2Z=G2Z+HS(I)*ZZ(NOD(I)) ENDDO ! G3X=G1Y*G2Z-G1Z*G2Y G3Y=G1Z*G2X-G1X*G2Z G3Z=G1X*G2Y-G1Y*G2X XSUM=DSQRT(G3X**2+G3Y**2+G3Z**2) DO I=1,4 VECT(NOD(I))=VECT(NOD(I))-XSUM*VAL*WGT(IG1)*WGT(IG2)*H(I) ENDDO ! enddo enddo ELSEIF( ISUF.EQ.2 ) THEN V1X=XX(NOD(2))-XX(NOD(1)) V1Y=YY(NOD(2))-YY(NOD(1)) V1Z=ZZ(NOD(2))-ZZ(NOD(1)) V2X=XX(NOD(3))-XX(NOD(1)) V2Y=YY(NOD(3))-YY(NOD(1)) V2Z=ZZ(NOD(3))-ZZ(NOD(1)) V3X= V1Y*V2Z-V1Z*V2Y V3Y= V1Z*V2X-V1X*V2Z V3Z= V1X*V2Y-V1Y*V2X AA=0.5*DSQRT( V3X*v3X + V3Y*V3Y + V3Z*V3Z ) VECT(NOD(1))= -VAL*AA/3.0 VECT(NOD(2))= -VAL*AA/3.0 VECT(NOD(3))= -VAL*AA/3.0 ENDIF ! !** VOLUME LOAD ! IF( IVOL.EQ.1 ) THEN DO I=1,6 VECT(I)=0.0 enddo ! DO LZ=1,2 ZI=XG(LZ) DO L12=1,3 X1=XG1(L12) X2=XG2(L12) H(1)=X1*(1.0-ZI)*0.5 H(2)=X2*(1.0-ZI)*0.5 H(3)=(1.0-X1-X2)*(1.0-ZI)*0.5 H(4)=X1*(1.0+ZI)*0.5 H(5)=X2*(1.0+ZI)*0.5 H(6)=(1.0-X1-X2)*(1.0+ZI)*0.5 HR(1)= (1.0-ZI)*0.5 HR(2)= 0.0 HR(3)=-(1.0-ZI)*0.5 HR(4)= (1.0+ZI)*0.5 HR(5)= 0.0 HR(6)=-(1.0+ZI)*0.5 HS(1)= 0.0 HS(2)= (1.0-ZI)*0.5 HS(3)= -(1.0-ZI)*0.5 HS(4)= 0.0 HS(5)= (1.0+ZI)*0.5 HS(6)= -(1.0+ZI)*0.5 HT(1)= -0.5*X1 HT(2)= -0.5*X2 HT(3)= -0.5*(1.0-X1-X2) HT(4)= 0.5*X1 HT(5)= 0.5*X2 HT(6)= 0.5*(1.0-X1-X2) ! XJ11=HR(1)*XX(1)+HR(2)*XX(2)+HR(3)*XX(3)+HR(4)*XX(4) & +HR(5)*XX(5)+HR(6)*XX(6) XJ21=HS(1)*XX(1)+HS(2)*XX(2)+HS(3)*XX(3)+HS(4)*XX(4) & +HS(5)*XX(5)+HS(6)*XX(6) XJ31=HT(1)*XX(1)+HT(2)*XX(2)+HT(3)*XX(3)+HT(4)*XX(4) & +HT(5)*XX(5)+HT(6)*XX(6) XJ12=HR(1)*YY(1)+HR(2)*YY(2)+HR(3)*YY(3)+HR(4)*YY(4) & +HR(5)*YY(5)+HR(6)*YY(6) XJ22=HS(1)*YY(1)+HS(2)*YY(2)+HS(3)*YY(3)+HS(4)*YY(4) & +HS(5)*YY(5)+HS(6)*YY(6) XJ32=HT(1)*YY(1)+HT(2)*YY(2)+HT(3)*YY(3)+HT(4)*YY(4) & +HT(5)*YY(5)+HT(6)*YY(6) XJ13=HR(1)*ZZ(1)+HR(2)*ZZ(2)+HR(3)*ZZ(3)+HR(4)*ZZ(4) & +HR(5)*ZZ(5)+HR(6)*ZZ(6) XJ23=HS(1)*ZZ(1)+HS(2)*ZZ(2)+HS(3)*ZZ(3)+HS(4)*ZZ(4) & +HS(5)*ZZ(5)+HS(6)*ZZ(6) XJ33=HT(1)*ZZ(1)+HT(2)*ZZ(2)+HT(3)*ZZ(3)+HT(4)*ZZ(4) & +HT(5)*ZZ(5)+HT(6)*ZZ(6) !*DETERMINANT OF JACOBIAN DET=XJ11*XJ22*XJ33 & +XJ12*XJ23*XJ31 & +XJ13*XJ21*XJ32 & -XJ13*XJ22*XJ31 & -XJ12*XJ21*XJ33 & -XJ11*XJ23*XJ32 WG=WGT1(L12)*WGT(LZ)*DET DO I=1,NN VECT(I)=VECT(I)-VAL*H(I)*WG ENDDO ! enddo enddo ENDIF ! RETURN end subroutine heat_DFLUX_351 !----------------------------------------------------------------------* SUBROUTINE heat_DFLUX_352(NN,XX,YY,ZZ,LTYPE,VAL,VECT) !----------------------------------------------------------------------* !** !** SET 352 DFLUX !** ! BF LTYPE=0 :BODY FLUX ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4 ! S5 LTYPE=5 :FLUX IN NORMAL-DIRECTION FOR FACE-5 use hecmw IMPLICIT REAL(kind=kreal)(A-H,O-Z) ! I/F VARIABLES INTEGER(kind=kint) NN, LTYPE REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),VAL,VECT(NN) ! LOCAL VARIABLES REAL(kind=kreal) H(15),HR(15),HS(15),HT(15),PL(15) REAL(kind=kreal) XG(3),WGT(3) INTEGER(kind=kint) NOD(8) !************************* ! GAUSS INTEGRATION POINT !************************* XG(1) = -0.7745966692 XG(2) = 0.0 XG(3) = 0.7745966692 WGT(1) = 0.5555555555 WGT(2) = 0.8888888888 WGT(3) = 0.5555555555 ! ! SELCTION OF LOAD TYPE ! IVOL=0 ISUF=0 IF( LTYPE.EQ.0 ) THEN IVOL=1 ELSE ISUF=1 IF( LTYPE.EQ.1 ) THEN NOD(1) = 1 NOD(2) = 2 NOD(3) = 3 NOD(4) = 7 NOD(5) = 8 NOD(6) = 9 ISUF=2 ELSE IF( LTYPE.EQ.2 ) THEN NOD(1) = 6 NOD(2) = 5 NOD(3) = 4 NOD(4) = 11 NOD(5) = 10 NOD(6) = 12 ISUF=2 ELSE IF( LTYPE.EQ.3 ) THEN NOD(1) = 4 NOD(2) = 5 NOD(3) = 2 NOD(4) = 1 NOD(5) = 10 NOD(6) = 14 NOD(7) = 7 NOD(8) = 13 ELSE IF( LTYPE.EQ.4 ) THEN NOD(1) = 5 NOD(2) = 6 NOD(3) = 3 NOD(4) = 2 NOD(5) = 11 NOD(6) = 15 NOD(7) = 8 NOD(8) = 14 ELSE IF( LTYPE.EQ.5 ) THEN NOD(1) = 6 NOD(2) = 4 NOD(3) = 1 NOD(4) = 3 NOD(5) = 12 NOD(6) = 13 NOD(7) = 9 NOD(8) = 15 ENDIF ENDIF DO I=1,NN VECT(I)=0.0 ENDDO ! !** SURFACE LOAD ! IF( ISUF.EQ.1 ) THEN DO IG2=1,3 SI=XG(IG2) DO IG1=1,3 RI=XG(IG1) RP=1.0+RI SP=1.0+SI RM=1.0-RI SM=1.0-SI H(1)=0.25*RM*SM*(-1.0-RI-SI) H(2)=0.25*RP*SM*(-1.0+RI-SI) H(3)=0.25*RP*SP*(-1.0+RI+SI) H(4)=0.25*RM*SP*(-1.0-RI+SI) H(5)=0.5*(1.0-RI*RI)*(1.0-SI) H(6)=0.5*(1.0-SI*SI)*(1.0+RI) H(7)=0.5*(1.0-RI*RI)*(1.0+SI) H(8)=0.5*(1.0-SI*SI)*(1.0-RI) HR(1)=-.25*SM*(-1.0-RI-SI)-0.25*RM*SM HR(2)= .25*SM*(-1.0+RI-SI)+0.25*RP*SM HR(3)= .25*SP*(-1.0+RI+SI)+0.25*RP*SP HR(4)=-.25*SP*(-1.0-RI+SI)-0.25*RM*SP HR(5)=-RI*(1.0-SI) HR(6)= .5*(1.0-SI*SI) HR(7)=-RI*(1.0+SI) HR(8)=-.5*(1.0-SI*SI) HS(1)=-.25*RM*(-1.0-RI-SI)-0.25*RM*SM HS(2)=-.25*RP*(-1.0+RI-SI)-0.25*RP*SM HS(3)= .25*RP*(-1.0+RI+SI)+0.25*RP*SP HS(4)= .25*RM*(-1.0-RI+SI)+0.25*RM*SP HS(5)=-.5*(1.0-RI*RI) HS(6)=-SI*(1.0+RI) HS(7)= .5*(1.0-RI*RI) HS(8)=-SI*(1.0-RI) G1X=0.0 G1Y=0.0 G1Z=0.0 G2X=0.0 G2Y=0.0 G2Z=0.0 G3X=0.0 G3Y=0.0 G3Z=0.0 DO I=1,8 G1X=G1X+HR(I)*XX(NOD(I)) G1Y=G1Y+HR(I)*YY(NOD(I)) G1Z=G1Z+HR(I)*ZZ(NOD(I)) G2X=G2X+HS(I)*XX(NOD(I)) G2Y=G2Y+HS(I)*YY(NOD(I)) G2Z=G2Z+HS(I)*ZZ(NOD(I)) G3X=G3X+HT(I)*XX(NOD(I)) G3Y=G3Y+HT(I)*YY(NOD(I)) G3Z=G3Z+HT(I)*ZZ(NOD(I)) ENDDO G3X=G1Y*G2Z-G1Z*G2Y G3Y=G1Z*G2X-G1X*G2Z G3Z=G1X*G2Y-G1Y*G2X XSUM=DSQRT(G3X**2+G3Y**2+G3Z**2) DO I=1,8 VECT(NOD(I))=VECT(NOD(I))-XSUM*VAL*WGT(IG1)*WGT(IG2)*H(I) ENDDO ENDDO ENDDO ELSEIF( ISUF.EQ.2 ) THEN DO IG2=1,3 SI=XG(IG2) DO IG1=1,3 RI=XG(IG1) RP=1.0+RI SP=1.0+SI RM=1.0-RI SM=1.0-SI H(1)=0.25*RM*SM*(-1.0-RI-SI) H(2)=0.25*RP*SM*(-1.0+RI-SI) H(3)=0.5*SP*SI H(4)=0.5*(1.0-RI*RI)*(1.0-SI) H(5)=0.5*(1.0-SI*SI)*(1.0+RI) H(6)=0.5*(1.0-SI*SI)*(1.0-RI) HR(1)=-.25*SM*(-1.0-RI-SI)-0.25*RM*SM HR(2)= .25*SM*(-1.0+RI-SI)+0.25*RP*SM HR(3)= 0.0 HR(4)=-RI*(1.0-SI) HR(5)= 0.5*(1.0-SI*SI) HR(6)=-0.5*(1.0-SI*SI) HS(1)=-.25*RM*(-1.0-RI-SI)-0.25*RM*SM HS(2)=-.25*RP*(-1.0+RI-SI)-0.25*RP*SM HS(3)=0.5*(1.0+2.0*SI) HS(4)=-0.5*(1.0-RI*RI) HS(5)=-SI *(1.0+RI) HS(6)=-SI*(1.0-RI) G1X=0.0 G1Y=0.0 G1Z=0.0 G2X=0.0 G2Y=0.0 G2Z=0.0 G3X=0.0 G3Y=0.0 G3Z=0.0 DO I=1, 6 G1X=G1X+HR(I)*XX(NOD(I)) G1Y=G1Y+HR(I)*YY(NOD(I)) G1Z=G1Z+HR(I)*ZZ(NOD(I)) G2X=G2X+HS(I)*XX(NOD(I)) G2Y=G2Y+HS(I)*YY(NOD(I)) G2Z=G2Z+HS(I)*ZZ(NOD(I)) G3X=G3X+HT(I)*XX(NOD(I)) G3Y=G3Y+HT(I)*YY(NOD(I)) G3Z=G3Z+HT(I)*ZZ(NOD(I)) ENDDO G3X=G1Y*G2Z-G1Z*G2Y G3Y=G1Z*G2X-G1X*G2Z G3Z=G1X*G2Y-G1Y*G2X XSUM=DSQRT(G3X**2+G3Y**2+G3Z**2) DO I=1,6 VECT(NOD(I))=VECT(NOD(I))-XSUM*VAL*WGT(IG1)*WGT(IG2)*H(I) ENDDO ENDDO ENDDO ENDIF ! !** VOLUME LOAD ! IF( IVOL.EQ.1 ) THEN DO LX=1,3 RI=XG(LX) DO LY=1,3 SI=XG(LY) DO LZ=1,3 TI=XG(LZ) RP=1.0+RI SP=1.0+SI TP=1.0+TI RM=1.0-RI SM=1.0-SI TM=1.0-TI ! INTERPOLATION FUNCTION H(1)=-0.125*RM*SM*TM*(2.0+RI+SI+TI) H(2)=-0.125*RP*SM*TM*(2.0-RI+SI+TI) H(3)=-0.25*SP*TM*(1.0-SI+TI) H(4)=-0.125*RM*SM*TP*(2.0+RI+SI-TI) H(5)=-0.125*RP*SM*TP*(2.0-RI+SI-TI) H(6)=0.25*SP*TP*(SI+TI-1.0) H(7)=0.25*(1.0-RI**2)*SM*TM H(8)=0.25*RP*(1.0-SI**2)*TM H(9)=0.25*RM*(1.0-SI**2)*TM H(10)=0.25*(1.0-RI**2)*SM*TP H(11)=0.25*RP*(1.0-SI**2)*TP H(12)=0.25*RM*(1.0-SI**2)*TP H(13)=0.25*RM*SM*(1.0-TI**2) H(14)=0.25*RP*SM*(1.0-TI**2) H(15)=0.5*SP*(1.0-TI**2) ! DERIVATIVE OF INTERPOLATION FUNCTION ! FOR R-COORDINATE HR(1)=-0.125*RM*SM*TM+0.125*SM*TM*(2.0+RI+SI+TI) HR(2)=+0.125*RP*SM*TM-0.125*SM*TM*(2.0-RI+SI+TI) HR(3)=0.0 HR(4)=-0.125*RM*SM*TP+0.125*SM*TP*(2.0+RI+SI-TI) HR(5)=+0.125*RP*SM*TP-0.125*SM*TP*(2.0-RI+SI-TI) HR(6)=0.0 HR(7)=-0.50*RI*SM*TM HR(8)=+0.25*(1.0-SI**2)*TM HR(9)=-0.25*(1.0-SI**2)*TM HR(10)=-0.50*RI*SM*TP HR(11)=+0.25*(1.0-SI**2)*TP HR(12)=-0.25*(1.0-SI**2)*TP HR(13)=-0.25*SM*(1.0-TI**2) HR(14)=+0.25*SM*(1.0-TI**2) HR(15)=0.0 ! FOR S-COORDINATE HS(1)=-0.125*RM*SM*TM+0.125*RM*TM*(2.0+RI+SI+TI) HS(2)=-0.125*RP*SM*TM+0.125*RP*TM*(2.0-RI+SI+TI) HS(3)=+0.25*TM*(2.0*SI-TI) HS(4)=-0.125*RM*SM*TP+0.125*RM*TP*(2.0+RI+SI-TI) HS(5)=-0.125*RP*SM*TP+0.125*RP*TP*(2.0-RI+SI-TI) HS(6)=+0.25*TP*(2.0*SI+TI) HS(7)=-0.25*(1.0-RI**2)*TM HS(8)=-0.50*RP*SI*TM HS(9)=-0.50*RM*SI*TM HS(10)=-0.25*(1.0-RI**2)*TP HS(11)=-0.50*RP*SI*TP HS(12)=-0.50*RM*SI*TP HS(13)=-0.25*RM*(1.0-TI**2) HS(14)=-0.25*RP*(1.0-TI**2) HS(15)=+0.50*(1.0-TI**2) ! FOR T-COORDINATE HT(1)=-0.125*RM*SM*TM+0.125*RM*SM*(2.0+RI+SI+TI) HT(2)=-0.125*RP*SM*TM+0.125*RP*SM*(2.0-RI+SI+TI) HT(3)=+0.25*SP*(2.0*TI-SI) HT(4)=+0.125*RM*SM*TP-0.125*RM*SM*(2.0+RI+SI-TI) HT(5)=+0.125*RP*SM*TP-0.125*RP*SM*(2.0-RI+SI-TI) HT(6)=+0.25*SP*(SI+2.0*TI) HT(7)=-0.25*(1.0-RI**2)*SM HT(8)=-0.25*RP*(1.0-SI**2) HT(9)=-0.25*RM*(1.0-SI**2) HT(10)=0.25*(1.0-RI**2)*SM HT(11)=0.25*RP*(1.0-SI**2) HT(12)=0.25*RM*(1.0-SI**2) HT(13)=-0.5*RM*SM*TI HT(14)=-0.5*RP*SM*TI HT(15)=-SP*TI ! JACOBI MATRIX XJ11=0.0 XJ21=0.0 XJ31=0.0 XJ12=0.0 XJ22=0.0 XJ32=0.0 XJ13=0.0 XJ23=0.0 XJ33=0.0 DO I=1,NN XJ11=XJ11+HR(I)*XX(I) XJ21=XJ21+HS(I)*XX(I) XJ31=XJ31+HT(I)*XX(I) XJ12=XJ12+HR(I)*YY(I) XJ22=XJ22+HS(I)*YY(I) XJ32=XJ32+HT(I)*YY(I) XJ13=XJ13+HR(I)*ZZ(I) XJ23=XJ23+HS(I)*ZZ(I) XJ33=XJ33+HT(I)*ZZ(I) ENDDO !DETERMINANT OF JACOBIAN DET=XJ11*XJ22*XJ33 & +XJ12*XJ23*XJ31 & +XJ13*XJ21*XJ32 & -XJ13*XJ22*XJ31 & -XJ12*XJ21*XJ33 & -XJ11*XJ23*XJ32 WG=DET*WGT(LX)*WGT(LY)*WGT(LZ) DO I=1,NN VECT(I)=VECT(I) - VAL*WG*H(I) ENDDO ENDDO ENDDO ENDDO ENDIF RETURN end subroutine heat_DFLUX_352 !----------------------------------------------------------------------* SUBROUTINE heat_DFLUX_361(NN,XX,YY,ZZ,LTYPE,VAL,VECT) !----------------------------------------------------------------------* !** !** SET 361 DFLUX !** ! BF LTYPE=0 :BODY FLUX ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4 ! S5 LTYPE=5 :FLUX IN NORMAL-DIRECTION FOR FACE-5 ! S6 LTYPE=6 :FLUX IN NORMAL-DIRECTION FOR FACE-6 use hecmw IMPLICIT NONE ! I/F VARIABLES INTEGER(kind=kint) NN, LTYPE REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),VAL,VECT(NN) ! LOCAL VARIABLES REAL(kind=kreal) H(8),HR(8),HS(8),HT(8),PL(8) REAL(kind=kreal) XG(2),WGT(2) REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33 ,DET,WG INTEGER(kind=kint) IVOL,ISUF INTEGER(kind=kint) NOD(4) INTEGER(kind=kint) IG1,IG2,LX,LY,LZ,I REAL(kind=kreal) VX,VY,VZ REAL(kind=kreal) G1X,G1Y,G1Z REAL(kind=kreal) G2X,G2Y,G2Z REAL(kind=kreal) G3X,G3Y,G3Z REAL(kind=kreal) XSUM !************************* ! GAUSS INTEGRATION POINT !************************* DATA XG/-0.5773502691896,0.5773502691896/ DATA WGT/1.0,1.0/ ! ! SELCTION OF LOAD TYPE ! IVOL=0 ISUF=0 IF( LTYPE.EQ.0 ) THEN IVOL=1 ELSE ISUF=1 IF( LTYPE.EQ.1 ) THEN NOD(1) = 1 NOD(2) = 2 NOD(3) = 3 NOD(4) = 4 ELSE IF( LTYPE.EQ.2 ) THEN NOD(1) = 8 NOD(2) = 7 NOD(3) = 6 NOD(4) = 5 ELSE IF( LTYPE.EQ.3 ) THEN NOD(1) = 5 NOD(2) = 6 NOD(3) = 2 NOD(4) = 1 ELSE IF( LTYPE.EQ.4 ) THEN NOD(1) = 6 NOD(2) = 7 NOD(3) = 3 NOD(4) = 2 ELSE IF( LTYPE.EQ.5 ) THEN NOD(1) = 7 NOD(2) = 8 NOD(3) = 4 NOD(4) = 3 ELSE IF( LTYPE.EQ.6 ) THEN NOD(1) = 8 NOD(2) = 5 NOD(3) = 1 NOD(4) = 4 ENDIF ENDIF ! CLEAR VECT DO I=1,NN VECT(I)=0.0 ENDDO ! !** SURFACE LOAD ! IF( ISUF.EQ.1 ) THEN ! INTEGRATION OVER SURFACE DO IG2=1,2 SI=XG(IG2) DO IG1=1,2 RI=XG(IG1) H(1)=0.25*(1.0-RI)*(1.0-SI) H(2)=0.25*(1.0+RI)*(1.0-SI) H(3)=0.25*(1.0+RI)*(1.0+SI) H(4)=0.25*(1.0-RI)*(1.0+SI) HR(1)=-.25*(1.0-SI) HR(2)= .25*(1.0-SI) HR(3)= .25*(1.0+SI) HR(4)=-.25*(1.0+SI) HS(1)=-.25*(1.0-RI) HS(2)=-.25*(1.0+RI) HS(3)= .25*(1.0+RI) HS(4)= .25*(1.0-RI) G1X=0.0 G1Y=0.0 G1Z=0.0 G2X=0.0 G2Y=0.0 G2Z=0.0 DO I=1,4 G1X=G1X+HR(I)*XX(NOD(I)) G1Y=G1Y+HR(I)*YY(NOD(I)) G1Z=G1Z+HR(I)*ZZ(NOD(I)) G2X=G2X+HS(I)*XX(NOD(I)) G2Y=G2Y+HS(I)*YY(NOD(I)) G2Z=G2Z+HS(I)*ZZ(NOD(I)) ENDDO G3X=G1Y*G2Z-G1Z*G2Y G3Y=G1Z*G2X-G1X*G2Z G3Z=G1X*G2Y-G1Y*G2X XSUM=DSQRT(G3X**2+G3Y**2+G3Z**2) DO I=1,4 VECT(NOD(I))=VECT(NOD(I))-XSUM*VAL*WGT(IG1)*WGT(IG2)*H(I) ENDDO ENDDO ENDDO ENDIF ! !** VOLUME LOAD ! IF( IVOL.EQ.1 ) THEN DO I=1,NN PL(I)=0.0 ENDDO ! LOOP FOR INTEGRATION POINTS DO LX=1,2 RI=XG(LX) DO LY=1,2 SI=XG(LY) DO LZ=1,2 TI=XG(LZ) RP=1.0+RI SP=1.0+SI TP=1.0+TI RM=1.0-RI SM=1.0-SI TM=1.0-TI ! INTERPOLATION FUNCTION H(1)=0.125*RM*SM*TM H(2)=0.125*RP*SM*TM H(3)=0.125*RP*SP*TM H(4)=0.125*RM*SP*TM H(5)=0.125*RM*SM*TP H(6)=0.125*RP*SM*TP H(7)=0.125*RP*SP*TP H(8)=0.125*RM*SP*TP ! DERIVATIVE OF INTERPOLATION FUNCTION ! FOR R-COORDINATE HR(1)=-.125*SM*TM HR(2)= .125*SM*TM HR(3)= .125*SP*TM HR(4)=-.125*SP*TM HR(5)=-.125*SM*TP HR(6)= .125*SM*TP HR(7)= .125*SP*TP HR(8)=-.125*SP*TP ! FOR S-COORDINATE HS(1)=-.125*RM*TM HS(2)=-.125*RP*TM HS(3)= .125*RP*TM HS(4)= .125*RM*TM HS(5)=-.125*RM*TP HS(6)=-.125*RP*TP HS(7)= .125*RP*TP HS(8)= .125*RM*TP ! FOR T-COORDINATE HT(1)=-.125*RM*SM HT(2)=-.125*RP*SM HT(3)=-.125*RP*SP HT(4)=-.125*RM*SP HT(5)= .125*RM*SM HT(6)= .125*RP*SM HT(7)= .125*RP*SP HT(8)= .125*RM*SP ! JACOBI MATRIX XJ11=0.0 XJ21=0.0 XJ31=0.0 XJ12=0.0 XJ22=0.0 XJ32=0.0 XJ13=0.0 XJ23=0.0 XJ33=0.0 DO I=1,NN XJ11=XJ11+HR(I)*XX(I) XJ21=XJ21+HS(I)*XX(I) XJ31=XJ31+HT(I)*XX(I) XJ12=XJ12+HR(I)*YY(I) XJ22=XJ22+HS(I)*YY(I) XJ32=XJ32+HT(I)*YY(I) XJ13=XJ13+HR(I)*ZZ(I) XJ23=XJ23+HS(I)*ZZ(I) XJ33=XJ33+HT(I)*ZZ(I) ENDDO !DETERMINANT OF JACOBIAN DET=XJ11*XJ22*XJ33 & +XJ12*XJ23*XJ31 & +XJ13*XJ21*XJ32 & -XJ13*XJ22*XJ31 & -XJ12*XJ21*XJ33 & -XJ11*XJ23*XJ32 WG=WGT(LX)*WGT(LY)*WGT(LZ)*DET DO I=1,NN VECT(I)=VECT(I)-H(I)*WG*VAL ENDDO ENDDO ENDDO ENDDO ENDIF !* RETURN end subroutine heat_DFLUX_361 !----------------------------------------------------------------------* SUBROUTINE heat_DFLUX_362(NN,XX,YY,ZZ,LTYPE,VAL,VECT) !----------------------------------------------------------------------* !** !** SET 362 DFLUX !** ! BF LTYPE=0 :BODY FLUX ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4 ! S5 LTYPE=5 :FLUX IN NORMAL-DIRECTION FOR FACE-5 ! S6 LTYPE=6 :FLUX IN NORMAL-DIRECTION FOR FACE-6 use hecmw IMPLICIT NONE ! I/F VARIABLES INTEGER(kind=kint) NN, LTYPE REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),VAL,VECT(NN) ! LOCAL VARIABLES REAL(kind=kreal) H(20),HR(20),HS(20),HT(20) REAL(kind=kreal) XG(3),WGT(3) REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33 REAL(kind=kreal) DET,WG INTEGER(kind=kint) IVOL,ISUF INTEGER(kind=kint) NOD(8) INTEGER(kind=kint) IG1,IG2,LX,LY,LZ,I REAL(kind=kreal) VX,VY,VZ REAL(kind=kreal) G1X,G1Y,G1Z REAL(kind=kreal) G2X,G2Y,G2Z REAL(kind=kreal) G3X,G3Y,G3Z REAL(kind=kreal) XSUM !************************* ! GAUSS INTEGRATION POINT !************************* XG(1) = -0.7745966692 XG(2) = 0.0 XG(3) = 0.7745966692 WGT(1) = 0.5555555555 WGT(2) = 0.8888888888 WGT(3) = 0.5555555555 ! ! SELCTION OF LOAD TYPE ! IVOL=0 ISUF=0 IF( LTYPE.EQ.0 ) THEN IVOL=1 ELSE ISUF=1 IF( LTYPE.EQ.1 ) THEN NOD(1) = 1 NOD(2) = 2 NOD(3) = 3 NOD(4) = 4 NOD(5) = 9 NOD(6) = 10 NOD(7) = 11 NOD(8) = 12 ELSE IF( LTYPE.EQ.2 ) THEN NOD(1) = 8 NOD(2) = 7 NOD(3) = 6 NOD(4) = 5 NOD(5) = 15 NOD(6) = 14 NOD(7) = 13 NOD(8) = 16 ELSE IF( LTYPE.EQ.3 ) THEN NOD(1) = 5 NOD(2) = 6 NOD(3) = 2 NOD(4) = 1 NOD(5) = 13 NOD(6) = 18 NOD(7) = 9 NOD(8) = 17 ELSE IF( LTYPE.EQ.4 ) THEN NOD(1) = 6 NOD(2) = 7 NOD(3) = 3 NOD(4) = 2 NOD(5) = 14 NOD(6) = 19 NOD(7) = 10 NOD(8) = 18 ELSE IF( LTYPE.EQ.5 ) THEN NOD(1) = 7 NOD(2) = 8 NOD(3) = 4 NOD(4) = 3 NOD(5) = 15 NOD(6) = 20 NOD(7) = 11 NOD(8) = 19 ELSE IF( LTYPE.EQ.6 ) THEN NOD(1) = 8 NOD(2) = 5 NOD(3) = 1 NOD(4) = 4 NOD(5) = 16 NOD(6) = 17 NOD(7) = 12 NOD(8) = 20 ENDIF ENDIF ! CLEAR VECT DO I=1,NN VECT(I)=0.0 ENDDO ! !** SURFACE LOAD ! IF( ISUF.EQ.1 ) THEN ! INTEGRATION OVER SURFACE DO IG2=1,3 SI=XG(IG2) DO IG1=1,3 RI=XG(IG1) RP=1.0+RI SP=1.0+SI RM=1.0-RI SM=1.0-SI H(1)=0.25*RM*SM*(-1.0-RI-SI) H(2)=0.25*RP*SM*(-1.0+RI-SI) H(3)=0.25*RP*SP*(-1.0+RI+SI) H(4)=0.25*RM*SP*(-1.0-RI+SI) H(5)=0.5*(1.0-RI*RI)*(1.0-SI) H(6)=0.5*(1.0-SI*SI)*(1.0+RI) H(7)=0.5*(1.0-RI*RI)*(1.0+SI) H(8)=0.5*(1.0-SI*SI)*(1.0-RI) HR(1)=-.25*SM*(-1.0-RI-SI)-0.25*RM*SM HR(2)= .25*SM*(-1.0+RI-SI)+0.25*RP*SM HR(3)= .25*SP*(-1.0+RI+SI)+0.25*RP*SP HR(4)=-.25*SP*(-1.0-RI+SI)-0.25*RM*SP HR(5)=-RI*(1.0-SI) HR(6)= .5*(1.0-SI*SI) HR(7)=-RI*(1.0+SI) HR(8)=-.5*(1.0-SI*SI) HS(1)=-.25*RM*(-1.0-RI-SI)-0.25*RM*SM HS(2)=-.25*RP*(-1.0+RI-SI)-0.25*RP*SM HS(3)= .25*RP*(-1.0+RI+SI)+0.25*RP*SP HS(4)= .25*RM*(-1.0-RI+SI)+0.25*RM*SP HS(5)=-.5*(1.0-RI*RI) HS(6)=-SI*(1.0+RI) HS(7)= .5*(1.0-RI*RI) HS(8)=-SI*(1.0-RI) G1X=0.0 G1Y=0.0 G1Z=0.0 G2X=0.0 G2Y=0.0 G2Z=0.0 DO I=1,8 G1X=G1X+HR(I)*XX(NOD(I)) G1Y=G1Y+HR(I)*YY(NOD(I)) G1Z=G1Z+HR(I)*ZZ(NOD(I)) G2X=G2X+HS(I)*XX(NOD(I)) G2Y=G2Y+HS(I)*YY(NOD(I)) G2Z=G2Z+HS(I)*ZZ(NOD(I)) ENDDO G3X=G1Y*G2Z-G1Z*G2Y G3Y=G1Z*G2X-G1X*G2Z G3Z=G1X*G2Y-G1Y*G2X XSUM=DSQRT(G3X**2+G3Y**2+G3Z**2) DO I=1,8 VECT(NOD(I))=VECT(NOD(I))-XSUM*VAL*WGT(IG1)*WGT(IG2)*H(I) ENDDO ENDDO ENDDO ENDIF ! !** VOLUME LOAD ! IF( IVOL.EQ.1 ) THEN ! LOOP FOR INTEGRATION POINTS DO LX=1,3 RI=XG(LX) DO LY=1,3 SI=XG(LY) DO LZ=1,3 TI=XG(LZ) RP=1.0+RI SP=1.0+SI TP=1.0+TI RM=1.0-RI SM=1.0-SI TM=1.0-TI ! INTERPOLATION FUNCTION H(1)=-0.125*RM*SM*TM*(2.0+RI+SI+TI) H(2)=-0.125*RP*SM*TM*(2.0-RI+SI+TI) H(3)=-0.125*RP*SP*TM*(2.0-RI-SI+TI) H(4)=-0.125*RM*SP*TM*(2.0+RI-SI+TI) H(5)=-0.125*RM*SM*TP*(2.0+RI+SI-TI) H(6)=-0.125*RP*SM*TP*(2.0-RI+SI-TI) H(7)=-0.125*RP*SP*TP*(2.0-RI-SI-TI) H(8)=-0.125*RM*SP*TP*(2.0+RI-SI-TI) H(9)=0.25*(1.0-RI**2)*SM*TM H(10)=0.25*RP*(1.0-SI**2)*TM H(11)=0.25*(1.0-RI**2)*SP*TM H(12)=0.25*RM*(1.0-SI**2)*TM H(13)=0.25*(1.0-RI**2)*SM*TP H(14)=0.25*RP*(1.0-SI**2)*TP H(15)=0.25*(1.0-RI**2)*SP*TP H(16)=0.25*RM*(1.0-SI**2)*TP H(17)=0.25*RM*SM*(1.0-TI**2) H(18)=0.25*RP*SM*(1.0-TI**2) H(19)=0.25*RP*SP*(1.0-TI**2) H(20)=0.25*RM*SP*(1.0-TI**2) ! DERIVATIVE OF INTERPOLATION FUNCTION ! FOR R-COORDINATE HR(1)=-0.125*RM*SM*TM+0.125*SM*TM*(2.0+RI+SI+TI) HR(2)=+0.125*RP*SM*TM-0.125*SM*TM*(2.0-RI+SI+TI) HR(3)=+0.125*RP*SP*TM-0.125*SP*TM*(2.0-RI-SI+TI) HR(4)=-0.125*RM*SP*TM+0.125*SP*TM*(2.0+RI-SI+TI) HR(5)=-0.125*RM*SM*TP+0.125*SM*TP*(2.0+RI+SI-TI) HR(6)=+0.125*RP*SM*TP-0.125*SM*TP*(2.0-RI+SI-TI) HR(7)=+0.125*RP*SP*TP-0.125*SP*TP*(2.0-RI-SI-TI) HR(8)=-0.125*RM*SP*TP+0.125*SP*TP*(2.0+RI-SI-TI) HR(9 )=-0.50*RI*SM*TM HR(10)=+0.25*(1.0-SI**2)*TM HR(11)=-0.50*RI*SP*TM HR(12)=-0.25*(1.0-SI**2)*TM HR(13)=-0.50*RI*SM*TP HR(14)=+0.25*(1.0-SI**2)*TP HR(15)=-0.50*RI*SP*TP HR(16)=-0.25*(1.0-SI**2)*TP HR(17)=-0.25*SM*(1.0-TI**2) HR(18)=+0.25*SM*(1.0-TI**2) HR(19)=+0.25*SP*(1.0-TI**2) HR(20)=-0.25*SP*(1.0-TI**2) ! FOR S-COORDINATE HS(1)=-0.125*RM*SM*TM+0.125*RM*TM*(2.0+RI+SI+TI) HS(2)=-0.125*RP*SM*TM+0.125*RP*TM*(2.0-RI+SI+TI) HS(3)=+0.125*RP*SP*TM-0.125*RP*TM*(2.0-RI-SI+TI) HS(4)=+0.125*RM*SP*TM-0.125*RM*TM*(2.0+RI-SI+TI) HS(5)=-0.125*RM*SM*TP+0.125*RM*TP*(2.0+RI+SI-TI) HS(6)=-0.125*RP*SM*TP+0.125*RP*TP*(2.0-RI+SI-TI) HS(7)=+0.125*RP*SP*TP-0.125*RP*TP*(2.0-RI-SI-TI) HS(8)=+0.125*RM*SP*TP-0.125*RM*TP*(2.0+RI-SI-TI) HS(9)=-0.25*(1.0-RI**2)*TM HS(10)=-0.50*RP*SI*TM HS(11)=+0.25*(1.0-RI**2)*TM HS(12)=-0.50*RM*SI*TM HS(13)=-0.25*(1.0-RI**2)*TP HS(14)=-0.50*RP*SI*TP HS(15)=+0.25*(1.0-RI**2)*TP HS(16)=-0.50*RM*SI*TP HS(17)=-0.25*RM*(1.0-TI**2) HS(18)=-0.25*RP*(1.0-TI**2) HS(19)=+0.25*RP*(1.0-TI**2) HS(20)=+0.25*RM*(1.0-TI**2) ! FOR T-COORDINATE HT(1)=-0.125*RM*SM*TM+0.125*RM*SM*(2.0+RI+SI+TI) HT(2)=-0.125*RP*SM*TM+0.125*RP*SM*(2.0-RI+SI+TI) HT(3)=-0.125*RP*SP*TM+0.125*RP*SP*(2.0-RI-SI+TI) HT(4)=-0.125*RM*SP*TM+0.125*RM*SP*(2.0+RI-SI+TI) HT(5)=+0.125*RM*SM*TP-0.125*RM*SM*(2.0+RI+SI-TI) HT(6)=+0.125*RP*SM*TP-0.125*RP*SM*(2.0-RI+SI-TI) HT(7)=+0.125*RP*SP*TP-0.125*RP*SP*(2.0-RI-SI-TI) HT(8)=+0.125*RM*SP*TP-0.125*RM*SP*(2.0+RI-SI-TI) HT(9)=-0.25*(1.0-RI**2)*SM HT(10)=-0.25*RP*(1.0-SI**2) HT(11)=-0.25*(1.0-RI**2)*SP HT(12)=-0.25*RM*(1.0-SI**2) HT(13)=0.25*(1.0-RI**2)*SM HT(14)=0.25*RP*(1.0-SI**2) HT(15)=0.25*(1.0-RI**2)*SP HT(16)=0.25*RM*(1.0-SI**2) HT(17)=-0.5*RM*SM*TI HT(18)=-0.5*RP*SM*TI HT(19)=-0.5*RP*SP*TI HT(20)=-0.5*RM*SP*TI ! JACOBI MATRIX XJ11=0.0 XJ21=0.0 XJ31=0.0 XJ12=0.0 XJ22=0.0 XJ32=0.0 XJ13=0.0 XJ23=0.0 XJ33=0.0 DO I=1,NN XJ11=XJ11+HR(I)*XX(I) XJ21=XJ21+HS(I)*XX(I) XJ31=XJ31+HT(I)*XX(I) XJ12=XJ12+HR(I)*YY(I) XJ22=XJ22+HS(I)*YY(I) XJ32=XJ32+HT(I)*YY(I) XJ13=XJ13+HR(I)*ZZ(I) XJ23=XJ23+HS(I)*ZZ(I) XJ33=XJ33+HT(I)*ZZ(I) ENDDO !DETERMINANT OF JACOBIAN DET=XJ11*XJ22*XJ33 & +XJ12*XJ23*XJ31 & +XJ13*XJ21*XJ32 & -XJ13*XJ22*XJ31 & -XJ12*XJ21*XJ33 & -XJ11*XJ23*XJ32 WG=WGT(LX)*WGT(LY)*WGT(LZ)*DET DO I=1,NN VECT(I)=VECT(I)-H(I)*WG*VAL ENDDO ENDDO ENDDO ENDDO ENDIF !* RETURN end subroutine heat_DFLUX_362 !----------------------------------------------------------------------* SUBROUTINE heat_DFLUX_731(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT) !----------------------------------------------------------------------* !** !** SET S3 DFLUX !** ! BF LTYPE=0 :BODY FLUX ! S1 LTYPE=1 :SURFACE FLUX use hecmw IMPLICIT REAL(kind=kreal)(A-H,O-Z) ! I/F VARIABLES INTEGER(kind=kint) NN, LTYPE REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),VAL,VECT(NN) ! IF( LTYPE.EQ.0 ) THEN THICK = THICK ELSEIF( LTYPE.EQ.1 ) THEN THICK = 1.0d0 ELSE THICK = 0.0d0 ENDIF DO I=1,NN VECT(I)=0.0 ENDDO ! V1X=XX(2)-XX(1) V1Y=YY(2)-YY(1) V1Z=ZZ(2)-ZZ(1) V2X=XX(3)-XX(1) V2Y=YY(3)-YY(1) V2Z=ZZ(3)-ZZ(1) V3X= V1Y*V2Z-V1Z*V2Y V3Y= V1Z*V2X-V1X*V2Z V3Z= V1X*V2Y-V1Y*V2X AA=0.5*DSQRT( V3X*v3X + V3Y*V3Y + V3Z*V3Z ) VECT(1)= -VAL*AA*THICK/3.0 VECT(2)= -VAL*AA*THICK/3.0 VECT(3)= -VAL*AA*THICK/3.0 ! RETURN end subroutine heat_DFLUX_731 !----------------------------------------------------------------------* SUBROUTINE heat_DFLUX_741(NN,XX,YY,ZZ,THICK,LTYPE,VAL,VECT) !----------------------------------------------------------------------* !** !** SET S4 DFLUX !** ! BF LTYPE=0 : BODY FLUX ! S1 LTYPE=1 : SURFACE FLUX use hecmw IMPLICIT REAL(kind=kreal)(A-H,O-Z) ! I/F VARIABLES INTEGER(kind=kint) NN, LTYPE REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),VAL,VECT(NN) ! LOCAL VARIABLES REAL(kind=kreal) H(4),HR(4),HS(4),HT(4),PL(4) REAL(kind=kreal) XG(2),WGT(2) !************************* ! GAUSS INTEGRATION POINT !************************* DATA XG / -0.5773502691896,0.5773502691896 / DATA WGT / 1.0, 1.0 / ! ! SELCTION OF LOAD TYPE ! IF ( LTYPE.EQ.0 ) THEN THICK = 1.0d0 * THICK ELSEIF( LTYPE.EQ.1 ) THEN THICK = 1.0d0 ELSE THICK = 0.0d0 ENDIF ! DO I=1,NN VECT(I)=0.0 ENDDO ! DO IG2=1,2 SI=XG(IG2) DO IG1=1,2 RI=XG(IG1) ! H(1)=0.25*(1.0-RI)*(1.0-SI) H(2)=0.25*(1.0+RI)*(1.0-SI) H(3)=0.25*(1.0+RI)*(1.0+SI) H(4)=0.25*(1.0-RI)*(1.0+SI) HR(1)=-.25*(1.0-SI) HR(2)= .25*(1.0-SI) HR(3)= .25*(1.0+SI) HR(4)=-.25*(1.0+SI) HS(1)=-.25*(1.0-RI) HS(2)=-.25*(1.0+RI) HS(3)= .25*(1.0+RI) HS(4)= .25*(1.0-RI) ! G1X=0.0 G1Y=0.0 G1Z=0.0 G2X=0.0 G2Y=0.0 G2Z=0.0 DO I=1,NN G1X=G1X+HR(I)*XX(I) G1Y=G1Y+HR(I)*YY(I) G1Z=G1Z+HR(I)*ZZ(I) G2X=G2X+HS(I)*XX(I) G2Y=G2Y+HS(I)*YY(I) G2Z=G2Z+HS(I)*ZZ(I) ENDDO ! G3X=G1Y*G2Z-G1Z*G2Y G3Y=G1Z*G2X-G1X*G2Z G3Z=G1X*G2Y-G1Y*G2X XSUM=DSQRT(G3X**2+G3Y**2+G3Z**2) DO I=1,NN VECT(I)=VECT(I)-XSUM*VAL*WGT(IG1)*WGT(IG2)*H(I)*THICK ENDDO ! enddo enddo RETURN end subroutine heat_DFLUX_741 !----------------------------------------------------------------------* SUBROUTINE NRCV_342CURE(ILOG,XX,YY,ZZ,CURE,ARRAY) !----------------------------------------------------------------------* ! ! NODAL RECOVERY FOR 3D 10 NODE PLANE ELEMENT /TOTAL LAGRANGE ! ! use m_fstr use hecmw use gauss_integration IMPLICIT NONE ! GAUSS POINT INTEGER(kind=kint) NG PARAMETER(NG=3) ! I/F VARIABLES INTEGER(kind=kint) ILOG REAL(kind=kreal) XX(*),YY(*),ZZ(*),ARRAY(13,*) REAL(kind=kreal) CURE(NG,NG,NG) REAL(kind=kreal) EE,PP,ALP ! LOCAL VARIABLES INTEGER(kind=kint) NN,NDOF PARAMETER(NN=10,NDOF=3) INTEGER(kind=kint) I,J,K,L INTEGER(kind=kint) LX,LY,LZ INTEGER(kind=kint) L1,L2,L3 REAL(kind=kreal) X1,X2,X3,X4 REAL(kind=kreal) XL1,XL2,XL3 REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM REAL(kind=kreal) H(NN),HL1(NN),HL2(NN),HL3(NN),HL4(NN) REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33,DET,WG,DUM REAL(kind=kreal) XJI11,XJI21,XJI31,XJI12,XJI22,XJI32,XJI13,XJI23,XJI33 REAL(kind=kreal) DNDX,DNDY,DNDZ REAL(kind=kreal) DUDX,DUDY,DUDZ REAL(kind=kreal) DVDX,DVDY,DVDZ REAL(kind=kreal) DWDX,DWDY,DWDZ REAL(kind=kreal) BL0(6,NN*NDOF),BL1(6,NN*NDOF) REAL(kind=kreal) UX(NN),UY(NN),UZ(NN) REAL(kind=kreal) S11,S22,S33,S23,S13,S12 REAL(kind=kreal) DETF,FMAT(3,3),SMAT(3,3),WK(3,3),TMAT(3,3),STRAIN_L(6) REAL(kind=kreal) SGM(6),EPS(6),PS,SMISES REAL(kind=kreal) ARRAY_TEMP(13,NG*NG*NG) INTEGER(kind=kint) IC REAL(kind=kreal) HG(NG*NG*NG),RX(NN),SX(NN),TX(NN) !************************* ! EVALUATION POINTS !************************* DATA RX/ 1., 0., 0.,-1., & 1., 0., 1., & 0.,-1.,-1./ DATA SX/-1., 1., 0.,-1., & 0., 1.,-1., & -1., 0.,-1./ DATA TX/-1.,-1., 1.,-1., & -1., 0., 0., & -1.,-1., 0./ !! IC=0 ! LOOP FOR INTEGRATION POINTS DO L3=1,NG XL3=XG(NG,L3) X3 =(XL3+1.0)*0.5 DO L2=1,NG XL2=XG(NG,L2) X2 =(1.0-X3)*(XL2+1.0)*0.5 DO L1=1,NG XL1=XG(NG,L1) X1=(1.0-X2-X3)*(XL1+1.0)*0.5 ! INTERPOLATION FUNCTION X4=1.0-X1-X2-X3 !! IC=IC+1 ARRAY_TEMP(1,IC)=CURE(L1,L2,L3) ENDDO ENDDO ENDDO !** !** NODAL VALUE !** DO J=1,NN RI=RX(J)/XG(NG,NG) SI=SX(J)/XG(NG,NG) TI=TX(J)/XG(NG,NG) RP=1.0+RI SP=1.0+SI TP=1.0+TI RM=1.0-RI SM=1.0-SI TM=1.0-TI ! INTERPOLATION FUNCTION (GAUSS POINT) HG( 1)=-0.125*RM*SM*TM*(2.0+RI+SI+TI) HG( 3)=-0.125*RP*SM*TM*(2.0-RI+SI+TI) HG( 9)=-0.125*RP*SP*TM*(2.0-RI-SI+TI) HG( 7)=-0.125*RM*SP*TM*(2.0+RI-SI+TI) HG(19)=-0.125*RM*SM*TP*(2.0+RI+SI-TI) HG(21)=-0.125*RP*SM*TP*(2.0-RI+SI-TI) HG(27)=-0.125*RP*SP*TP*(2.0-RI-SI-TI) HG(25)=-0.125*RM*SP*TP*(2.0+RI-SI-TI) HG( 2)=0.25*(1.0-RI**2)*SM*TM HG( 6)=0.25*RP*(1.0-SI**2)*TM HG( 8)=0.25*(1.0-RI**2)*SP*TM HG( 4)=0.25*RM*(1.0-SI**2)*TM HG(20)=0.25*(1.0-RI**2)*SM*TP HG(24)=0.25*RP*(1.0-SI**2)*TP HG(26)=0.25*(1.0-RI**2)*SP*TP HG(22)=0.25*RM*(1.0-SI**2)*TP HG(10)=0.25*RM*SM*(1.0-TI**2) HG(12)=0.25*RP*SM*(1.0-TI**2) HG(18)=0.25*RP*SP*(1.0-TI**2) HG(16)=0.25*RM*SP*(1.0-TI**2) HG( 5)=0.0 HG(11)=0.0 HG(13)=0.0 HG(14)=0.0 HG(15)=0.0 HG(17)=0.0 HG(23)=0.0 !! DO I=1,1 ARRAY(I,J)=0.0 DO K=1,NG*NG*NG ARRAY(I,J)=ARRAY(I,J)+HG(K)*ARRAY_TEMP(I,K) ENDDO ENDDO ENDDO RETURN END SUBROUTINE NRCV_342CURE !----------------------------------------------------------------------* SUBROUTINE cure_DFLUX_342(NN,XX,YY,ZZ,LTYPE,VAL,VECT & ,TT,IMAT,DTIME,k_elm,jj & ,ntab2,temp2,funcA2,funcB2 & ,idx_cure & ,ccc_rl & ,CURE3D_nod,idbg,incr,iterCUR,D_TT) !----------------------------------------------------------------------* !** !** SET 342 DFLUX FOR CURE !** ! BF LTYPE=0 :BODY FLUX ! S1 LTYPE=1 :FLUX IN NORMAL-DIRECTION FOR FACE-1 ! S2 LTYPE=2 :FLUX IN NORMAL-DIRECTION FOR FACE-2 ! S3 LTYPE=3 :FLUX IN NORMAL-DIRECTION FOR FACE-3 ! S4 LTYPE=4 :FLUX IN NORMAL-DIRECTION FOR FACE-4 use hecmw use m_heat_LIB_CAPACITY IMPLICIT NONE ! I/F VARIABLES INTEGER(kind=kint) NN,LTYPE, IMAT,ntab2,idx_cure,k_elm,jj,idbg, incr,iterCUR REAL(kind=kreal) XX(NN),YY(NN),ZZ(NN),VECT(NN),TT(NN),D_TT(NN) REAL(kind=kreal) temp2(ntab2),funcA2(ntab2+1),funcB2(ntab2+1) REAL(kind=kreal) DTIME, ccc_rl(*), Hr0,A1,A2,E1,E2,R0,m0,n0 REAL(kind=kreal) CURE3D_nod(NN),CURE3D_nod_I(NN),CURE3D_itg ! LOCAL VARIABLES REAL(kind=kreal) H(10) REAL(kind=kreal) HL1(10),HL2(10),HL3(10),HL4(10) REAL(kind=kreal) XG(3),WGT(3) REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33 REAL(kind=kreal) DET,WG INTEGER(kind=kint) IVOL,ISUF INTEGER(kind=kint) NOD(10) INTEGER(kind=kint) LX,LY,LZ,I,IG1,IG2,L1,L2,L3 REAL(kind=kreal) VAL,XSUM REAL(kind=kreal) G1X,G1Y,G1Z REAL(kind=kreal) G2X,G2Y,G2Z REAL(kind=kreal) G3X,G3Y,G3Z REAL(kind=kreal) XL1,XL2,XL3 REAL(kind=kreal) X1,X2,X3,X4 ! REAL(kind=kreal) DEN(3),CTEMP,D_CTEMP,dadt,K01,K02,K11,K22 REAL(kind=kreal) C3,C4,C40,D3,D4,D40,dtemp_dt,dtemp_dt2,LIM_d1 REAL(kind=kreal) time_reac1, REAC2(10) ! !!!! cure !!!!!!!!!!! Hr0 = ccc_rl( 1) A1 = ccc_rl( 2) A2 = ccc_rl( 3) E1 = ccc_rl( 4) E2 = ccc_rl( 5) R0 = ccc_rl( 6) m0 = ccc_rl( 7) n0 = ccc_rl( 8) !!!! cure !!!!!!!!!!! !************************* ! GAUSS INTEGRATION POINT !************************* XG(1) = -0.7745966692 XG(2) = 0.0 XG(3) = 0.7745966692 WGT(1) = 0.5555555555 WGT(2) = 0.8888888888 WGT(3) = 0.5555555555 ! IVOL=0 ISUF=0 IF( LTYPE.EQ.0 ) THEN IVOL=1 ELSE ISUF=1 IF(LTYPE.EQ.1) THEN NOD(1)=1 NOD(2)=2 NOD(3)=3 NOD(4)=5 NOD(5)=6 NOD(6)=7 ELSE IF(LTYPE.EQ.2) THEN NOD(1)=4 NOD(2)=2 NOD(3)=1 NOD(4)=9 NOD(5)=5 NOD(6)=8 ELSE IF(LTYPE.EQ.3) THEN NOD(1)=4 NOD(2)=3 NOD(3)=2 NOD(4)=10 NOD(5)=6 NOD(6)=9 ELSE IF(LTYPE.EQ.4) THEN NOD(1)=4 NOD(2)=1 NOD(3)=3 NOD(4)=8 NOD(5)=7 NOD(6)=10 ENDIF ENDIF ! CLEAR VECT DO I=1,10 VECT(I)=0.0 ENDDO !** SURFACE LOAD IF( ISUF.EQ.1 ) THEN ! INTEGRATION OVER SURFACE DO L2=1,3 XL2=XG(L2) X2 =(XL2+1.0)*0.5 DO L1=1,3 XL1=XG(L1) X1=0.5*(1.0-X2)*(XL1+1.0) ! INTERPOLATION FUNCTION X3=1.0-X1-X2 H(1)= X1*(2.0*X1-1.) H(2)= X2*(2.0*X2-1.) H(3)= X3*(2.0*X3-1.) H(4)= 4.0*X1*X2 H(5)= 4.0*X2*X3 H(6)= 4.0*X1*X3 ! DERIVATIVE OF INTERPOLATION FUNCTION ! FOR L1-COORDINATE HL1(1)=4.0*X1-1.0 HL1(2)= 0.0 HL1(3)= 0.0 HL1(4)= 4.0*X2 HL1(5)= 0.0 HL1(6)= 4.0*X3 ! FOR L2-COORDINATE HL2(1)= 0.0 HL2(2)= 4.0*X2-1.0 HL2(3)= 0.0 HL2(4)= 4.0*X1 HL2(5)= 4.0*X3 HL2(6)= 0.0 ! FOR L3-COORDINATE HL3(1)= 0.0 HL3(2)= 0.0 HL3(3)= 4.0*X3-1.0 HL3(4)= 0.0 HL3(5)= 4.0*X2 HL3(6)= 4.0*X1 ! JACOBI MATRIX G1X=0.0 G1Y=0.0 G1Z=0.0 G2X=0.0 G2Y=0.0 G2Z=0.0 DO I=1,6 G1X=G1X+(HL1(I)-HL3(I))*XX(NOD(I)) G1Y=G1Y+(HL1(I)-HL3(I))*YY(NOD(I)) G1Z=G1Z+(HL1(I)-HL3(I))*ZZ(NOD(I)) G2X=G2X+(HL2(I)-HL3(I))*XX(NOD(I)) G2Y=G2Y+(HL2(I)-HL3(I))*YY(NOD(I)) G2Z=G2Z+(HL2(I)-HL3(I))*ZZ(NOD(I)) ENDDO G3X=G1Y*G2Z-G1Z*G2Y G3Y=G1Z*G2X-G1X*G2Z G3Z=G1X*G2Y-G1Y*G2X XSUM=DSQRT(G3X**2+G3Y**2+G3Z**2) G3X=G3X/XSUM G3Y=G3Y/XSUM G3Z=G3Z/XSUM ! JACOBI MATRIX XJ11=G1X XJ12=G1Y XJ13=G1Z XJ21=G2X XJ22=G2Y XJ23=G2Z XJ31=G3X XJ32=G3Y XJ33=G3Z !DETERMINANT OF JACOBIAN DET=XJ11*XJ22*XJ33 & +XJ12*XJ23*XJ31 & +XJ13*XJ21*XJ32 & -XJ13*XJ22*XJ31 & -XJ12*XJ21*XJ33 & -XJ11*XJ23*XJ32 WG=WGT(L1)*WGT(L2)*DET*(1.0-X2)*0.25 DO I=1,6 VECT(NOD(I))=VECT(NOD(I))-VAL*WG*H(I) ENDDO ! ENDDO ENDDO ENDIF !** VOLUME LOAD IF( IVOL.EQ.1 ) THEN ! LOOP FOR INTEGRATION POINTS !DETERMINANT OF JACOBIAN DO L3=1,3 XL3=XG(L3) X3 =(XL3+1.0)*0.5 DO L2=1,3 XL2=XG(L2) X2 =(1.0-X3)*(XL2+1.0)*0.5 DO L1=1,3 XL1=XG(L1) X1=(1.0-X2-X3)*(XL1+1.0)*0.5 !C* INTERPOLATION FUNCTION X4=1.0-X1-X2-X3 H(1)=X1*(2.0*X1-1.0) H(2)=X2*(2.0*X2-1.0) H(3)=X3*(2.0*X3-1.0) H(4)=X4*(2.0*X4-1.0) H(5)=4.0*X1*X2 H(6)=4.0*X2*X3 H(7)=4.0*X1*X3 H(8)=4.0*X1*X4 H(9)=4.0*X2*X4 H(10)=4.0*X3*X4 !C* DERIVATIVE OF INTERPOLATION FUNCTION !C* FOR L1-COORDINATE HL1(1)= 4.0*X1-1.0 HL1(2)= 0.0 HL1(3)= 0.0 HL1(4)= 0.0 HL1(5)= 4.0*X2 HL1(6)= 0.0 HL1(7)= 4.0*X3 HL1(8)= 4.0*X4 HL1(9)= 0.0 HL1(10)= 0.0 !C* FOR L2-COORDINATE HL2(1)= 0.0 HL2(2)= 4.0*X2-1.0 HL2(3)= 0.0 HL2(4)= 0.0 HL2(5)= 4.0*X1 HL2(6)= 4.0*X3 HL2(7)= 0.0 HL2(8)= 0.0 HL2(9)= 4.0*X4 HL2(10)= 0.0 !C* FOR L3-COORDINATE HL3(1)= 0.0 HL3(2)= 0.0 HL3(3)= 4.0*X3-1.0 HL3(4)= 0.0 HL3(5)= 0.0 HL3(6)= 4.0*X2 HL3(7)= 4.0*X1 HL3(8)= 0.0 HL3(9)= 0.0 HL3(10)= 4.0*X4 !C* FOR L4-COORDINATE HL4(1)= 0.0 HL4(2)= 0.0 HL4(3)= 0.0 HL4(4)= 4.0*X4-1.0 HL4(5)= 0.0 HL4(6)= 0.0 HL4(7)= 0.0 HL4(8)= 4.0*X1 HL4(9)= 4.0*X2 HL4(10)= 4.0*X3 !C* JACOBI MATRIX XJ11=0.0 XJ21=0.0 XJ31=0.0 XJ12=0.0 XJ22=0.0 XJ32=0.0 XJ13=0.0 XJ23=0.0 XJ33=0.0 DO I=1,NN XJ11=XJ11+(HL1(I)-HL4(I))*XX(I) XJ21=XJ21+(HL2(I)-HL4(I))*XX(I) XJ31=XJ31+(HL3(I)-HL4(I))*XX(I) XJ12=XJ12+(HL1(I)-HL4(I))*YY(I) XJ22=XJ22+(HL2(I)-HL4(I))*YY(I) XJ32=XJ32+(HL3(I)-HL4(I))*YY(I) XJ13=XJ13+(HL1(I)-HL4(I))*ZZ(I) XJ23=XJ23+(HL2(I)-HL4(I))*ZZ(I) XJ33=XJ33+(HL3(I)-HL4(I))*ZZ(I) ENDDO DET=XJ11*XJ22*XJ33 & +XJ12*XJ23*XJ31 & +XJ13*XJ21*XJ32 & -XJ13*XJ22*XJ31 & -XJ12*XJ21*XJ33 & -XJ11*XJ23*XJ32 DET=-DET WG=DET*WGT(L1)*WGT(L2)*WGT(L3)*(1.-X3)*(1.0-X2-X3)*0.125 !!----- oga2009 ------------------------------------------------- if(idx_cure .eq. 0) then VAL = 0.0D+00 DO I = 1, NN VECT(I) = VECT(I) - VAL*WG*H(I) ENDDO else if(idx_cure .eq. 1) then if(IMAT .ne. 2) then dadt = 0.0D+00 DEN(1) = 1.0D+00 VAL = DEN(1)*dadt DO I = 1, NN VECT(I) = VECT(I) - VAL*WG*H(I) ENDDO else if(IMAT .eq. 2) then ! CTEMP =0.0D+00 D_CTEMP =0.0D+00 CURE3D_itg=0.0D+00 DO I=1,NN CTEMP = CTEMP + TT(I)*H(I) D_CTEMP =D_CTEMP +D_TT(I)*H(I) CURE3D_itg=CURE3D_itg +CURE3D_nod(I)*H(I) END DO !!! if(CURE3D_itg .lt. 0.0D+00) then if(k_elm .eq. -999) then write(idbg,*) 'stop : CURE3D_itg .lt. 0.0 ' write(idbg,*) 'elem = ',k_elm,' CURE3D_itg= ',CURE3D_itg call flush(idbg) ! call hecmw_abort( hecmw_comm_get_comm()) end if CURE3D_itg = 0.0D+00 else if(CURE3D_itg .gt. 1.0D+00) then if(k_elm .eq. -999) then write(idbg,*) 'stop : CURE3D_itg .gt. 1.0 ' write(idbg,*) 'elem = ',k_elm,' CURE3D_itg= ',CURE3D_itg call flush(idbg) ! call hecmw_abort( hecmw_comm_get_comm()) end if CURE3D_itg = 1.0D+00 end if ! CALL heat_GET_DENSITY( CTEMP,IMAT,DEN,ntab2,temp2,funcA2,funcB2 ) K01 = A1*dexp(-E1/R0/(273.15D+00+CTEMP)) K02 = A2*dexp(-E2/R0/(273.15D+00+CTEMP)) dadt = (K01+K02*CURE3D_itg**m0)*(1.0D+00-CURE3D_itg)**n0 VAL = DEN(1)*dadt DO I = 1, NN VECT(I) = VECT(I) - VAL*WG*H(I) ENDDO ! end if end if !!----- oga2009 ------------------------------------------------- ENDDO ENDDO ENDDO ENDIF ! RETURN end subroutine cure_DFLUX_342 end module m_heat_LIB_DFLUX