1 : C INOSC2 SOURCE CHAT 05/01/13 00:38:35 5004 2 : CREAD INOSC1 SOURCE A1 GIBI 06/12/89 12:32:19 3 : C INOSC1 SOURCE CHAT 89/06/01 22:06:52 4 : SUBROUTINE INOSC2(IPT,IPG,DFREQ,XSI,RMAX,TMAX, 5 : 1 DIMALI,DUCTIL, 6 : 2 ALFA,BETA,GAMA,TETA,imod,idema,epse) 7 : IMPLICIT INTEGER(I-N) 8 : IMPLICIT REAL*8(A-H,O-Z) 9 : C 10 : C ===================================================== 11 : C = 12 : C CALCUL DU MAXIMUM DE LA REPONSE A UN SIGNAL D'UN = 13 : C OSCILLATEUR NON LINEAIRE AVEC CONDITIONS INITIALES = 14 : C EN VITESSE ET EN DEPLACEMENT NULLES. = 15 : C = 16 : C CREATION : 15/06/90 = 17 : C PROGRAMMEUR : A.PINTO AND P.PEGON [ISPRA(I)] = 18 : C INSPIREE SUR LA SUBROUTINE 'INOSC1' DE 03/06/87 = 19 : c 20 : c I.P 05/98 ajout du modele bilineaire, iterations sur 21 : c la ductilite, changement du modèle Takeda 22 : C = 23 : c imod = 1 : modele takeda 24 : c 2 elastoplastique ecrouissage isotrope 25 : c 3 elastoplastique ecrouissage cinematique 26 : c 4 elastique nl 27 : C ===================================================== 28 : C 29 : -INC CCOPTIO 30 : -INC SMLREEL 31 : -INC CCREEL 32 : 33 : * quelques tableaux pour les variables internes du modèle Takeda 34 : REAL*8 VAR0(21),VARF(21) 35 : * qques tableaux mour utiliser le meme s_p que DYNE (dychec) 36 : REAL*8 XABSCI(1,4),XORDON(1,4) 37 : 38 : 39 : 40 : MLREE1=IPT 41 : SEGACT MLREE1 42 : MLREE2=IPG 43 : SEGACT MLREE2 44 : NN=MLREE1.PROG(/1) 45 : 46 : C 47 : C DIFFERENCE ENTRE LES 2 PREMIERS TEMPS SUCCESSIFS : DLTF 48 : C (PAS DE TEMP CONSTANT) 49 : C 50 : DLTF=MLREE1.PROG(2)-MLREE1.PROG(1) 51 : ** precision 52 : epsi = .05d0 53 : prec = 1d-3 54 : iprec = 0 55 : nitmax = 30 56 : 57 : C 58 : C CALCUL DE CONSTANTES UTILES 59 : C 60 : W=2*XPI*DFREQ 61 : W2=W*W 62 : 63 : C 64 : C CALCULO DOS COEFICIENTES DE IWAN 65 : C ================================ 66 : DW=XSI*W 67 : D2=XSI*XSI 68 : A0=EXP(-DW*DLTF) 69 : A1=W*SQRT(1.0d0-D2) 70 : AD1=A1*DLTF 71 : A2=SIN(AD1) 72 : A3=COS(AD1) 73 : A4=(2.0d0*D2-1.0d0)/W2 74 : A5=XSI/W 75 : A6=2.0d0*A5/W2 76 : A7=1.0d0/W2 77 : A8=(A1*A3-DW*A2)*A0 78 : A9=-(A1*A2+DW*A3)*A0 79 : A10=A8/A1 80 : A30=A0/A1 81 : A12=A30*A2 82 : A13=A0*A3 83 : A14=A10*A4 84 : A15=A12*A4 85 : A16=A6*A13 86 : A17=A9*A6 87 : A11=A0*(DW*A2/A1+A3) 88 : A21=A10*DW+A9 89 : A22=A10 90 : B11=(-A15-A16+A6)/DLTF-A12*A5-A7*A13 91 : B12=(A15+A16-A6)/DLTF+A7 92 : B21=(-A14-A17-A7)/DLTF-A10*A5-A9*A7 93 : B22=(A14+A17+A7)/DLTF 94 : 95 : ELAS0=W2 96 : ELAS=ELAS0 97 : ELAS1=ELAS0 98 : if (imod.eq.1 .and. teta.le.0d0) teta = .0001d0 99 : ELAS2=ELAS1*TETA 100 : 101 : * qques initialisations pour etre compatible avec dychec 102 : xabsci(1,1) = 0.d0 103 : xordon(1,1) = 0.d0 104 : xjeu = 0d0 105 : iperm = imod 106 : if (imod.eq.4) iperm = -2 107 : nl = 1 108 : xamo = 0d0 109 : iannul = 0 110 : xvit = 0d0 111 : nip = 4 112 : 113 : 400 continue 114 : 115 : iterm = 0 116 : iterm2 = 0 117 : ibisec = 0 118 : 200 continue 119 : iterm = iterm + 1 120 : iterm2 = iterm2 + 1 121 : C 122 : C INICIALIZACAO DE VARIAVEIS 123 : C ========================== 124 : F00=0.D0 125 : F11=0.D0 126 : D0=0.D0 127 : V0=0.D0 128 : FNMDV1=0.D0 129 : FNMDV=0.D0 130 : DES0=0.D0 131 : FORC0=0.D0 132 : RMAX=0.D0 133 : epse= 0.d0 134 : DY=DIMALI/DUCTIL 135 : fy = dy*elas0 136 : ITER=0 137 : * qques parametres de Takeda 138 : dyn = -dy 139 : fyn = -fy 140 : fc = .999d0*fy 141 : fcn = .999d0*fyn 142 : pinp = beta*fy 143 : pinn = beta*fyn 144 : bp = gama/fy 145 : * qques parametres encore pour dychec 146 : if (imod.ne.1) then 147 : xabsci(1,2) = dy 148 : xordon(1,2) = fy 149 : if (elas2.ge.0d0) then 150 : xabsci(1,3) = dy + dy 151 : xordon(1,3) = fy + elas2*dy 152 : xabsci(1,4) = 3d0*dy 153 : xordon(1,4) = fy + elas2*dy*2d0 154 : else 155 : * si raideur negative on considere un plateau a .01fy 156 : xabsci(1,3) = -.9d0*fy/elas2 + dy 157 : xordon(1,3) = .1d0*fy 158 : xabsci(1,4) = 2d0*xabsci(1,3) 159 : xordon(1,4) = .01d0*fy 160 : endif 161 : xdplas = 0d0 162 : xdplac = 0d0 163 : endif 164 : C 165 : C CREATION DE LA BOUCLE 166 : C 167 : DO 100 I=1,NN 168 : C 169 : C INICIO DAS INTEGRACOES 170 : C ====================== 171 : * force exterieure a la fin du pas precedent et courant 172 : F00=F11 173 : F11=-MLREE2.PROG(I) 174 : C INICIO DA CONVERGENCIA 175 : C ====================== 176 : 999 ITER=ITER+1 177 : * on ajoute les forces nl au 2nd membre 178 : F0=F00-FNMDV 179 : F1=F11-FNMDV1 180 : D1=A11*D0+A12*V0+B11*F0+B12*F1 181 : V1=A21*D0+A22*V0+B21*F0+B22*F1 182 : 183 : C 184 : C CALCULO DOS TERMOS NAO-LINEARES NO PASSO 185 : C ======================================== 186 : DELTAD = D1 - D0 187 : 188 : if (imod.ne.1) then 189 : 189 : call DYCHEC(d1,XDPLAS,XDPLAC,XJEU,IPERM,XABSCI, 190 : & XORDON,nl,xfla,nl,NIP,XVIT,XAMO,iannul) 191 : forca = -xfla 192 : 193 : else 194 : * modele Takeda 195 : ** quelques initialisations 196 : forca = forc0 197 : if (iter.eq.1) then 198 : if (i.eq.1) then 199 : do 1 ii = 1,21 200 : var0(ii) = 0.d0 201 : 1 continue 202 : else 203 : do 2 m = 1,21 204 : var0(m) = VARF(m) 205 : 2 continue 206 : endif 207 : endif 208 : 209 : DO 3 k=1,21 210 : VARF(k)=VAR0(k) 211 : 3 CONTINUE 212 : IFR1=INT(VARF(1)) 213 : 214 : IF(IFR1.EQ.0) THEN 215 : IFR1 = 1 216 : VARF(2)=elas0 217 : ENDIF 218 : DP1 = DELTAD * VARF(2) 219 : 219 : CALL DDNSTH(IFR1,IFC1,DP1,FORCA,VARF(21), 220 : * VARF(2),elas0,elas2,elas2,fc,fcn,fy,fyn, 221 : * dy,dyn,alfa,alfa,pinp,pinn,bp,bp, 222 : * VARF(3),VARF(4),VARF(5), 223 : * VARF(6),VARF(7),VARF(8),VARF(9), 224 : * VARF(10),VARF(11),VARF(12),VARF(13), 225 : * VARF(14),VARF(15),VARF(16),VARF(17), 226 : * VARF(18),VARF(19),VARF(20)) 227 : 228 : VARF(1)=DBLE(IFR1) 229 : 230 : endif 231 : 232 : C ANALISE DA CONVERGENCIA EM FORCA 233 : C ================================ 234 : * forces non lineaires à la fin de l'iteration courante 235 : FICT=-D1*ELAS0+FORCA 236 : ERRO=ABS((d1-des0)/deltad) 237 : FNMDV1=FICT 238 : DES0=D1 239 : 240 : C 241 : C ANALISE DE CONVERGENCIA GLOBAL 242 : C ============================== 243 : IF(ERRO.gt.prec .and. ITER.lt.200) goto 999 244 : 245 : * valeurs a la fin du pas 246 : ITER=0 247 : depse = abs(deltad - (forca - forc0)/elas0) 248 : epse = epse + depse 249 : FNMDV=FNMDV1 250 : D0=D1 251 : V0=V1 252 : forc0 = forca 253 : 254 : IF(ABS(D1).GT.RMAX)THEN 255 : RMAX=ABS(D1) 256 : TMAX=MLREE1.PROG(I) 257 : ENDIF 258 : 259 : IF(ITER.GT.200)THEN 260 : print*,'SPON: temps :',MLREE1.PROG(I),'nb d iterations > 200' 261 : ENDIF 262 : C 263 : 100 CONTINUE 264 : 265 : if (idema.eq.0) goto 300 266 : 267 : duct2 = rmax/dy 268 : * 269 : ** petite correction de la limite elastique si elle est trop petite 270 : ** (ca arrive pour des frequences elevees) 271 : ** pour une frequence tres grande ca revient a commencer les iterations 272 : ** avec une ductilite de 1./.9 273 : 274 : c if (iterm2.eq.1) print*,'dfreq',dfreq 275 : c print*,'duct2=',duct2,'iterm=',iterm,'iterm2=',iterm2, 276 : c & 'freq=',dfreq 277 : 278 : c print*,dimali,duct2 279 : 280 : 281 : if (duct2.gt. (5d0*ductil).and.iterm.eq.1 282 : & .and. (abs(teta)).lt.1d0) then 283 : if (iterm2.eq.1 .or. 284 : & (duct2.gt. (6d0*ductil).and.iterm2.le.3)) then 285 : iterm = iterm - 1 286 : dimali = dimali*(.9d0*ductil - 287 : & ((.9d0*ductil - 1)* exp(-.02*(duct2-5d0*ductil)))) 288 : goto 200 289 : endif 290 : endif 291 : 292 : diaux1 = dimali 293 : 294 : if ((abs(duct2 - ductil)).gt. 295 : & (epsi*ductil).and.iterm2.le.nitmax) then 296 : if (iterm.eq.1) dimali = rmax 297 : if (iterm.ge.2) then 298 : aa = (duct2 - ductil)*(duct1-ductil) 299 : if (aa.gt.0.d0. and. ibisec.eq.0) then 300 : dimali = rmax 301 : else if (ibisec.eq.0) then 302 : * premiere dichotomie 303 : ibisec = 1 304 : dimali = (dimali + dimal1)/2d0 305 : ducaux = duct1 306 : diaux2 = dimal1 307 : else 308 : * dichotomies suivantes 309 : aa2 = (duct2 - ductil)*(duct3-ductil) 310 : if (aa.lt.0.d0) then 311 : dimali = (dimali + dimal1)/2d0 312 : ducaux = duct1 313 : diaux2 = dimal1 314 : else if (aa2.lt.0.d0) then 315 : dimali = (dimali + dimal3)/2d0 316 : ducaux = duct3 317 : diaux2 = dimal3 318 : else 319 : * dans le cas ou la dichotomie ne marche pas 320 : dimali = rmax 321 : endif 322 : endif 323 : endif 324 : duct1 = duct2 325 : dimal1 = diaux1 326 : if (ibisec.eq.1) then 327 : duct3 = ducaux 328 : dimal3 = diaux2 329 : endif 330 : goto200 331 : else if (iterm2.gt.nitmax) then 332 : * on essaie 2 fois encore avec un critere plus petit 333 : if (iprec.le.1) then 334 : prec = prec/10d0 335 : iprec = iprec + 1 336 : goto 400 337 : endif 338 : REAERR(1)=dfreq 339 : REAERR(2)=duct2 339 : CALL ERREUR(-326) 340 : endif 341 : rmax = dy 342 : 343 : 300 continue 344 : SEGDES MLREE1 345 : SEGDES MLREE2 346 : 347 : 348 : RETURN 349 : END 350 : 351 : 352 : 353 : 354 :
© Cast3M 2003 - All rights reserved.
Disclaimer