Télécharger bcoq4o.eso

Retour à la liste

Numérotation des lignes :

  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