Télécharger la plaquette Cast3M


fr en

Télécharger bcoq4o.eso

Retour à la liste

Numérotation des lignes : 

bcoq4o

  1. C BCOQ4O    SOURCE    CHAT      06/03/29    21:15:44     5360
  2.       SUBROUTINE BCOQ4O
  3.      &          (IGAU,XEL,SHPTOT,SHP,BGENE,DJAC,EXCEN,NOPLAN,IARR,IFS)
  4.       IMPLICIT INTEGER(I-N)
  5.       IMPLICIT REAL*8(A-H,O-Z)
  6. ************************************************************************
  7. *
  8. *                             B C O Q 4
  9. *                             ---------
  10. *
  11. * FONCTION:
  12. * ---------
  13. *
  14. *    Calcul de la matrice "B" pour un COQ4 de comprtement non-isotrope,
  15. *                                                    excentré ou non.
  16. *
  17. * PARAMETRES:   (E)=ENTREE   (S)=SORTIE
  18. * -----------
  19. *
  20. *     IGAU    (E)  Numéro du point d'intégration.
  21. *     XEL     (E)  Coordonnées locales des noeuds de l'élément.
  22. *     SHPTOT  (E)  Fonctions de forme et dérivées dans l'espace de
  23. *                  référence.
  24. *                  SHPTOT(1,... = fonction,
  25. *                  SHPTOT(2,... = dérivée par rapport à Qsi,
  26. *                  SHPTOT(3,... = dérivée par rapport à Eta.
  27. *     EXCEN   (E)  Excentrement.
  28. *     NOPLAN  (E)  = 1 si élément non plan,
  29. *                  = 0 sinon.
  30. *     IFS     (E)  = 1 si demande de Matrice "B" complète
  31. *                  (pour SIGMA et EPSILON),
  32. *                  = 0 sinon.
  33. *     SHP     (S)  Fonctions de forme et dérivées dans l'espace
  34. *                  géométrique.
  35. *                  SHP(1,... = fonction,
  36. *                  SHP(2,... = dérivée par rapport à "X" local,
  37. *                  SHP(3,... = dérivée par rapport à "Y" local.
  38. *                  - au point d'intégration n.5 si IFS=1,
  39. *                  - au point d'intégration IGAU courant sinon.
  40. *     DJAC    (S)  Jacobien au point d'intégration IGAU.
  41. *     BGENE   (S)  Matrice "B" au point d'intégration IGAU.
  42. *     IARR    (S)  Différent de 0  si erreur.
  43. *
  44.       PARAMETER (LRE=24,NST=8,NBNO=4)
  45.       REAL*8 XEL(3,*),BGENE(NST,*),SHP(6,*),SHPTOT(6,NBNO,*)
  46. *
  47. * CONSTANTES:
  48. * -----------
  49. *
  50. *     NBNO   = Nombre de noeuds.
  51. *     LRE    = Nombre de colonnes de la Matrice "B".
  52. *     NST    = Nombre de composantes de contraintes.
  53. *
  54. * VARIABLES:
  55. * ----------
  56. *
  57. *     FULL   = .TRUE. si on veut une Matrice "B" complète.
  58. *     GAUSS  = .TRUE. si utilisation d'un point d'intégration n.1 à 4.
  59. *     SHP7   = Fonctions d'interpolation mixte pour le cisaillement X-Z
  60. *     SHP8   = Fonctions d'interpolation mixte pour le cisaillement Y-Z
  61. *
  62.       REAL*8 SHP7(4:5,NBNO),SHP8(4:5,NBNO)
  63.       LOGICAL FULL,GAUSS
  64. *
  65. * MODE DE FONCTIONNEMENT:
  66. * -----------------------
  67. *
  68. * POUR LA RIGIDITE:
  69. *
  70. *     Pour la membrane,
  71. *     - intégration pleine si élément plan,
  72. *     - intégration réduite si élément non plan.
  73. *
  74. *     Pour le cisaillement transverse, voir:
  75. *     DONEA / LAMAIN (1987) - ISPRA
  76. *     (élément de type "A.N.S.")
  77. *
  78. *     Pour le couplage membrane-flexion (excentrement), intégration
  79. *     réduite analogue à membrane, non justifiée au niveau de la
  80. *     Rigidité, mais indispensable pour le calcul pratique des forces
  81. *     nodales à partir des contraintes, par "BSIGMA".
  82. *
  83. * POUR LES CONTRAINTES:
  84. *
  85. *     Cependant, pour que le champ de contraintes produit soit facile
  86. *     à lire par l'utilisateur et à traiter par les opérateurs tels
  87. *     que:
  88. *     - CHANGER CHPOINT,
  89. *     - TRACER, ...
  90. *     on inscrit aux points de Gauss 1 à 4 les valeurs de contrainte
  91. *     au point d'intégration n.5 pour les parties concernées.
  92. *     De plus, on remplit complètement le champ de valeurs au point
  93. *     d'intégration n.5 .
  94. *
  95. * AUTEUR, DATE DE CREATION:
  96. * -------------------------
  97. *
  98. *     Michel BULIK       26 Novembre 1996 d'après BCOQ4 écrit par
  99. *     Pascal MANIGOT     17 Juin 1991
  100. *
  101. ************************************************************************
  102. *
  103.       IARR=0
  104.       FULL = IFS .EQ. 1
  105.       GAUSS = IGAU .LE. 4
  106. *
  107. *
  108. *     Matrice jacobienne:
  109. *     -------------------
  110. *
  111.       DO 110 I=1,NBNO
  112.          SHP(1,I)=SHPTOT(1,I,IGAU)
  113.          SHP(2,I)=SHPTOT(2,I,IGAU)
  114.          SHP(3,I)=SHPTOT(3,I,IGAU)
  115.   110    CONTINUE
  116. *     END DO
  117. *
  118.       DXDQSI = 0.D0
  119.       DXDETA = 0.D0
  120.       DYDQSI = 0.D0
  121.       DYDETA = 0.D0
  122.       DO 120 I=1,NBNO
  123.          DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I)
  124.          DXDETA=DXDETA+SHP(3,I)*XEL(1,I)
  125.          DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I)
  126.          DYDETA=DYDETA+SHP(3,I)*XEL(2,I)
  127.  120     CONTINUE
  128. *     END DO
  129.       DJAC=DXDQSI*DYDETA-DXDETA*DYDQSI
  130.       IF(DJAC.LE. 0.D0) THEN
  131.          IARR=1
  132.          RETURN
  133.       END IF
  134. *
  135. *     Par inversion analytique de la matrice J nous avons les dérivées
  136. *
  137.       DQSIDX= DYDETA/DJAC
  138.       DQSIDY=-DXDETA/DJAC
  139.       DETADX=-DYDQSI/DJAC
  140.       DETADY= DXDQSI/DJAC
  141. *
  142. *     puis les dérivées des fonctions de forme dans l'espace géométrique:
  143. *
  144.       DO 130 I=1,NBNO
  145.          TT=SHP(2,I)*DQSIDX+SHP(3,I)*DETADX
  146.          SHP(3,I)=SHP(2,I)*DQSIDY+SHP(3,I)*DETADY
  147.          SHP(2,I)=TT
  148.   130    CONTINUE
  149. *     END DO
  150. *
  151. *
  152. *     Matrice "B":
  153. *     ------------
  154. *
  155.       CALL ZERO(BGENE,NST,LRE)
  156. *
  157.       IF (GAUSS .OR. FULL) THEN
  158. *        Préparation pour cisaillement transverse:
  159. *
  160.          N1 = 1
  161.          N2 = 2
  162.          DX1 = SHPTOT(2,N1,IGAU)*XEL(1,N1)+SHPTOT(2,N2,IGAU)*XEL(1,N2)
  163.          DY1 = SHPTOT(2,N1,IGAU)*XEL(2,N1)+SHPTOT(2,N2,IGAU)*XEL(2,N2)
  164.          N1 = 3
  165.          N2 = 4
  166.          DX2 = SHPTOT(2,N1,IGAU)*XEL(1,N1)+SHPTOT(2,N2,IGAU)*XEL(1,N2)
  167.          DY2 = SHPTOT(2,N1,IGAU)*XEL(2,N1)+SHPTOT(2,N2,IGAU)*XEL(2,N2)
  168.          N1 = 1
  169.          N2 = 4
  170.          DX3 = SHPTOT(3,N1,IGAU)*XEL(1,N1)+SHPTOT(3,N2,IGAU)*XEL(1,N2)
  171.          DY3 = SHPTOT(3,N1,IGAU)*XEL(2,N1)+SHPTOT(3,N2,IGAU)*XEL(2,N2)
  172.          N1 = 2
  173.          N2 = 3
  174.          DX4 = SHPTOT(3,N1,IGAU)*XEL(1,N1)+SHPTOT(3,N2,IGAU)*XEL(1,N2)
  175.          DY4 = SHPTOT(3,N1,IGAU)*XEL(2,N1)+SHPTOT(3,N2,IGAU)*XEL(2,N2)
  176. *
  177.          D2JAC = 2.D0 * DJAC
  178. *
  179.          SHP7(4,1) = (- DYDETA*DY1 + DYDQSI*DY3) / D2JAC
  180.          SHP7(4,2) = (- DYDETA*DY1 + DYDQSI*DY4) / D2JAC
  181.          SHP7(4,3) = (- DYDETA*DY2 + DYDQSI*DY4) / D2JAC
  182.          SHP7(4,4) = (- DYDETA*DY2 + DYDQSI*DY3) / D2JAC
  183.          SHP7(5,1) = (  DYDETA*DX1 - DYDQSI*DX3) / D2JAC
  184.          SHP7(5,2) = (  DYDETA*DX1 - DYDQSI*DX4) / D2JAC
  185.          SHP7(5,3) = (  DYDETA*DX2 - DYDQSI*DX4) / D2JAC
  186.          SHP7(5,4) = (  DYDETA*DX2 - DYDQSI*DX3) / D2JAC
  187. *
  188.          SHP8(4,1) = (  DXDETA*DY1 - DXDQSI*DY3) / D2JAC
  189.          SHP8(4,2) = (  DXDETA*DY1 - DXDQSI*DY4) / D2JAC
  190.          SHP8(4,3) = (  DXDETA*DY2 - DXDQSI*DY4) / D2JAC
  191.          SHP8(4,4) = (  DXDETA*DY2 - DXDQSI*DY3) / D2JAC
  192.          SHP8(5,1) = (- DXDETA*DX1 + DXDQSI*DX3) / D2JAC
  193.          SHP8(5,2) = (- DXDETA*DX1 + DXDQSI*DX4) / D2JAC
  194.          SHP8(5,3) = (- DXDETA*DX2 + DXDQSI*DX4) / D2JAC
  195.          SHP8(5,4) = (- DXDETA*DX2 + DXDQSI*DX3) / D2JAC
  196. *
  197.       END IF
  198. *
  199.       DO 300 I=1,NBNO
  200.          K = 6*(I-1)
  201. *
  202.          IF (GAUSS .OR. FULL) THEN
  203.             IF (NOPLAN .EQ. 0) THEN
  204. *              Membrane:
  205.                BGENE(1,K+1)=SHP(2,I)
  206.                BGENE(2,K+2)=SHP(3,I)
  207. *              Couplage membrane-flexion:
  208.                BGENE(1,K+5) = SHP(2,I)  * EXCEN
  209.                BGENE(2,K+4) = -SHP(3,I) * EXCEN
  210. *              Cisaillement de membrane et de couplage membrane-flexion:
  211.                BGENE(3,K+1)=SHP(3,I)
  212.                BGENE(3,K+2)=SHP(2,I)
  213.                BGENE(3,K+4) = -SHP(2,I) * EXCEN
  214.                BGENE(3,K+5) = SHP(3,I)  * EXCEN
  215.             ENDIF
  216. *           Flexion:
  217.             BGENE(4,K+5)=SHP(2,I)
  218.             BGENE(5,K+4)=-SHP(3,I)
  219.             BGENE(6,K+4)=-SHP(2,I)
  220.             BGENE(6,K+5)=SHP(3,I)
  221. *           Cisaillement transverse:
  222.             BGENE(7,K+3)=SHP(2,I)
  223.             BGENE(7,K+4) = SHP7(4,I)
  224.             BGENE(7,K+5) = SHP7(5,I)
  225.             BGENE(8,K+3)=SHP(3,I)
  226.             BGENE(8,K+4) = SHP8(4,I)
  227.             BGENE(8,K+5) = SHP8(5,I)
  228.          ENDIF
  229. *
  230.          IF (.NOT. GAUSS) THEN
  231.             IF (NOPLAN.EQ.1 .OR. FULL) THEN
  232. *              Membrane:
  233.                BGENE(1,K+1)=SHP(2,I)
  234.                BGENE(2,K+2)=SHP(3,I)
  235. *              Couplage membrane-flexion:
  236.                BGENE(1,K+5) = SHP(2,I)  * EXCEN
  237.                BGENE(2,K+4) = -SHP(3,I) * EXCEN
  238. *              Cisaillement de membrane et de couplage membrane-flexion:
  239.                BGENE(3,K+1)=SHP(3,I)
  240.                BGENE(3,K+2)=SHP(2,I)
  241.                BGENE(3,K+4) = -SHP(2,I) * EXCEN
  242.                BGENE(3,K+5) = SHP(3,I)  * EXCEN
  243.             ENDIF
  244.          ENDIF
  245. *
  246.   300    CONTINUE
  247. *     END DO
  248. *
  249. *
  250.       IF (GAUSS .AND. FULL .AND. NOPLAN .EQ. 1) THEN
  251. *        On finit de compléter la Matrice "B", pour les contraintes ou
  252. *        les déformations, par des valeurs calculées au centre:
  253. *
  254.          DO 500 I=1,NBNO
  255.             SHP(1,I)=SHPTOT(1,I, 5  )
  256.             SHP(2,I)=SHPTOT(2,I, 5  )
  257.             SHP(3,I)=SHPTOT(3,I, 5  )
  258.   500       CONTINUE
  259. *        END DO
  260. *        FONCTIONS DE FORME DANS LA GÉOMÉTRIE RÉELLE:
  261.          CALL JACOBI(XEL,SHP,2,NBNO,DJAC5)
  262.          IF(DJAC5 .LE. 0.D0) THEN
  263.             IARR=1
  264.             RETURN
  265.          END IF
  266. *
  267.          DO 520 I=1,NBNO
  268.             K = 6*(I-1)
  269. *           Membrane:
  270.             BGENE(1,K+1)=SHP(2,I)
  271.             BGENE(2,K+2)=SHP(3,I)
  272. *           Couplage membrane-flexion:
  273.             BGENE(1,K+5) = SHP(2,I)  * EXCEN
  274.             BGENE(2,K+4) = -SHP(3,I) * EXCEN
  275. *           Cisaillement de membrane et de couplage membrane-flexion:
  276.             BGENE(3,K+1)=SHP(3,I)
  277.             BGENE(3,K+2)=SHP(2,I)
  278.             BGENE(3,K+4) = -SHP(2,I) * EXCEN
  279.             BGENE(3,K+5) = SHP(3,I)  * EXCEN
  280.   520       CONTINUE
  281. *        END DO
  282. *
  283.       ENDIF
  284. *
  285. ccc      write(*,*) 'BCOQ4O: Matrice B ='
  286. ccc      write(*,'(8(1X,1P,G10.3))') ((bgene(iii,jjj),jjj=1,lre),iii=1,nst)
  287.       END
  288.  
  289.  
  290.  
  291.  

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