Télécharger la plaquette Cast3M


fr en

Télécharger issgra.eso

Retour à la liste

Numérotation des lignes : 

issgra

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales