C ISSGRA SOURCE PV 17/12/08 21:17:29 9660 C SUBROUTINE ISSGRA(WRK52,WRK53,WRK54,WRK27,IB,IGAU,NBPGAU) C C______________________________________________________________________ C C MODELE ISSGRA C______________________________________________________________________ C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) c inutile? -INC PPARAM -INC CCOPTIO -INC CCREEL -INC DECHE * C GLOBAL MODEL DIMENSION XMKel(5,5) C DIMENSION XMCel(5,5) DIMENSION uo(5), duo(5), Fo(5) DIMENSION XMK(5,5), Ftr(5), force(5) * DIMENSION upl(5), uel(5) DIMENSION uoi(5), duoi(5), dueli(5), force_preci(5) C DIMENSION ueli(5) C Modèle de Plasticité DIMENSION ddlambda(2), yield(2) DIMENSION Ho(2,2), H(2,2), Hinv(2,2) DIMENSION duP(5) * C Modèle de plasticité 1 DIMENSION q1(5), upl1(5) C DIMENSION q1_prec(5) DIMENSION Gradgs(5),Gradfs(5),Gradfq(5) DIMENSION hh1(5), XKGradgs(5) * C Modèle de plasticité vertical DIMENSION upl2(5), Gradf2s(5), XKGradf2s(5) * C Modèle de Décollement DIMENSION dforce(5), forcepred(5), force_prec(5), duel(5) * C Amortissment DIMENSION upl_preci(5),q1_preci(5) * DIMENSION da(20) C DIMENSION TIME(2) * C----------------------------------------------------------------------- C----------------------------------------------------------------------- C Affectation des données du modèle da(1) =XMAT(13) da(2) =XMAT(14) da(3) =XMAT(15) da(4) =XMAT(7) da(5) =XMAT(8) da(6) =XMAT(12) da(7) =XMAT(11) da(20)=XMAT(10) da(8) =XMAT(16) da(9) =XMAT(17) da(10)=XMAT(18) da(11)=XMAT(19) da(12)=XMAT(20) da(13)=XMAT(21) da(14)=XMAT(22) da(15)=XMAT(23) da(16)=XMAT(24) da(17)=XMAT(25) da(18)=XMAT(26) da(19)=XMAT(27) * Diam=da(1) XLx=da(2) XLy=da(3) XKelz=da(4) XKelh=da(5) XKelry=da(6) XKelrx=da(7) C a=da(8) C b=da(9) C c=da(10) C d=da(11) C e=da(12) C f=da(13) qmax=da(14) a9=da(15) C a6=da(16) eta3=da(17) dtime=da(18) a8=da(19) XKelt=da(20) * * upl1(1)=VAR0(2) upl1(2)=VAR0(3) upl1(3)=VAR0(4) upl1(4)=VAR0(5) upl1(5)=VAR0(6) q1(1)=VAR0(7) q1(2)=VAR0(8) q1(3)=VAR0(9) q1(4)=VAR0(10) q1(5)=VAR0(11) force_prec(1)=-SIG0(1) force_prec(2)=-SIG0(2) force_prec(3)=SIG0(6) force_prec(4)=SIG0(3) force_prec(5)=-SIG0(5) delta_prec=VAR0(12) Xkapa_prec=VAR0(13) Xksi_prec=VAR0(14) distmin_prec=VAR0(15) upl2(1)=VAR0(17) upl2(2)=VAR0(18) upl2(3)=VAR0(19) upl2(4)=VAR0(20) upl2(5)=VAR0(21) q2=VAR0(22) DO I=1,5 upl(I)=upl1(I)+upl2(I) ENDDO * *deplx joint --> -deplz ISS uo(1)=-epst0(1) *deply joint --> -deplx ISS uo(2)=-epst0(2) *rotz joint --> -roty ISS uo(3)=-epst0(6) *deplz joint --> deply ISS uo(4)=epst0(3) *roty joint --> -rotx ISS uo(5)=-epst0(5) * duo(1)=-DEPST(1) duo(2)=-DEPST(2) duo(3)=-DEPST(6) duo(4)=DEPST(3) duo(5)=-DEPST(5) C Correction données d'entrées Castem DO I=1,5 uo(I)=uo(I)+duo(I) * epstf(I) = epst0(I)+DEPST(I) ENDDO * C on adimensionne les déplacements if (a9.eq.2) then DO I=1,5 uo(I) = uo(I)*sqrt(XLx**2+XLy**2)/(XLx*XLy) ENDDO uo(3) = uo(3)*XLx uo(5) = uo(5)*XLy DO I=1,5 duo(I) = duo(I)*sqrt(XLx**2+XLy**2)/(XLx*XLy) ENDDO duo(3)= duo(3)*XLx duo(5)= duo(5)*XLy else DO I=1,5 uo(I) = uo(I)/Diam ENDDO uo(3) = uo(3)*Diam uo(5) = uo(5)*Diam DO I=1,5 duo(I) = duo(I)/Diam ENDDO duo(3)= duo(3)*Diam duo(5)= duo(5)*Diam endif * C Adimensionnement des efforts if (a9.eq.1) then C semelle filante DO I=1,5 force_prec(I)=force_prec(I)/(Diam*qmax) ENDDO force_prec(3)=force_prec(3)/Diam force_prec(5)=force_prec(5)/Diam else if (a9.eq.2) then C semelle rectangulaire DO I=1,5 force_prec(I)=force_prec(I)/(XLx*XLy*qmax) ENDDO force_prec(3)=force_prec(3)/XLx force_prec(5)=force_prec(5)/XLy else C semelle circulaire DO I=1,5 force_prec(I)=force_prec(I)/(XPI*Diam**2/4*qmax) ENDDO force_prec(3)=force_prec(3)/Diam force_prec(5)=force_prec(5)/Diam endif endif C *********************************************************************** C Macro-élément d'ISS C adimensionnement de la matrice de rigidité élastique suivant la forme de la fondation (paramètre a9) call MKADIM(da,XMKel) do i=1,5 do j=1,5 XMK(i,j)=XMKel(i,j) enddo enddo C Calcul de l'effort de prédiction élastique en considérant le décollement élastique non-linéaire. duP(1)=0 duP(2)=0 duP(3)=0 duP(4)=0 duP(5)=0 C Réduction du pas de calcul dtimei=dtime/10. DO I=1,5 duoi(I)=1./10.*duo(I) force_preci(I)=force_prec(I) ENDDO delta_preci=delta_prec Xkapa_preci=Xkapa_prec Xksi_preci=Xksi_prec distmin_preci=distmin_prec do k=1,10 q2_preci=q2 * Xkapa_k=Xkapa_preci Xksi_k=Xksi_preci DO I=1,5 q1_preci(I)=q1(I) upl_preci(I)=upl(I) uoi(I)=uo(I)-(10.-k)*duoi(I) uel(I)=uoi(I)-upl(I) duel(I)=duoi(I)-duP(I) ENDDO *si a8=1 on désactive le décollement if (a8.eq.1) then call MATVE1 (XMK, uel,5,5,Ftr,2) delta=delta_prec deltamax=1 else call MATVE1 (XMK, duel,5,5,dforce,2) DO I=1,5 force(I)=force_preci(I)+dforce(I) ENDDO call ELNLIN (da, uel, force, XMKel, & delta_preci, Ftr, delta, deltamax) endif C Calcul des fonctions de charge au point de charge élastique call LOADSF (da, Ftr,q1,yield1) call LOASFV (Ftr,q2,yield2) DO I=1,5 force(I)=Ftr(I) ENDDO C Boucle de plasticité (Return Mapping) epsilon=1E-9 if ((yield1.le.epsilon) .and. (yield2.le.epsilon)) then DO I=1,5 force(I)=Ftr(I) ENDDO Xkapa=Xkapa_k Xksi=Xksi_k distmin=distmin_preci else m=0 do while ((yield1.gt.epsilon) .or. (yield2.gt.epsilon)) m=m+1 IF (m.gt.1000) THEN KERRE = 22 RETURN ENDIF IF ((force(1).LT.0).or.(force(1).GE.1)) THEN KERRE = 23 RETURN ENDIF call GRADIF(da,force,q1,Gradfs,Gradfq) call GRADIG(da,force, q1, force_preci, Xkapa_k, & Xksi_k, Xkapa_preci, Xksi_preci, Gradgs, Xkapa, Xksi) call FVEH(da,force,q1, XMKel, & Gradgs,distmin_preci,hh1,distmin) call GRADFV(force,q2,Gradf2s,Gradf2q) call FVEHV(force,q2, XMKel,Gradf2s,hh2) call MATVE1(XMK,Gradgs,5,5,XKGradgs,2) call MATVE1(XMK,Gradf2s,5,5,XKGradf2s,2) call SCALPR(5,Gradfs,XKGradgs,Ho11) call SCALPR(5,Gradfs,XKGradf2s,Ho12) call SCALPR(5,Gradf2s,XKGradgs,Ho21) call SCALPR(5,Gradf2s,XKGradf2s,Ho22) call SCALPR(5,Gradfq,hh1,H1) H2=Gradf2q*hh2 H(1,1)=Ho11+H1 H(1,2)=Ho12 H(2,1)=Ho21 H(2,2)=Ho22+H2 yield(1)=yield1 yield(2)=yield2 *initialisation de ces variables necessaire pour les tests sur les mécanismes ddlambda(1)=1 ddlambda(2)=1 if ((yield(1).gt.epsilon). & and.(yield(2).gt.epsilon)) then call HINV22 (H,Hinv) call MATVE1(Hinv,yield,2,2,ddlambda,2) ddlambda1=ddlambda(1) ddlambda2=ddlambda(2) endif C Tests sur les mécanismes if ((yield(1).le.epsilon).or.(ddlambda(1).le.0)) then ddlambda1=0 ddlambda2=yield2/H(2,2) if (ddlambda2.le.0) then ddlambda2=0 endif endif if ((yield(2).le.epsilon).or.(ddlambda(2).le.0)) then ddlambda2=0 ddlambda1=yield1/H(1,1) if (ddlambda1.le.0) then ddlambda1=0 endif endif ddlambda(1)=ddlambda1 ddlambda(2)=ddlambda2 q2=q2-ddlambda(2)*hh2 DO I=1,5 q1(I)=q1(I)-ddlambda(1)*hh1(I) upl1(I)=upl1(I)+ddlambda(1)*Gradgs(I) upl2(I)=upl2(I)+ddlambda(2)*Gradf2s(I) upl(I) =upl1(I)+upl2(I) * uel(I)=uoi(I)-upl(I) duP(I)=upl(I)-upl_preci(I) duel(I)=duoi(I)-duP(I) ENDDO *si a8=1 on désactive le décollement if (a8.eq.1) then call MATVE1 (XMK, uel,5,5, force,2) delta=delta_prec else call MATVE1 (XMK, duel, 5,5, dforce,2) DO I=1,5 forcepred(I)=force_preci(I)+dforce(I) ENDDO *forcepred et delta_prec sont là pour une initialisation à la boucle de Newton call ELNLIN (da, uel, forcepred, XMKel, & delta_preci, force, delta, deltamax) endif Xkapa_k=Xkapa Xksi_k=Xksi call LOADSF (da,force,q1,yield1) call LOASFV (force,q2,yield2) enddo endif * C Correction visqueuse * calcul de la norme de la force dans l'hyperplan (pour voir si on est dans la phase de chargement statique poids propre) Xnorm_f25=sqrt(force(2)**2+force(3)**2+force(4)**2+force(5)**2) if (dtime.ne.0) then if ((eta3.ne.0) .and. (Xnorm_f25.gt.0.001)) then DO I=1,5 force(I) = (Ftr(I)+dtime/eta3*force(I))/(1+dtime/eta3) upl(I) = (upl_preci(I)+dtime/eta3*upl(I))/(1+dtime/eta3) q1(I) = (q1_preci(I)+dtime/eta3*q1(I))/(1+dtime/eta3) ENDDO q2 = (q2_preci+dtime/eta3*q2)/(1+dtime/eta3) endif endif * IF (force(1).GE.1) THEN KERRE = 25 RETURN ENDIF delta_preci=delta DO I=1,5 force_preci(I)=force(I) ENDDO Xkapa_preci=Xkapa Xksi_preci=Xksi distmin_preci=distmin * enddo * call SCALPR(5,upl,upl,EPSE) EPSE=sqrt(EPSE) * VARF(1)=EPSE VARF(2)=upl1(1) VARF(3)=upl1(2) VARF(4)=upl1(3) VARF(5)=upl1(4) VARF(6)=upl1(5) VARF(7)=q1(1) VARF(8)=q1(2) VARF(9)=q1(3) VARF(10)=q1(4) VARF(11)=q1(5) VARF(12)=delta VARF(13)=Xkapa VARF(14)=Xksi VARF(15)=distmin *deltamax ne sert pas, on la garde juste pour la tracer (provisoire) VARF(16)=deltamax VARF(17)=upl2(1) VARF(18)=upl2(2) VARF(19)=upl2(3) VARF(20)=upl2(4) VARF(21)=upl2(5) VARF(22)=q2 C Redimensionnement des variables if (a9.eq.1) then C semelle filante DO I=1,5 force(I)=force(I)*Diam*qmax ENDDO force(3)=force(3)*Diam force(5)=force(5)*Diam else if (a9.eq.2) then C semelle rectangulaire DO I=1,5 force(I)=force(I)*XLx*XLy*qmax ENDDO force(3)=force(3)*XLx force(5)=force(5)*XLy else C semelle circulaire DO I=1,5 force(I)=force(I)*XPI*Diam**2/4*qmax ENDDO force(3)=force(3)*Diam force(5)=force(5)*Diam endif endif * call MKREDI(da,XMKel) do i=1,5 do j=1,5 XMK(i,j)=XMKel(i,j) enddo enddo * *Fz ISS --> -Fx joint SIGF(1)=-force(1) *Fx ISS --> -Fy joint SIGF(2)=-force(2) *Fy ISS --> Fz joint SIGF(3)=force(4) *Mz ISS --> -Mx joint (torsion linéaire) SIGF(4)=-XKelt*epst0(4) *Mx ISS --> -My joint SIGF(5)=-force(5) *My ISS --> Mz joint SIGF(6)=force(3) * return end