SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,PROPS, 1 NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME,KSTEP,KINC, 2 JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF,LFLAGS, 3 MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD) C INCLUDE 'aba_param.inc' C DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),SVARS(NSVARS), 1 ENERGY(8),PROPS(*),COORDS(MCRD,NNODE), 2 U(NDOFEL),DU(MLVARX,*),V(NDOFEL),A(NDOFEL),TIME(2), 3 PARAMS(3),JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*), 4 DDLMAG(MDLOAD,*),PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*) C DIMENSION GPX(9),GPY(9),GWEI(9),PHI(8),PHIX(8),PHIY(8),PHIC(8), 1 PHIE(8),IFACE(9),GWE(3) C PARAMETER(ZERO=0.D0,TWOHUN=200.D0,FIVHUN=500.D0,CONDUC=204.D0) DATA IFACE/1,5,2,6,3,7,4,8,1/ C C 8 NODE CONTINUUM UEL FOR TRANSIENT HEAT TRANSFER ANALYSIS; C TEMPERATURE DEPENDENT CONDUCTIVITY AND THE UNSYMMETRIC C CONTRIBUTION OF THIS TERM TO THE JACOBIAN INCLUDED C C VARIABLE DECLARATION C C THICK : THICKNESS C RHO : DENSITY C SPEC : SPECIFIC HEAT C COND : THERMAL CONDUCTIVITY (MAY BE TEMP DEPENDENT) C HCOEF : FILM COEFFICIENT C T : TEMPERATURE AT TIME T + DELTA T C TOLD : TEMPERATURE AT TIME T C TSINK : SINK TEMPERATURE C DTDX : DERIVATIVE OF TEMPERATURE WRT X C DTDY : DERIVATIVE OF TEMPERATURE WRT Y C DTDT : DERIVATIVE OF TEMPERATURE WRT TO TIME C DCDT : DERIVATIVE OF CONDUCTIVITY WRT TO TEMPERATURE C C : ISOPARAMETRIC COORDINATE, XI C E : ISOPARAMETRIC COORDINATE, ETA C WE : GAUSS WEIGHT MULTIPLIED BY JACOBIAN OF TRANSFORMATION C PHI : INTERPOLATION FUNCTIONS C PHIX : DERIVATIVE OF PHI WRT X C PHIY : DERIVATIVE OF PHI WRT Y C PHIC : DERIVATIVE OF PHI WRT XI C PHIE : DERIVATIVE OF PHI WRT ETA C C MATERIAL PROPERTY DEFINITION C THICK = PROPS(1) RHO = PROPS(2) SPEC = PROPS(3) C C INITIALIZATION (NRHS=1) C DO 6 K1=1,NDOFEL RHS(K1,NRHS)=ZERO DO 4 K2=1,NDOFEL AMATRX(K2,K1)=ZERO 4 CONTINUE 6 CONTINUE C IF (LFLAGS(3).EQ.4) RETURN C C TRANSIENT ANALYSIS C IF (LFLAGS(1).EQ.33) THEN C C DETERMINE GAUSS POINT LOCATIONS C CALL GSPT(GPX,GPY) C C DETERMINE GAUSS WEIGHTS C CALL GSWT(GWEI,GWE) C C ASSEMBLE AMATRX AND RHS C DO 300 K=1,9 C LOOP THROUGH GAUSS PTS C=GPX(K) E=GPY(K) CALL DER(C,E,GPX,GPY,GWEI,PHI,PHIX,PHIY,PHIC,PHIE 1 ,DXDC,DXDE,DYDC,DYDE,AJACOB,COORDS,MCRD,NNODE) DTDX=ZERO DTDY=ZERO T =ZERO TOLD=ZERO DO I=1,8 DTDX=U(I)*PHIX(I)+DTDX DTDY=U(I)*PHIY(I)+DTDY T =U(I)*PHI(I)+T TOLD=(U(I)-DU(I,NRHS))*PHI(I)+TOLD END DO C CHECK ON TEMPERATURE DEPENDENCE OF THERMAL CONDUCTIVITY IF (NPROPS.EQ.4) THEN COND=CONDUC DCDT=ZERO ELSE CALL UCOND(T,COND,DCDT) END IF DTDT=(T-TOLD)/DTIME WE=GWEI(K)*AJACOB DO KI=1,8 C LOOP OVER NODES RHS(KI,NRHS) = RHS(KI,NRHS) - 1 WE*(PHI(KI)*RHO*SPEC*DTDT + 2 COND*(PHIX(KI)*DTDX + PHIY(KI)*DTDY)) DO KJ=1,8 AMATRX(KI,KJ)= AMATRX(KI,KJ) + WE*(PHI(KI)*PHI(KJ)*RHO* 1 SPEC/DTIME + COND*(PHIX(KI)*PHIX(KJ) + PHIY(KI)* 2 PHIY(KJ)) + DCDT*PHI(KJ)*(PHIX(KI)*DTDX + 3 PHIY(KI)*DTDY)) END DO END DO 300 CONTINUE C C FLUX BOUNDARY CONDITIONS C C APPLIED DISTRIBUTED FLUX C IF (JELEM.EQ.101.AND.JDLTYP(1,1).EQ.4) THEN C=-1. DO KI=1,3 C LOOP THROUGH GAUSS PTS E=GPY(KI) CALL DER(C,E,GPX,GPY,GWEI,PHI,PHIX,PHIY,PHIC,PHIE 1 ,DXDC,DXDE,DYDC,DYDE,AJACOB,COORDS,MCRD,NNODE) DS=SQRT(DXDE*DXDE + DYDE*DYDE) DO KJ=7,9 C LOOP THROUGH NODES RHS(IFACE(KJ),NRHS) = RHS(IFACE(KJ),NRHS) 1 + GWE(KI)*DS*PHI(IFACE(KJ))*ADLMAG(1,1) END DO END DO END IF C C FILM CONDITION C IF (JELEM.EQ.102) THEN TSINK=TWOHUN HCOEF=FIVHUN C=1. DO KI=1,3 C LOOP THROUGH GAUSS PTS E=GPY(KI) CALL DER(C,E,GPX,GPY,GWEI,PHI,PHIX,PHIY,PHIC,PHIE 1 ,DXDC,DXDE,DYDC,DYDE,AJACOB,COORDS,MCRD,NNODE) T=0.0 DO I=1,8 T=U(I)*PHI(I)+T END DO DS=SQRT(DXDE*DXDE + DYDE*DYDE) DO KJ=3,5 C LOOP THROUGH NODES RHS(IFACE(KJ),NRHS) = RHS(IFACE(KJ),NRHS) - 1 GWE(KI)*DS*PHI(IFACE(KJ))*HCOEF*(T-TSINK) DO KK=1,8 C LOOP THROUGH NODES AMATRX(IFACE(KJ),KK)= AMATRX(IFACE(KJ),KK) + 1 GWE(KI)*DS*PHI(IFACE(KJ))*HCOEF*PHI(KK) END DO END DO END DO END IF END IF RETURN END C C SUBROUTINE UCOND(T,C,DCDT) INCLUDE 'aba_param.inc' C C RETURNS CONDUCTIVITY (C) AND DERIVATIVE OF CONDUCTIVITY WITH C RESPECT TO TEMPERATURE (DCDT) FOR A GIVEN TEMPERATURE (T) C DIMENSION TABLE(2,5) C PARAMETER(ZERO=0.D0) DATA TABLE/204.D0,250.D0,225.D0,275.D0,235.D0,300.D0, 1 245.D0,325.D0,255.D0,350.D0/ C C TABLE(1,N): CONDUCTIVITY (N DATA PAIRS) C TABLE(2,N): TEMPERATURE (N DATA PAIRS) C C RESET INC=1 WHEN C AND DCDT ARE FOUND FOR GIVEN T C INC=0 C C ASSIGN C AND DCDT IF T IS LESS THAN MIN. TEMPERATURE IN TABLE C IF (T.LT.TABLE(2,1)) THEN C=TABLE(1,1) DCDT=ZERO RETURN END IF C C ASSIGN C AND DCDT IF T IS GREATER THAN MAX. TEMPERATURE IN TABLE C IF (T.GT.TABLE(2,5)) THEN C=TABLE(1,5) DCDT=ZERO RETURN END IF DO 10 K1=1,4 TL1=TABLE(2,K1+1) IF(T.LT.TL1.AND.INC.EQ.0) THEN TL0=TABLE(2,K1) DT=TL1-TL0 C0=TABLE(1,K1) C1=TABLE(1,K1+1) DC=C1-C0 DCDT=DC/DT C=DCDT*(T-TL0)+C0 INC=1 ENDIF 10 CONTINUE RETURN END C C SUBROUTINE GSPT(GPX,GPY) INCLUDE 'aba_param.inc' DIMENSION AR(3),GPX(9),GPY(9) C PARAMETER(ZERO=0.D0,ONENEG=-1.D0,ONE=1.D0,SIX=6.D0,TEN=10.D0) C C GPX: X COORDINATE OF GAUSS PT C GPY: Y COORDINATE OF GAUSS PT C R=SQRT(SIX/TEN) AR(1)=ONENEG AR(2)=ZERO AR(3)=ONE DO 10 I=1,3 DO 10 J=1,3 NUMGP=(I-1)*3+J GPX(NUMGP)=AR(I)*R GPY(NUMGP)=AR(J)*R 10 CONTINUE RETURN END C C SUBROUTINE GSWT(GWEI,GWE) INCLUDE 'aba_param.inc' DIMENSION GWEI(9),GWE(3) C PARAMETER(FIVE=5.D0,EIGHT=8.D0,NINE=9.D0) C C GWEI : GAUSS WEIGHT C GWE(1)=FIVE/NINE GWE(2)=EIGHT/NINE GWE(3)=FIVE/NINE DO 10 I=1,3 DO 10 J=1,3 NUMGP=(I-1)*3+J GWEI(NUMGP)=GWE(I)*GWE(J) 10 CONTINUE RETURN END C C SUBROUTINE DER(C,E,GPX,GPY,GWEI,PHI,PHIX,PHIY,PHIC,PHIE, 1 DXDC,DXDE,DYDC,DYDE,AJACOB,COORDS,MCRD,NNODE) INCLUDE 'aba_param.inc' DIMENSION PHI(8),PHIX(8),PHIY(8),PHIC(8),PHIE(8), 1 COORDS(MCRD,NNODE) C PARAMETER(ZERO=0.D0,FOURTH=0.25D0,HALF=0.5D0,ONE=1.D0,TWO=2.D0) C C INTERPOLATION FUNCTIONS C PHI(1) = FOURTH*(ONE-C)*(ONE-E)*(-C-E-ONE) PHI(2) = FOURTH*(ONE+C)*(ONE-E)*(C-E-ONE) PHI(3) = FOURTH*(ONE+C)*(ONE+E)*(C+E-ONE) PHI(4) = FOURTH*(ONE-C)*(ONE+E)*(-C+E-ONE) PHI(5) = HALF*(ONE-C*C)*(ONE-E) PHI(6) = HALF*(ONE+C)*(ONE-E*E) PHI(7) = HALF*(ONE-C*C)*(ONE+E) PHI(8) = HALF*(ONE-C)*(ONE-E*E) C C DERIVATIVES WRT TO C C PHIC(1) = FOURTH*(ONE-E)*(TWO*C+E) PHIC(2) = FOURTH*(ONE-E)*(TWO*C-E) PHIC(3) = FOURTH*(ONE+E)*(TWO*C+E) PHIC(4) = FOURTH*(ONE+E)*(TWO*C-E) PHIC(5) = -C*(ONE-E) PHIC(6) = HALF*(ONE-E*E) PHIC(7) = -C*(ONE+E) PHIC(8) = -HALF*(ONE-E*E) C C DERIVATIVES WRT TO E C PHIE(1) = FOURTH*(ONE-C)*(TWO*E+C) PHIE(2) = FOURTH*(ONE+C)*(TWO*E-C) PHIE(3) = FOURTH*(ONE+C)*(TWO*E+C) PHIE(4) = FOURTH*(ONE-C)*(TWO*E-C) PHIE(5) = -HALF*(ONE-C*C) PHIE(6) = -E*(ONE+C) PHIE(7) = HALF*(ONE-C*C) PHIE(8) = -E*(ONE-C) DXDC=ZERO DXDE=ZERO DYDC=ZERO DYDE=ZERO DO 3 I=1,8 DXDC=DXDC+COORDS(1,I)*PHIC(I) DXDE=DXDE+COORDS(1,I)*PHIE(I) DYDC=DYDC+COORDS(2,I)*PHIC(I) DYDE=DYDE+COORDS(2,I)*PHIE(I) 3 CONTINUE C C CALCULATION OF JACOBIAN C AJACOB=(DXDC*DYDE-DXDE*DYDC) C C DERIVATIVES WRT TO X AND Y C DO 5 I=1,8 PHIX(I)=(PHIC(I)*DYDE-PHIE(I)*DYDC)/AJACOB PHIY(I)=(PHIE(I)*DXDC-PHIC(I)*DXDE)/AJACOB 5 CONTINUE RETURN END