Esope source files

   1 : C LABORD    SOURCE    CHAT      12/04/06    21:15:21     7348           
   2 : C     ================================================================== 
   3 : C     =       SUBROUTINE INTEGRANT LE MODELE D'ENDOMMAGEMENT           = 
   4 : C     =     UNILATERAL EN UNIAXIAL AVEC PILOTAGE EN DEFORMATION        = 
   5 : C     =              FG + FR (LMT-Cachan, fevrier 2001)                = 
   6 : C     ================================================================== 
   7 :       SUBROUTINE LABORD(XMAT,DEPS,SIG0,VAR0,SIGF,VARF) 
   8 : 
   9 :       IMPLICIT INTEGER(I-N)
  10 :       IMPLICIT REAL*8 (a-h,o-z)
  11 : 
  12 :       REAL*8 XMAT(13),DEPS(3),SIG0(3),VAR0(8),SIGF(3),VARF(8)
  13 : 
  14 : C--------------------------------------------------------------- 
  15 : C     MODELE DE CLB (CLB_UNI) SPECIALISE POUR LE MODELE A FIBRE 
  16 : C--------------------------------------------------------------- 
  17 : C 
  18 : C     variables en entree 
  19 : C 
  20 : C     XMAT : Caracteristique du materiaux 
  21 : C 
  22 : C     EPS : Deformation totale 
  23 : C 
  24 : C     DEPS : Increment de deformations 
  25 : C 
  26 : C     SIG0 : Contraintes initiales 
  27 : C 
  28 : C     VAR0 : Variables internes initiales 
  29 : C 
  30 : C 
  31 : C     variables en sortie 
  32 : C 
  33 : C     SIGF : Contraintes finales 
  34 : C 
  35 : C     VARF : Variables internes finales 
  36 : C 
  37 : C     declaration des variables 
  38 : C 
  39 :       logical change 
  40 : 
  41 : C     ================================================================== 
  42 : C     =     Initialisation des Constantes                              = 
  43 : C     ================================================================== 
  44 :       E=XMAT(1) 
  45 :       XNU=XMAT(2) 
  46 :       Y01=XMAT(5) 
  47 :       Y02=XMAT(6) 
  48 :       A1=XMAT(7) 
  49 :       A2=XMAT(8) 
  50 :       B1=XMAT(9) 
  51 :       B2=XMAT(10) 
  52 :       BETA1=XMAT(11) 
  53 :       BETA2=XMAT(12) 
  54 :       SIGF1=XMAT(13) 
  55 : C     ================================================================== 
  56 : C     =     Initialisation des VARIABLES HISTORIQUES                   = 
  57 : C     ================================================================== 
  58 :       ESEC = VAR0(1) 
  59 :       D1 = VAR0(2) 
  60 :       D2 = VAR0(3) 
  61 :       YLIM1 = VAR0(4) 
  62 :       YLIM2 = VAR0(5) 
  63 :       EPS10 = VAR0(6)
  64 : C 
  65 :       EPS = EPS10 + DEPS(1) 
  66 :       y1=YLIM1-1.D0 
  67 :       y2=YLIM2-1.D0 
  68 : C     ================================================================== 
  69 : *** 
  70 :       change=.true. 
  71 :       do while (change) 
  72 : ***       
  73 :       X1=(BETA1*D1/(1.D0-D1)+BETA2*D2/(1.D0-D2))/E 
  74 :       X2=(BETA2*D2-SIGF1)/(1.D0-D2)/E 
  75 :       if ( X1 .LE. Eps ) then 
  76 :         EPSF=Eps-BETA2*D2/(1.0D0-D2)/E  
  76 :         CALL CAS13(EPSF,SIG,ESEC,E,D1,BETA1,A1,B1,Y1,YLIM1,Y01)
  77 :         icas1=1
  78 :       else if ( Eps .LE. X2 ) then 
  78 :         CALL CAS13(Eps,SIG,ESEC,E,D2,BETA2,A2,B2,Y2,YLIM2,Y02) 
  79 :         icas1=3 
  80 :       else  
  80 :         call CAS2(Eps,BETA1,E,D1,SIGF1,BETA2,D2,SIG,ESEC) 
  81 :         icas1=2 
  82 :       endif 
  83 :       X1=(BETA1*D1/(1.D0-D1)+BETA2*D2/(1.D0-D2))/E
  84 :       X2=(BETA2*D2-SIGF1)/(1.D0-D2)/E 
  85 :       if ( X1 .LE. Eps ) then 
  86 :         icas2=1 
  87 :       else if ( Eps .LE. X2 ) then 
  88 :         icas2=3 
  89 :       else  
  90 :         icas2=2 
  91 :       endif 
  92 :       if (icas1 .eq. icas2) then 
  93 :          change=.false. 
  94 :       endif 
  95 : ***       
  96 :       enddo 
  97 : ***       
  98 : c     Si on garde le module du beton constant 
  99 :       ESEC=E 
 100 :  
 101 : C     ================================================================== 
 102 : C     =     reactualisation des seuils                                 = 
 103 : C     ================================================================== 
 104 :       YLIM1=DMAX1(Y1,YLIM1) 
 105 :       YLIM2=DMAX1(Y2,YLIM2) 
 106 : 
 107 : C     ================================================================== 
 108 : C     =     reactualisation des VARIABLES HISTORIQUES                  = 
 109 : C     ================================================================== 
 110 :       VARF(1) = ESEC 
 111 :       VARF(2) = D1 
 112 :       VARF(3) = D2 
 113 :       VARF(4) = YLIM1 
 114 :       VARF(5) = YLIM2 
 115 :       VARF(6) = EPS 
 116 :       VARF(7) = VAR0(7) + DEPS(2) 
 117 :       VARF(8) = VAR0(8) + DEPS(3) 
 118 : C     ================================================================== 
 119 : C     =     Calcul des contraintes finales                             = 
 120 : C     ================================================================== 
 121 :       SIGF(1) = SIG 
 122 :       r_z = 0.5D0 * E / (1.D0+XNU)
 123 :       SIGF(2) = r_z * VARF(7)
 124 :       SIGF(3) = r_z * VARF(8)
 125 :             
 126 :       RETURN 
 127 :       END                                  
 128 :  
 129 : C     ================================================================== 
 130 : C     =               CALCUL DE SIGMA DANS LE CAS 1 & 3                = 
 131 : C     =             ( le calcul se fait par dichotomie )               = 
 132 : C     ================================================================== 
 133 :       SUBROUTINE CAS13(EPS,SIG,ESEC,E,D,BETA,A,B,Y,YLIM,Y00) 
 134 :  
 135 :       IMPLICIT INTEGER(I-N)
 136 :       IMPLICIT REAL*8(A-H,O-Z)
 137 : 
 138 :       REAL*8 m0,m
 139 : 
 140 :       dt = 1.0d0
 141 :       m0 = 0.0d0
 142 :       alpha = 0.0d0
 143 : 
 144 : c     Erreur relative par rapport a D 
 145 :       Erreur_max=1.0D-6 
 146 :       Nbpas=50 
 147 : 
 148 : c     Dans le cas No 1 ou 3, on recherche D13 dans [0,1] 
 149 :       X1= d-0.01D0 
 150 :       X2= 1.0D0  
 151 : 
 152 :       pas=(X2-X1)/Nbpas 
 153 :       m=m0 
 154 :       if ( EPS .ge. 0.D0) then 
 155 :          m=m0/650.D0
 156 :       endif 
 157 : 
 158 : C     ------------------------------------------------------------------ 
 159 : C     -   Avant de faire les calculs, on verifie que l'on depasse le   - 
 160 : C     -   seuil YLIM                                                   - 
 161 : C     ------------------------------------------------------------------ 
 162 :       D13=D 
 163 :       y0=y00 
 164 :       YY = (E*EPS+BETA)*(E*EPS+BETA)-BETA*BETA/(1.0D0-D13)/(1.0D0-D13) 
 165 :       YY = YY/(2.0D0*E)              
 166 : 
 167 : ***   Prise en compte de la viscosite 
 168 :  
 169 : ****  coefficient de viscosite  
 170 : c      m=3.d4 
 171 : c      alpha=0.18d0 
 172 : ****      
 173 : **      Si yy est inferieur a y0 on sort sans evolution de d 
 174 :       If ( YY .LE. Y0 ) Then 
 175 :         y=yy 
 176 : c       Module Secant 
 177 :         ESEC=E*(1.0D0-D) 
 178 : c       Valeur de la contrainte 
 179 :         SIG=EPS*ESEC-BETA*D   
 180 :         return 
 181 :       endif 
 182 :  
 183 : *********** 
 184 : **      Calcul du seuil 
 185 : c      y=(E*EPS+BETA)*(E*EPS+BETA)/(2.D0*E) 
 186 :        
 187 : **      Calcul de l endommagement 
 188 :       Erreur=1.D0 
 189 :       XG=d 
 190 :       Xd=X1 
 191 :       Do while ( (Erreur .GT. Erreur_max) .and. (pas .ge. 1.d-8) )  
 192 :  
 193 : c       Valeur de la fonction a la borne de droite D13=XD 
 194 :         D13=XD 
 195 :         YY = (E*EPS+BETA)*(E*EPS+BETA)/(2.0D0*E) 
 196 :      &      -(BETA*BETA/(1.D0-D13)/(1.D0-D13))/(2.D0*E) 
 197 : c        if (yy .lt. ylim) then 
 198 : c          yy=ylim 
 199 : c        endif 
 200 :         y0=y00+m*( (abs(d13-d)+(d13-d))/(2.D0*dt) )**alpha 
 201 :         If ( YY .GT. Y0 ) Then 
 202 : **      si Dconverge est > a D13 alors fctd est < 0 
 203 :           Fctd =  1.0D0 - (1.0D0+(A*(YY-Y0))**B)*(1.0D0-D13) 
 204 :         Else 
 205 :           Fctd =  1.d0 
 206 :         Endif 
 207 :  
 208 :         If ( Fctd .lt. 0.D0) then 
 209 :           Xg=Xd 
 210 :           Xd=Xd+pas 
 211 :         Else 
 212 :           Xd=Xg 
 213 :           Pas=Pas * 0.5D0
 214 :           Xd=Xd+pas 
 215 :         Endif 
 216 :  
 217 :         If ( Xg .gt. X2 ) Then 
 218 :           Stop ' . DANS CAS13 . Pas de Zero.' 
 219 :         Endif 
 220 :  
 221 :         Erreur=abs(Fctd)  
 222 :       Enddo 
 223 :  
 224 : c     On verifie que d13 > d . Si ce n'est pas le cas rien n'evolue. 
 225 : c      print*, (abs(d13-d)+(d13-d))/(2.D0*dt) 
 226 : *** 
 227 :       y=y00+1.D0/a*(d13/(1.D0-d13))**(1.D0/b) 
 228 : *** 
 229 :       If ( d13 .Ge. d ) then 
 230 :         D = D13 
 231 :       Endif  
 232 :       if (d .gt. 0.9999999D0) then 
 233 :          d=0.9999999D0 
 234 :       endif    
 235 :              
 236 : c     Module Secant 
 237 :       ESEC=E*(1.0D0-D) 
 238 : c     Valeur de la contrainte 
 239 :       SIG=EPS*ESEC-BETA*D 
 240 : 
 241 :       Return
 242 :       End 
 243 :  
 244 : C     ================================================================== 
 245 : C     =                 CALCUL DE SIGMA DANS LE CAS 2                  = 
 246 : C     =             ( le calcul se fait par dichotomie )               = 
 247 : C     ================================================================== 
 248 :       SUBROUTINE CAS2(EPS,BETA1,E,D1,SIGF1,BETA2,D2,SIG,ESEC) 
 249 : 
 250 :       IMPLICIT INTEGER(I-N)
 251 :       IMPLICIT REAL*8(A-H,O-Z)
 252 : 
 253 :       fprim=  (e*eps*(1.D0-d2)-beta2*d2+sigf1)
 254 :      &      / (sigf1+beta1*d1*(1.D0-d2)/(1.D0-d1)) 
 255 : 
 256 : c     La solution est la borne de droite. 
 257 :       SIG=e*eps*(1.D0-d2)-beta2*d2-beta1*d1*(1.D0-d2)*fprim/(1.D0-d1)  
 258 : 
 259 : c     ESEC est le module du beton de compression 
 260 :       ESEC= E*(1.0D0-D2) 
 261 : 
 262 :       Return 
 263 :       End 
 264 : 
 265 : C=======================================================================
 266 : C REMARQUE : LES SOUS-PROGRAMMES RAIDTANi (avec i=1,2,3) NE SONT
 267 : C            PAS UTILISES !
 268 : C=======================================================================
 269 :       SUBROUTINE RAIDTAN1 (EPS,SIG,ETAN,E,D1,BETA1,D2,BETA2,A,B,Y,Y0) 
 270 : 
 271 :       IMPLICIT INTEGER(I-N)
 272 :       IMPLICIT REAL*8(A-H,O-Z)
 273 :  
 274 : **      derive de d par rapport a y 
 275 :         dddy= b*(a**b)*(y-y0)**(b-1.D0)/(1.D0+(a*(y-y0))**b)**2 
 276 : 
 277 : **      derive de d par rapport a epsilon 
 278 :         dddeps=dddy*( e*eps+beta1-beta2*d2/(1.D0-d2) ) 
 279 :         dddeps=dddeps/( 1.D0+dddy*beta1**2/(e*(1.D0-d1)**3) ) 
 280 : 
 281 : **      derive de sigma par rapport a epsilon -- calcul de E tangent 
 282 :         etan=e*(1.D0-d1)-e*eps*dddeps-beta1*dddeps+beta2*d2*dddeps/
 283 :      $  (1.D0-d2) 
 284 : 
 285 :         return  
 286 :         end 
 287 : C---------
 288 :       SUBROUTINE RAIDTAN2 (EPS,BETA1,E,D1,SIGF1,BETA2,D2,SIG,ETAN) 
 289 :  
 290 :       IMPLICIT INTEGER(I-N)
 291 :       IMPLICIT REAL*8(A-H,O-Z)
 292 : 
 293 :         etan=E*(1.D0-d2)/( sigf1+beta1*d1*(1.D0-d2)/(1.D0-d1) ) 
 294 :         etan=etan*beta1*d1*(1.D0-d2)/(1.D0-d1) 
 295 :         etan=e*(1.D0-d2)-etan 
 296 :  
 297 :         return 
 298 :         end 
 299 : C---------
 300 :       SUBROUTINE RAIDTAN3 (EPS,SIG,ETAN,E,D1,BETA1,D2,BETA2,A,B,Y,Y0) 
 301 :  
 302 :       IMPLICIT INTEGER(I-N)
 303 :       IMPLICIT REAL*8(A-H,O-Z)
 304 :  
 305 : **      derive de d par rapport a y 
 306 :         dddy= b*(a**b)*(y-y0)**(b-1.D0)/(1.D0+(a*(y-y0))**b)**2 
 307 :  
 308 : **      derive de d par rapport a epsilon 
 309 :         dddeps=dddy*(e*eps+beta2) 
 310 :         dddeps=dddeps/(1.D0+dddy*beta2**2/(e*(1.D0-d2)**3)) 
 311 :  
 312 : **      derive de sigma par rapport a epsilon -- calcul de E tangent 
 313 :         etan=e*(1.D0-d2)-e*eps*dddeps-beta2*dddeps 
 314 :  
 315 :         return  
 316 :         end 
 317 :  
 318 : C     ================================================================== 
 319 : C     =                 CALCUL DE SIGMA DANS LE CAS 2                  = 
 320 : C     =             ( le calcul se fait par dichotomie )               = 
 321 : C     ================================================================== 
 322 :       SUBROUTINE CAS21(EPS,BETA1,E,D1,SIGF1,BETA2,D2,SIG,ESEC) 
 323 :  
 324 :       IMPLICIT INTEGER(I-N)
 325 :       IMPLICIT REAL*8(A-H,O-Z)
 326 : 
 327 : c     Erreur relative par rapport a Sigma 
 328 :       Erreur_max=1.0D-6 
 329 :       Nbpas=50 
 330 :       itermax=1000 
 331 :       iter=0 
 332 : c     Dans le cas No 2, on recherche Y dans [-1,0]   Y=sig/sigf1 
 333 :       X1= -1.0D0
 334 :       X2=  0.0D0
 335 :  
 336 :       pas=(X2-X1)/Nbpas
 337 :  
 338 : c     Valeur initiale des bornes de gauche et droite 
 339 :       XG=X1 
 340 :       XD=XG+PAS 
 341 :  
 342 : c     Valeur de la fonction a la borne de gauche y=XG 
 343 :       Y=XG 
 344 :       Fctg    =  EPS*E - BETA1*D1*(1.0D0-3.0D0*Y*Y-2.0D0*Y*Y*Y)/
 345 :      $  (1.0D0-D1) 
 346 :      &        - BETA2*D2/(1.0D0-D2) - Y*SIGF1/(1.0D0-D1) 
 347 :      &        - SIGF1*(1.0D0/(1.0D0-D1)-1.0D0/(1.0D0-D2))*(2.0D0+Y)*Y*Y 
 348 :       If (Fctg .GE. 0.0D0) then 
 349 :         isigg=1 
 350 :       else 
 351 :         isigg=0 
 352 :       endif 
 353 :       Erreur=abs(Fctg) 
 354 : C
 355 :       Do while ( (Erreur .GT. Erreur_max) .and. (iter .le. itermax) )  
 356 :          iter=iter+1 
 357 : c       Valeur de la fonction a la borne de droite y=XD 
 358 :         Y=XD 
 359 :         Fctd    =  EPS*E - BETA1*D1*(1.0D0-3.0D0*Y*Y-2.0D0*Y*Y*Y)/ 
 360 :      &   (1.0D0-D1)    - BETA2*D2/(1.0D0-D2) - Y*SIGF1/(1.0D0-D1) 
 361 :      &       - SIGF1*(1.0D0/(1.0D0-D1)-1.0D0/(1.0D0-D2))*(2.0D0+Y)*Y*Y 
 362 :         If (Fctd .GE. 0.0D0) then 
 363 :           isigd=1 
 364 :         else 
 365 :           isigd=0 
 366 :         endif 
 367 :         Erreur=abs(Fctd) 
 368 :         If ( isigg .EQ. isigd) then 
 369 :           Fctg=Fctd 
 370 :           isigg=isigd
 371 :           Xg=Xd 
 372 :           Xd=Xd+pas 
 373 :         Else 
 374 :           Pas=Pas * 0.5D0
 375 :           Xd=Xg+pas 
 376 :         Endif 
 377 :         If ( Xg .gt. X2 ) Then 
 378 :           Stop ' . DANS CAS2 . Pas de Zero.' 
 379 :         Endif 
 380 :       Enddo 
 381 : c     La solution est la borne de droite. 
 382 :       SIG=Y*Sigf1 
 383 : c     Si on fait varier lineairement ESEC en fonction de Sig 
 384 :       ESEC= E*(1.0D0 - D1 + Y*(D2-D1)) 
 385 :       Return 
 386 :       End 
 387 : 
 388 :  
 389 :