1 : C ISSGRA SOURCE BP208322 13/11/13 21:15:08 7865 2 : C 3 : SUBROUTINE ISSGRA(WRK52,WRK53,WRK54,WRK27,IB,IGAU,NBPGAU) 4 : C 5 : C______________________________________________________________________ 6 : C 7 : C MODELE ISSGRA 8 : C______________________________________________________________________ 9 : C 10 : IMPLICIT INTEGER(I-N) 11 : IMPLICIT REAL*8(A-H,O-Z) 12 : c inutile? -INC CCOPTIO 13 : -INC CCREEL 14 : -INC DECHE 15 : * 16 : C GLOBAL MODEL 17 : DIMENSION XMKel(5,5) 18 : C DIMENSION XMCel(5,5) 19 : DIMENSION uo(5), duo(5), Fo(5) 20 : DIMENSION XMK(5,5), Ftr(5), force(5) 21 : * 22 : DIMENSION upl(5), uel(5) 23 : DIMENSION uoi(5), duoi(5), dueli(5), force_preci(5) 24 : C DIMENSION ueli(5) 25 : 26 : C Modèle de Plasticité 27 : DIMENSION ddlambda(2), yield(2) 28 : DIMENSION Ho(2,2), H(2,2), Hinv(2,2) 29 : DIMENSION duP(5) 30 : * 31 : C Modèle de plasticité 1 32 : DIMENSION q1(5), upl1(5) 33 : C DIMENSION q1_prec(5) 34 : DIMENSION Gradgs(5),Gradfs(5),Gradfq(5) 35 : DIMENSION hh1(5), XKGradgs(5) 36 : * 37 : C Modèle de plasticité vertical 38 : DIMENSION upl2(5), Gradf2s(5), XKGradf2s(5) 39 : * 40 : C Modèle de Décollement 41 : DIMENSION dforce(5), forcepred(5), force_prec(5), duel(5) 42 : * 43 : C Amortissment 44 : DIMENSION upl_preci(5),q1_preci(5) 45 : * 46 : DIMENSION da(20) 47 : C DIMENSION TIME(2) 48 : * 49 : C----------------------------------------------------------------------- 50 : C----------------------------------------------------------------------- 51 : C Affectation des données du modèle 52 : da(1) =XMAT(13) 53 : da(2) =XMAT(14) 54 : da(3) =XMAT(15) 55 : da(4) =XMAT(7) 56 : da(5) =XMAT(8) 57 : da(6) =XMAT(12) 58 : da(7) =XMAT(11) 59 : da(20)=XMAT(10) 60 : da(8) =XMAT(16) 61 : da(9) =XMAT(17) 62 : da(10)=XMAT(18) 63 : da(11)=XMAT(19) 64 : da(12)=XMAT(20) 65 : da(13)=XMAT(21) 66 : da(14)=XMAT(22) 67 : da(15)=XMAT(23) 68 : da(16)=XMAT(24) 69 : da(17)=XMAT(25) 70 : da(18)=XMAT(26) 71 : da(19)=XMAT(27) 72 : * 73 : Diam=da(1) 74 : XLx=da(2) 75 : XLy=da(3) 76 : XKelz=da(4) 77 : XKelh=da(5) 78 : XKelry=da(6) 79 : XKelrx=da(7) 80 : C a=da(8) 81 : C b=da(9) 82 : C c=da(10) 83 : C d=da(11) 84 : C e=da(12) 85 : C f=da(13) 86 : qmax=da(14) 87 : a9=da(15) 88 : C a6=da(16) 89 : eta3=da(17) 90 : dtime=da(18) 91 : a8=da(19) 92 : XKelt=da(20) 93 : * 94 : * 95 : upl1(1)=VAR0(2) 96 : upl1(2)=VAR0(3) 97 : upl1(3)=VAR0(4) 98 : upl1(4)=VAR0(5) 99 : upl1(5)=VAR0(6) 100 : q1(1)=VAR0(7) 101 : q1(2)=VAR0(8) 102 : q1(3)=VAR0(9) 103 : q1(4)=VAR0(10) 104 : q1(5)=VAR0(11) 105 : force_prec(1)=-SIG0(1) 106 : force_prec(2)=-SIG0(2) 107 : force_prec(3)=SIG0(6) 108 : force_prec(4)=SIG0(3) 109 : force_prec(5)=-SIG0(5) 110 : delta_prec=VAR0(12) 111 : Xkapa_prec=VAR0(13) 112 : Xksi_prec=VAR0(14) 113 : distmin_prec=VAR0(15) 114 : upl2(1)=VAR0(17) 115 : upl2(2)=VAR0(18) 116 : upl2(3)=VAR0(19) 117 : upl2(4)=VAR0(20) 118 : upl2(5)=VAR0(21) 119 : q2=VAR0(22) 120 : DO I=1,5 121 : upl(I)=upl1(I)+upl2(I) 122 : ENDDO 123 : * 124 : *deplx joint --> -deplz ISS 125 : uo(1)=-epst0(1) 126 : *deply joint --> -deplx ISS 127 : uo(2)=-epst0(2) 128 : *rotz joint --> -roty ISS 129 : uo(3)=-epst0(6) 130 : *deplz joint --> deply ISS 131 : uo(4)=epst0(3) 132 : *roty joint --> -rotx ISS 133 : uo(5)=-epst0(5) 134 : * 135 : duo(1)=-DEPST(1) 136 : duo(2)=-DEPST(2) 137 : duo(3)=-DEPST(6) 138 : duo(4)=DEPST(3) 139 : duo(5)=-DEPST(5) 140 : 141 : C Correction données d'entrées Castem 142 : DO I=1,5 143 : uo(I)=uo(I)+duo(I) 144 : * epstf(I) = epst0(I)+DEPST(I) 145 : ENDDO 146 : * 147 : C on adimensionne les déplacements 148 : if (a9.eq.2) then 149 : DO I=1,5 150 : uo(I) = uo(I)*sqrt(XLx**2+XLy**2)/(XLx*XLy) 151 : ENDDO 152 : uo(3) = uo(3)*XLx 153 : uo(5) = uo(5)*XLy 154 : DO I=1,5 155 : duo(I) = duo(I)*sqrt(XLx**2+XLy**2)/(XLx*XLy) 156 : ENDDO 157 : duo(3)= duo(3)*XLx 158 : duo(5)= duo(5)*XLy 159 : else 160 : DO I=1,5 161 : uo(I) = uo(I)/Diam 162 : ENDDO 163 : uo(3) = uo(3)*Diam 164 : uo(5) = uo(5)*Diam 165 : DO I=1,5 166 : duo(I) = duo(I)/Diam 167 : ENDDO 168 : duo(3)= duo(3)*Diam 169 : duo(5)= duo(5)*Diam 170 : endif 171 : * 172 : C Adimensionnement des efforts 173 : if (a9.eq.1) then 174 : C semelle filante 175 : DO I=1,5 176 : force_prec(I)=force_prec(I)/(Diam*qmax) 177 : ENDDO 178 : force_prec(3)=force_prec(3)/Diam 179 : force_prec(5)=force_prec(5)/Diam 180 : else 181 : if (a9.eq.2) then 182 : C semelle rectangulaire 183 : DO I=1,5 184 : force_prec(I)=force_prec(I)/(XLx*XLy*qmax) 185 : ENDDO 186 : force_prec(3)=force_prec(3)/XLx 187 : force_prec(5)=force_prec(5)/XLy 188 : else 189 : C semelle circulaire 190 : DO I=1,5 191 : force_prec(I)=force_prec(I)/(XPI*Diam**2/4*qmax) 192 : ENDDO 193 : force_prec(3)=force_prec(3)/Diam 194 : force_prec(5)=force_prec(5)/Diam 195 : endif 196 : endif 197 : C *********************************************************************** 198 : C Macro-élément d'ISS 199 : C adimensionnement de la matrice de rigidité élastique suivant la forme de la fondation (paramètre a9) 199 : call MKADIM(da,XMKel) 200 : do i=1,5 201 : do j=1,5 202 : XMK(i,j)=XMKel(i,j) 203 : enddo 204 : enddo 205 : 206 : C Calcul de l'effort de prédiction élastique en considérant le décollement élastique non-linéaire. 207 : duP(1)=0 208 : duP(2)=0 209 : duP(3)=0 210 : duP(4)=0 211 : duP(5)=0 212 : 213 : C Réduction du pas de calcul 214 : dtimei=dtime/10. 215 : DO I=1,5 216 : duoi(I)=1./10.*duo(I) 217 : force_preci(I)=force_prec(I) 218 : ENDDO 219 : delta_preci=delta_prec 220 : Xkapa_preci=Xkapa_prec 221 : Xksi_preci=Xksi_prec 222 : distmin_preci=distmin_prec 223 : do k=1,10 224 : q2_preci=q2 225 : * 226 : Xkapa_k=Xkapa_preci 227 : Xksi_k=Xksi_preci 228 : DO I=1,5 229 : q1_preci(I)=q1(I) 230 : upl_preci(I)=upl(I) 231 : uoi(I)=uo(I)-(10.-k)*duoi(I) 232 : uel(I)=uoi(I)-upl(I) 233 : duel(I)=duoi(I)-duP(I) 234 : ENDDO 235 : *si a8=1 on désactive le décollement 236 : if (a8.eq.1) then 236 : call MATVE1 (XMK, uel,5,5,Ftr,2) 237 : delta=delta_prec 238 : deltamax=1 239 : else 239 : call MATVE1 (XMK, duel,5,5,dforce,2) 240 : DO I=1,5 241 : force(I)=force_preci(I)+dforce(I) 242 : ENDDO 242 : call ELNLIN (da, uel, force, XMKel, 243 : & delta_preci, Ftr, delta, deltamax) 244 : endif 245 : C Calcul des fonctions de charge au point de charge élastique 245 : call LOADSF (da, Ftr,q1,yield1) 245 : call LOASFV (Ftr,q2,yield2) 246 : DO I=1,5 247 : force(I)=Ftr(I) 248 : ENDDO 249 : C Boucle de plasticité (Return Mapping) 250 : epsilon=1E-9 251 : if ((yield1.le.epsilon) .and. (yield2.le.epsilon)) then 252 : DO I=1,5 253 : force(I)=Ftr(I) 254 : ENDDO 255 : Xkapa=Xkapa_k 256 : Xksi=Xksi_k 257 : distmin=distmin_preci 258 : else 259 : m=0 260 : do while ((yield1.gt.epsilon) .or. (yield2.gt.epsilon)) 261 : m=m+1 262 : IF (m.gt.1000) THEN 263 : KERRE = 22 264 : RETURN 265 : ENDIF 266 : IF ((force(1).LT.0).or.(force(1).GE.1)) THEN 267 : KERRE = 23 268 : RETURN 269 : ENDIF 269 : call GRADIF(da,force,q1,Gradfs,Gradfq) 269 : call GRADIG(da,force, q1, force_preci, Xkapa_k, 270 : & Xksi_k, Xkapa_preci, Xksi_preci, Gradgs, Xkapa, Xksi) 270 : call FVEH(da,force,q1, XMKel, 271 : & Gradgs,distmin_preci,hh1,distmin) 271 : call GRADFV(force,q2,Gradf2s,Gradf2q) 271 : call FVEHV(force,q2, XMKel,Gradf2s,hh2) 272 : 272 : call MATVE1(XMK,Gradgs,5,5,XKGradgs,2) 272 : call MATVE1(XMK,Gradf2s,5,5,XKGradf2s,2) 273 : 273 : call SCALPR(5,Gradfs,XKGradgs,Ho11) 273 : call SCALPR(5,Gradfs,XKGradf2s,Ho12) 273 : call SCALPR(5,Gradf2s,XKGradgs,Ho21) 273 : call SCALPR(5,Gradf2s,XKGradf2s,Ho22) 274 : 274 : call SCALPR(5,Gradfq,hh1,H1) 275 : H2=Gradf2q*hh2 276 : H(1,1)=Ho11+H1 277 : H(1,2)=Ho12 278 : H(2,1)=Ho21 279 : H(2,2)=Ho22+H2 280 : yield(1)=yield1 281 : yield(2)=yield2 282 : *initialisation de ces variables necessaire pour les tests sur les mécanismes 283 : ddlambda(1)=1 284 : ddlambda(2)=1 285 : if ((yield(1).gt.epsilon). 286 : & and.(yield(2).gt.epsilon)) then 286 : call HINV22 (H,Hinv) 286 : call MATVE1(Hinv,yield,2,2,ddlambda,2) 287 : ddlambda1=ddlambda(1) 288 : ddlambda2=ddlambda(2) 289 : endif 290 : C Tests sur les mécanismes 291 : if ((yield(1).le.epsilon).or.(ddlambda(1).le.0)) then 292 : ddlambda1=0 293 : ddlambda2=yield2/H(2,2) 294 : if (ddlambda2.le.0) then 295 : ddlambda2=0 296 : endif 297 : endif 298 : if ((yield(2).le.epsilon).or.(ddlambda(2).le.0)) then 299 : ddlambda2=0 300 : ddlambda1=yield1/H(1,1) 301 : if (ddlambda1.le.0) then 302 : ddlambda1=0 303 : endif 304 : endif 305 : 306 : ddlambda(1)=ddlambda1 307 : ddlambda(2)=ddlambda2 308 : 309 : 310 : q2=q2-ddlambda(2)*hh2 311 : DO I=1,5 312 : q1(I)=q1(I)-ddlambda(1)*hh1(I) 313 : upl1(I)=upl1(I)+ddlambda(1)*Gradgs(I) 314 : upl2(I)=upl2(I)+ddlambda(2)*Gradf2s(I) 315 : upl(I) =upl1(I)+upl2(I) 316 : * 317 : uel(I)=uoi(I)-upl(I) 318 : duP(I)=upl(I)-upl_preci(I) 319 : duel(I)=duoi(I)-duP(I) 320 : ENDDO 321 : *si a8=1 on désactive le décollement 322 : if (a8.eq.1) then 322 : call MATVE1 (XMK, uel,5,5, force,2) 323 : delta=delta_prec 324 : else 324 : call MATVE1 (XMK, duel, 5,5, dforce,2) 325 : DO I=1,5 326 : forcepred(I)=force_preci(I)+dforce(I) 327 : ENDDO 328 : *forcepred et delta_prec sont là pour une initialisation à la boucle de Newton 328 : call ELNLIN (da, uel, forcepred, XMKel, 329 : & delta_preci, force, delta, deltamax) 330 : endif 331 : Xkapa_k=Xkapa 332 : Xksi_k=Xksi 332 : call LOADSF (da,force,q1,yield1) 332 : call LOASFV (force,q2,yield2) 333 : enddo 334 : endif 335 : * 336 : C Correction visqueuse 337 : * calcul de la norme de la force dans l'hyperplan (pour voir si on est dans la phase de chargement statique poids propre) 338 : Xnorm_f25=sqrt(force(2)**2+force(3)**2+force(4)**2+force(5)**2) 339 : if (dtime.ne.0) then 340 : if ((eta3.ne.0) .and. (Xnorm_f25.gt.0.001)) then 341 : DO I=1,5 342 : force(I) = (Ftr(I)+dtime/eta3*force(I))/(1+dtime/eta3) 343 : upl(I) = (upl_preci(I)+dtime/eta3*upl(I))/(1+dtime/eta3) 344 : q1(I) = (q1_preci(I)+dtime/eta3*q1(I))/(1+dtime/eta3) 345 : ENDDO 346 : q2 = (q2_preci+dtime/eta3*q2)/(1+dtime/eta3) 347 : endif 348 : endif 349 : * 350 : IF (force(1).GE.1) THEN 351 : KERRE = 25 352 : RETURN 353 : ENDIF 354 : delta_preci=delta 355 : DO I=1,5 356 : force_preci(I)=force(I) 357 : ENDDO 358 : Xkapa_preci=Xkapa 359 : Xksi_preci=Xksi 360 : distmin_preci=distmin 361 : * 362 : enddo 363 : * 363 : call SCALPR(5,upl,upl,EPSE) 364 : EPSE=sqrt(EPSE) 365 : * 366 : VARF(1)=EPSE 367 : VARF(2)=upl1(1) 368 : VARF(3)=upl1(2) 369 : VARF(4)=upl1(3) 370 : VARF(5)=upl1(4) 371 : VARF(6)=upl1(5) 372 : VARF(7)=q1(1) 373 : VARF(8)=q1(2) 374 : VARF(9)=q1(3) 375 : VARF(10)=q1(4) 376 : VARF(11)=q1(5) 377 : VARF(12)=delta 378 : VARF(13)=Xkapa 379 : VARF(14)=Xksi 380 : VARF(15)=distmin 381 : *deltamax ne sert pas, on la garde juste pour la tracer (provisoire) 382 : VARF(16)=deltamax 383 : VARF(17)=upl2(1) 384 : VARF(18)=upl2(2) 385 : VARF(19)=upl2(3) 386 : VARF(20)=upl2(4) 387 : VARF(21)=upl2(5) 388 : VARF(22)=q2 389 : 390 : C Redimensionnement des variables 391 : if (a9.eq.1) then 392 : C semelle filante 393 : DO I=1,5 394 : force(I)=force(I)*Diam*qmax 395 : ENDDO 396 : force(3)=force(3)*Diam 397 : force(5)=force(5)*Diam 398 : 399 : else 400 : if (a9.eq.2) then 401 : C semelle rectangulaire 402 : DO I=1,5 403 : force(I)=force(I)*XLx*XLy*qmax 404 : ENDDO 405 : force(3)=force(3)*XLx 406 : force(5)=force(5)*XLy 407 : else 408 : C semelle circulaire 409 : DO I=1,5 410 : force(I)=force(I)*XPI*Diam**2/4*qmax 411 : ENDDO 412 : force(3)=force(3)*Diam 413 : force(5)=force(5)*Diam 414 : endif 415 : endif 416 : * 416 : call MKREDI(da,XMKel) 417 : do i=1,5 418 : do j=1,5 419 : XMK(i,j)=XMKel(i,j) 420 : enddo 421 : enddo 422 : * 423 : *Fz ISS --> -Fx joint 424 : SIGF(1)=-force(1) 425 : *Fx ISS --> -Fy joint 426 : SIGF(2)=-force(2) 427 : *Fy ISS --> Fz joint 428 : SIGF(3)=force(4) 429 : *Mz ISS --> -Mx joint (torsion linéaire) 430 : SIGF(4)=-XKelt*epst0(4) 431 : *Mx ISS --> -My joint 432 : SIGF(5)=-force(5) 433 : *My ISS --> Mz joint 434 : SIGF(6)=force(3) 435 : * 436 : return 437 : end 438 : 439 : 440 :
© Cast3M 2003 - All rights reserved.
Disclaimer