C KELVIN SOURCE PV 17/12/08 21:17:31 9660 C KELVIN SOURCE MOUSTACAS 21/07/15 01:38:12 5004 SUBROUTINE KELVIN (wrk52,wrk53,wrk54,IB,IGAU,NBPGAU) C C============================================================== C chaine de KELVIN C============================================================== C---------------------------------------------------------------------- C-ENTREES: C------- C-6 = NBR. DE COMPOSANTES DES CONTR. OU DES DEFORM. C-SIG0(6) = CONTR. AU DEBUT DU PAS D'INTEGRATION C-DEPST(6) = INCREMENT DES DEFORM. TOTALES C-NVARI = NBR. DE VARIABLES INTERNES C-VAR0(NVARI) = VARIABLES INTERNES AU DEBUT DU PAS D'INTEGRATION C C-NCOMAT = NBR. DE CARACTERISTIQUES MECANIQUES DU MATERIAU C-IVAL(NCOMAT) = INDICE DES COMPOSANTES DE MATERIAU C-XMAT(NCOMAT) = CARACTERISTIQUES MECANIQUES DU MATERIAU C-MFR = INDICE DE LA FORMULATION MECANIQUE C-ICARA = NBR. DE CARACT. GEOMETRIQUES DES ELEMENTS FINIS C-XCAR(ICARA) = CARACT. GEOMETRIQUES DES ELEMENTS FINIS C-IFOUR = DIMENSION DU MODÈLE 3D = 2 C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCREEL -INC DECHE PARAMETER (UNDEMI=.5D0,UN=1.D0,DEUX=2.D0,TROIS=3.D0) PARAMETER (QUATRE=4.D0) PARAMETER (HEURE= 24.D0,SECONDE=3600.D0) REAL*8 K,K1,K2,K3,TR1,TR2,TR3,ACN REAL*8 RMF1,RMF2,RMF3 REAL*8 JPF1,JPF2,JPF3 DIMENSION AN(6,6),BN(6,6) DIMENSION SR0(6,6),SR1(6,6),SR2(6,6),SR3(6,6) DIMENSION EF0(6,6),EF1(6,6),EF2(6,6),EF3(6,6) DIMENSION DEFTOT(6) DIMENSION SIGF(6) NVARI = VAR0(/1) c WRITE(*,*) 'NVARI', NVARI C TESTER IFOUR POUR VÉRIFIER SI LE MODELE EST BIEN 3D IF (INPLAS.EQ.174) THEN c WRITE(*,*) 'INPLAS', INPLAS c WRITE(*,*) 'CALCUL FLUAGE CHAINE DE KELVIN' c WRITE(*,*) 'ENTRER TAU EN JOURS' ENDIF c WRITE(*,*) 'MATEPL kelvin.eso ', MATEPL c write(*,*) ' inplas jnppla ',inplas,jnpla IF (IFOUR.NE.2) THEN c WRITE(*,*) 'IFOUR', IFOUR c WRITE(*,*) 'L OPTION DE CALCUL N EST PAS IMPLEMENTE' c WRITE(*,*) 'SIMULATION MASSIF 3D ISOTROPE SEULEMENT' KERRE = 99 RETURN ENDIF C---------------------------------------------------------------------- C CARACTERISTIQUES C---------------------------------------------------------------------- KERRE=0 C list xmax pour verifier la positon de rho et alpha epf1=epstf(1) epf2=epstf(2) epf3=epstf(3) epf4=epstf(4) epf5=epstf(5) epf6=epstf(6) YOUN=XMAT(1) XNU=XMAT(2) RHO=XMAT(3) ALP=XMAT(4) YF1=XMAT(5) TF1=XMAT(6) YF2=XMAT(7) TF2=XMAT(8) YF3=XMAT(9) TF3=XMAT(10) SSXX = SIG0(1) SSYY = SIG0(2) SSZZ = SIG0(3) SSXY = SIG0(4) SSXZ = SIG0(5) SSYZ = SIG0(6) TINI = TEMP0 TIFI = TEMPF DELT = (TEMPF-TEMP0) EPSFL01 = VAR0(1) EPSFL02 = VAR0(2) EPSFL03 = VAR0(3) EPSFL04 = VAR0(4) EPSFL05 = VAR0(5) EPSFL06 = VAR0(6) SIGRE01 = VAR0(7) SIGRE02 = VAR0(8) SIGRE03 = VAR0(9) SIGRE04 = VAR0(10) SIGRE05 = VAR0(11) SIGRE06 = VAR0(12) ELAS01 = VAR0(13) ELAS02 = VAR0(14) ELAS03 = VAR0(15) ELAS04 = VAR0(16) ELAS05 = VAR0(17) ELAS06 = VAR0(18) C-MODULE D YOUNG A LA FIN DE LA RELAXATION = C-SOMME DE TOUS LES RESSORT EN PARALLELE K = UN/((UN/YOUN)+(UN/YF1)+(UN/YF2)+(UN/YF3)) C-MODULE D YOUNG RESSORT 1 K1 = (YOUN*YF1)/(YOUN+YF1) K2 = (YOUN*YF2)/(YOUN+YF2) K3 = (YOUN*YF3)/(YOUN+YF3) C -(TEMPF)AU STEME PARALLELE 1 TR1 = (YF1*TF1)/(YOUN+YF1) TR2 = (YF2*TF2)/(YOUN+YF2) TR3 = (YF3*TF3)/(YOUN+YF3) ACN = UN/((UN+XNU)*(UN-(DEUX*XNU))) C-CALCUL TRIDIMENSIONEL C-CREATION DE LA MATRICE DE RIGIDITE ELASTIQUE DO I=1,6 DO J=1,6 AN(I,J)=XZERO BN(I,J)=XZERO SR0(I,J)=XZERO SR1(I,J)=XZERO SR2(I,J)=XZERO SR3(I,J)=XZERO EF0(I,J)=XZERO EF1(I,J)=XZERO EF2(I,J)=XZERO EF3(I,J)=XZERO ENDDO ENDDO IF(IFOUR.EQ.2) THEN AN(1,1)=1.D0 AN(1,2)=-XNU AN(1,3)=-XNU AN(2,1)=-XNU AN(2,2)=1.D0 AN(2,3)=-XNU AN(3,1)=-XNU AN(3,2)=-XNU AN(3,3)=1.D0 AN(4,4)=1.D0+XNU AN(5,5)=1.D0+XNU AN(6,6)=1.D0+XNU BN(1,1)=1.D0-XNU BN(1,2)=XNU BN(1,3)=XNU BN(2,1)=XNU BN(2,2)=1.D0-XNU BN(2,3)=XNU BN(3,1)=XNU BN(3,2)=XNU BN(3,3)=1.D0-XNU BN(4,4)=(1.D0-(2.D0*XNU))/2.D0 BN(5,5)=(1.D0-(2.D0*XNU))/2.D0 BN(6,6)=(1.D0-(2.D0*XNU))/2.D0 C MATRICE RIGIDITÉ SANS YOUNG ENDIF DO I=1,6 DO J=1,6 AN(I,J) = AN(I,J) BN(I,J) = ACN*BN(I,J) ENDDO ENDDO C===========RELAXATION : CALCUL DES CONTRAINTES======================= C===================================================================== RMF1 = (YOUN-K1)*EXP(-(TEMPF/TR1)) RMF2 = (YOUN-K2)*EXP(-(TEMPF/TR2)) RMF3 = (YOUN-K3)*EXP(-(TEMPF/TR3)) DO I=1,6 DO J=1,6 SR0(I,J) = BN(I,J)*K SR1(I,J) = BN(I,J)*RMF1 SR2(I,J) = BN(I,J)*RMF2 SR3(I,J) = BN(I,J)*RMF3 ENDDO ENDDO C======CAS ISOTROPE 3D C======RESSORT ELASTICITE INSTANTANNÉE S0XX = SR0(1,1)+SR0(2,1)+SR0(3,1) C WRITE(*,*) 'SOXX', S0XX S0YY = SR0(1,2)+SR0(2,2)+SR0(3,2) S0ZZ = SR0(1,3)+SR0(2,3)+SR0(3,3) S0XY = SR0(4,4) S0XZ = SR0(5,5) S0YZ = SR0(6,6) C======SYSTEME KV1 S1XX = SR1(1,1)+SR1(2,1)+SR1(3,1) S1YY = SR1(1,2)+SR1(2,2)+SR1(3,2) S1ZZ = SR1(1,3)+SR1(2,3)+SR1(3,3) S1XY = SR1(4,4) S1XZ = SR1(5,5) S1YZ = SR1(6,6) C======SYSTEME KV2 S2XX = SR2(1,1)+SR2(2,1)+SR2(3,1) S2YY = SR2(1,2)+SR2(2,2)+SR2(3,2) S2ZZ = SR2(1,3)+SR2(2,3)+SR2(3,3) S2XY = SR2(4,4) S2XZ = SR2(5,5) S2YZ = SR2(6,6) C======SYSTEME KV3 S3XX = SR3(1,1)+SR3(2,1)+SR3(3,1) S3YY = SR3(1,2)+SR3(2,2)+SR3(3,2) S3ZZ = SR3(1,3)+SR3(2,3)+SR3(3,3) S3XY = SR3(4,4) S3XZ = SR3(5,5) S3YZ = SR3(6,6) C======TOTAL DES CONTRAINTES SIGREF1 = (S0XX+S1XX+S2XX+S3XX)*epf1 SIGREF2 = (S0YY+S1YY+S2YY+S3YY)*epf2 SIGREF3 = (S0ZZ+S1ZZ+S2ZZ+S3ZZ)*epf3 SIGREF4 = (S0XY+S1XY+S2XY+S3XY)*epf4 SIGREF5 = (S0XZ+S1XZ+S2XZ+S3XZ)*epf5 SIGREF6 = (S0YZ+S1YZ+S2YZ+S3YZ)*epf6 SIGF(1) = SIGREF1 SIGF(2) = SIGREF2 SIGF(3) = SIGREF3 SIGF(4) = SIGREF4 SIGF(5) = SIGREF5 SIGF(6) = SIGREF6 c DO I=1,6 c RELL=SIGF(I) c WRITE(6,*) ' I= ',I, 'SIGF = ',SIGF(I), 'RELL', RELL c ENDDO C=============FLUAGE : CALCUL DES DÉFORMATIONS========================= C====================================================================== JPF1 = (UN/YF1)*(UN-EXP(-(TEMPF/TF1))) JPF2 = (UN/YF2)*(UN-EXP(-(TEMPF/TF2))) JPF3 = (UN/YF3)*(UN-EXP(-(TEMPF/TF3))) DO I=1,6 DO J=1,6 EF0(I,J)=AN(I,J)/YOUN EF1(I,J)=AN(I,J)*JPF1 EF2(I,J)=AN(I,J)*JPF2 EF3(I,J)=AN(I,J)*JPF3 ENDDO ENDDO C======RESSORT ELASTIQUE E0XX = EF0(1,1)+EF0(2,1)+EF0(3,1) E0YY = EF0(1,2)+EF0(2,2)+EF0(3,2) E0ZZ = EF0(1,3)+EF0(2,3)+EF0(3,3) E0XY = EF0(4,4) E0XZ = EF0(5,5) E0YZ = EF0(6,6) C======SYSTEME KV1 E1XX = EF1(1,1)+EF1(2,1)+EF1(3,1) E1YY = EF1(1,2)+EF1(2,2)+EF1(3,2) E1ZZ = EF1(1,3)+EF1(2,3)+EF1(3,3) E1XY = EF1(4,4) E1XZ = EF1(5,5) E1YZ = EF1(6,6) C======SYSTEME KV2 E2XX = EF2(1,1)+EF2(2,1)+EF2(3,1) E2YY = EF2(1,2)+EF2(2,2)+EF2(3,2) E2ZZ = EF2(1,3)+EF2(2,3)+EF2(3,3) E2XY = EF2(4,4) E2XZ = EF2(5,5) E2YZ = EF2(6,6) C======SYSTEME KV3 E3XX = EF3(1,1)+EF3(2,1)+EF3(3,1) E3YY = EF3(1,2)+EF3(2,2)+EF3(3,2) E3ZZ = EF3(1,3)+EF3(2,3)+EF3(3,3) E3XY = EF3(4,4) E3XZ = EF3(5,5) E3YZ = EF3(6,6) C======DÉFORMATIONS DUES AU FLUAGE EPSFLF1 = (E1XX+E2XX+E3XX)*SIGREF1 EPSFLF2 = (E1YY+E2YY+E3YY)*SIGREF2 EPSFLF3 = (E1ZZ+E2ZZ+E3ZZ)*SIGREF3 EPSFLF4 = (E1XY+E2XY+E3XY)*SIGREF4 EPSFLF5 = (E1XZ+E2XZ+E3XZ)*SIGREF5 EPSFLF6 = (E1YZ+E2YZ+E3YZ)*SIGREF6 c======DEFORMATIONS ELASTIQUES ELAS1 = E0XX*SIGREF1 ELAS2 = E0YY*SIGREF2 ELAS3 = E0ZZ*SIGREF3 ELAS4 = E0XY*SIGREF4 ELAS5 = E0XZ*SIGREF5 ELAS6 = E0YZ*SIGREF6 C======TOTAL DES DÉFORMATION : DEFORMATIONS DUES AU FLUAGE +DEFORMATIONS C======DUES A L'ELASTICITÉ INSTANTANNÉE DU RESSORT DEF1 = EPSFLF1+ ELAS1 DEF2 = EPSFLF2+ ELAS2 DEF3 = EPSFLF3+ ELAS3 DEF4 = EPSFLF4+ ELAS4 DEF5 = EPSFLF5+ ELAS5 DEF6 = EPSFLF6+ ELAS6 DEFTOT(1) = DEF1 DEFTOT(2) = DEF2 DEFTOT(3) = DEF3 DEFTOT(4) = DEF4 DEFTOT(5) = DEF5 DEFTOT(6) = DEF6 c DO I=1,6 c IF(EPSTF(I).NE.0.D0) THEN c RELL=(DEFTOT(I)-EPSTF(I))/EPSTF(I) c ELSE c RELL=(DEFTOT(I)-EPSTF(I)) c ENDIF c WRITE(6,*) ' I= ',I, 'EPSTF = ',EPSTF(I), c & 'DEFTOT = ', DEFTOT(I), 'RELL', RELL c ENDDO C===== VARIABLE INTERNE AU LA FIN DU PAS VARF(1)= EPSFLF1 VARF(2)= EPSFLF2 VARF(3)= EPSFLF3 VARF(4)= EPSFLF4 VARF(5)= EPSFLF5 VARF(6)= EPSFLF6 VARF(7)= SIGREF1 VARF(8)= SIGREF2 VARF(9)= SIGREF3 VARF(10)= SIGREF4 VARF(11)= SIGREF5 VARF(12)= SIGREF6 VARF(13)=ELAS1 VARF(14)=ELAS2 VARF(15)=ELAS3 VARF(16)=ELAS4 VARF(17)=ELAS5 VARF(18)=ELAS6 END