Télécharger hbmz.eso

Retour à la liste

Numérotation des lignes :

hbmz
  1. C HBMZ SOURCE CB215821 23/01/25 21:15:23 11573
  2.  
  3. SUBROUTINE HBMZ(NT,NHBM,NDDL,KTQ,KTKAM,LOGAMO)
  4.  
  5. *=======================================================================
  6. * Calcule la matrice de raideur dynamique Z
  7. *
  8. * Z = [ Z0 0 0 ... ]
  9. * [ Z1 0 ... ]
  10. * [ ... ]
  11. * [ Zh ]
  12. *
  13. * avec Zj = | [ K-(jw)²M jwC ] si LOGAMO
  14. * | [ -jwC K-(jw)²M ]
  15. * |
  16. * | [ K-(jw)²M ] si .not.LOGAMO
  17. * | [ K-(jw)²M ]
  18. *
  19. * TODO : prevoir le cas des matrices non diagonales
  20. *=======================================================================
  21. *
  22. *----- Declarations ----------------------------------------------------
  23. *
  24. IMPLICIT INTEGER(I-N)
  25. IMPLICIT REAL*8(A-H,O-Z)
  26.  
  27. INTEGER NT,NHBM,NDDL,I,J,K,L,M,P,PK,KM,KK,KASM
  28. REAL*8 AA,BB
  29. LOGICAL LOGAMO
  30. *
  31. *-INC TMDYNC.INC
  32. ************************** debut TMDYNC.INC ****************************
  33.  
  34. * TMDYNC : FUTUR INCLUDE POUR LES SEGMENTS DE L'OPERATEUR DYNC
  35. * TODO : a extraire dans un include des que stabilise
  36. *
  37. * Segment des variables generalisees:
  38. * -----------------------------------
  39. SEGMENT MTQ
  40. REAL*8 Q1(NT1)
  41. REAL*8 OMEG,XPARA
  42. REAL*8 JAC(NT1,NT1),ZZ(NT1,NT1),RX(NT1,NT1)
  43. REAL*8 dX(NT1), dw, dv
  44. ENDSEGMENT
  45. * Q1 : vecteur des inconnues frequentielles de dimension (2h+1)*n
  46. * Q1 = {q_0 q_c1 q_s1 ... q_sh}
  47. * avec q_i vecteur de dimension n ou n=nombre de modes
  48. * OMEG : frequence fondamentale de l'approximation
  49. * XPARA: parametre de continuation (par defaut la frequence)
  50. * \in [PARINI,PARFIN]
  51. * RX : matrice jacobienne = ZZ + dFnl/dX
  52. * JAC : jacobienne des efforts non-lineaires = dFnl/dX
  53. * ZZ : matrice dynamique associee aux matrices modales K, M et C
  54. * lineaires et constantes
  55. * {dX,dw,(dv)} : vecteur tangent utilise pour la prediction
  56. *
  57. *
  58. * Segment contenant les matrices XK, XASM et XM:
  59. * ---------------------------------------------
  60. SEGMENT MTKAM
  61. REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  62. REAL*8 GAM(NPC1,nl1),IGAM(nl1,NPC1),DL(nl1)
  63. * REAL*8 GAMFIN(NPC2,nl1)
  64. ENDSEGMENT
  65. * XK,XASM et XM : matrices de raideur, amortissement et masse
  66. * GAM et IGAM : matrices pour la FFT et son inverse
  67. * GAMFIN :
  68. *
  69. * Segment des deformees modales:
  70. * ------------------------------
  71. * (idem DYNE)
  72. SEGMENT MTPHI
  73. INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  74. INTEGER IAROTA(NSB)
  75. REAL*8 XPHILB(NSB,NPLSB,NA2,IDIMB)
  76. ENDSEGMENT
  77. *
  78. * Segment descriptif des liaisons en base A:
  79. * ------------------------------------------
  80. * (idem DYNE)
  81. SEGMENT MTLIAA
  82. INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
  83. REAL*8 XPALA(NLIAA,NXPALA)
  84. ENDSEGMENT
  85. *
  86. * Segment descriptif des liaisons en base B:
  87. * ------------------------------------------
  88. * (idem DYNE)
  89. SEGMENT MTLIAB
  90. INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  91. REAL*8 XPALB(NLIAB,NXPALB)
  92. REAL*8 XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  93. ENDSEGMENT
  94. *
  95. * Segment representant les chargements exterieurs:
  96. * -----------------------------------------------
  97. SEGMENT MTFEX
  98. REAL*8 FEXA(NT1)
  99. REAL*8 FEXPSM(NPLB,NPC1,2,IDIMB)
  100. INTEGER BAL
  101. ENDSEGMENT
  102. * FEXA : Vecteur des efforts ext. sous la forme de coefficients de
  103. * Fourier et exprimes en base A
  104. * FEXPSM: chargement/deplacement statique lie aux modes negliges
  105. * (neglige aussi les Fnl). Dans DYNC toujours =0, cree pour
  106. * compatibilite avec calcul des Fnl.
  107. * BAL : indique s'il s'agit d'un chargement de type balourd
  108. * (cad proportionnel a OMEG**2)
  109. *
  110. * Segment "local" pour DEVLFA:
  111. * ----------------------------
  112. SEGMENT LOCLFA
  113. REAL*8 FTEST(NA1,4)
  114. ENDSEGMENT
  115. *
  116. * Segment "local" pour DEVLB1:
  117. * ----------------------------
  118. SEGMENT LOCLB1
  119. REAL*8 FTEST2(NPLB,6)
  120. ENDSEGMENT
  121. *
  122. * Segment contenant les variables au cours d un pas de temps:
  123. * ----------------------------------------------------------
  124. SEGMENT MTPAS
  125. REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
  126. REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4)
  127. REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
  128. REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
  129. REAL*8 KTOTXA(NA1,NA1),KTOTVA(NA1,NA1)
  130. REAL*8 KTOTXB(NPLB,IDIMB,IDIMB), KTOTVB(NPLB,IDIMB,IDIMB)
  131. ENDSEGMENT
  132. * FTOTA/B/BA : forces sur base A, B et B projetees sur A
  133. * XPTB : deplacement du point d'une liaison en base B
  134. * XVALA/B : grandeurs de la liaison en base A/B a stocker
  135. * FEXB : forces exterieures en base B (a priori uniquement
  136. * pour les moments appliques aux rotations rigides ?)
  137. * XCHPFB : forces de contact en base B (lorsqu'on considere un
  138. * maillage de contact dans certaines liaisons)
  139. * KTOTXA/XB/VA/VB : Jacobienne par rapport au deplacement/vitesse en
  140. * base A/B (= contributions a dFnl/dX)
  141. *
  142. *
  143. * Segment des points de reference des modes (base A):
  144. * --------------------------------------------------
  145. SEGMENT MPREF
  146. INTEGER IPOREF(NPREF)
  147. ENDSEGMENT
  148. *
  149. * Segment des points en base B:
  150. * -----------------------------
  151. SEGMENT NCPR(NBPTS)
  152. * NCRP(#global) = #local dans XPTB (1er indice)
  153. *
  154. * Segment des parametres numeriques pour la continuation:
  155. * ------------------------------------------------------
  156. SEGMENT PARNUM
  157. CHARACTER*4 TYPS
  158. REAL*8 DS,DSMAX,DSMIN,ANGMIN,ANGMAX,ITERMOY,ISENS,TOLMIN
  159. REAL*8 PARINI,PARFIN
  160. INTEGER ITERMAX,NBPAS
  161. LOGICAL JANAL
  162. ENDSEGMENT
  163. *
  164. * Segment des resultats:
  165. * ---------------------
  166. SEGMENT PSORT
  167. REAL*8 QSAVE(NT1,NPAS),WSAVE(NPAS),LSAVE(2,2*NA1,NPAS)
  168. REAL*8 VSAVE(NPAS)
  169. LOGICAL ZSAVE(NPAS)
  170. CHARACTER*2 TYPBIF(NBIFU)
  171. REAL*8 QBIFU(NT1,NBIFU),WBIFU(NBIFU),WBIF2(NBIFU)
  172. REAL*8 QPSIR(NT1,NBIFU),QPSII(NT1,NBIFU)
  173. INTEGER CBIF
  174. ENDSEGMENT
  175. * QSAVE(i,j) = Q harmonique i au pas j
  176. * VSAVE(j) = parametre de continuation (si non w) au j-eme pas
  177. * ZSAVE(j) = stabilite au j-eme pas
  178. * LSAVE(1,j) : partie reelle de l'exposant de Floquet
  179. * LSAVE(2,j) : partie imaginaire de l'exposant de Floquet
  180. * TYPBIF = {LimitPoint, BranchPoint, NeimarkSacker, PeriodDoubling}
  181. * QBIFU,WBIFU : vecteur Q et w au point de bifurcation
  182. * WBIF2 : partie imaginaire de l'exposant de Floquet
  183. * QPSIR,QPSII : vecteur propre au point de bifurcation
  184.  
  185. * Segment des tableaux de travail:
  186. * -------------------------------
  187. SEGMENT MTEMP
  188. REAL*8 RW(NT1),A,T0(NT1+1),TP(NT1+1),AMPX,AUX
  189. REAL*8 T02(NT1+2), TP2(NT1+2)
  190. INTEGER IPIV(NT1),IPIV2(NT1+1),IPIV3(NT1+2)
  191. REAL*8 res
  192. REAL*8 RHS(NT1+1),Ja(NT1+1,NT1+1)
  193. REAL*8 QOLD(NT1),OMEGOLD
  194. REAL*8 MATJA(NT1+1,NT1+1),Rw2(NT1)
  195. REAL*8 Jaa(NT1+2,NT1+2),RHS2(NT1+2),Ra(NT1),VV,VVOLD
  196. ENDSEGMENT
  197. * Jacobiennes augmentees
  198. * Ja : [ RX Rw ; dX dw]
  199. * Jaa: [ RX Rw Ra; gx 0 0; dX dw da]
  200.  
  201. * SEGMENT NNNN
  202. * REAL*8 IGAM2(nl1,NPC2),DL2(nl1)
  203. * ENDSEGMENT
  204.  
  205. *************************** fin TMDYNC.INC *****************************
  206.  
  207. *
  208. * Segment des variables
  209. MTQ=KTQ
  210. * Segment des proprietes mecaniques
  211. MTKAM=KTKAM
  212.  
  213. *----- Initialisation --------------------------------------------------
  214. *
  215. DO J = 1,NT
  216. DO I = 1,NT
  217. ZZ(I,J)=0.D0
  218. ENDDO
  219. ENDDO
  220.  
  221.  
  222. *----- Remplissage (cas K, C et M diagonales uniquement) ---------------
  223.  
  224. * test ici mais a replacer an amont + tard
  225. if(XK(/2).ne.1) call erreur(5)
  226.  
  227. c harmonique 0
  228. DO M=1,NDDL
  229. ZZ(M,M) = XK(M,1)
  230. ENDDO
  231.  
  232. c harmoniques > 0
  233. c -cas non conservatif (i.e. avec amortissement)
  234. IF (LOGAMO) THEN
  235. DO J=2,2*NHBM,2
  236. DO I=1,NDDL
  237. AA = XK(I,1) - ((OMEG*J/2)**2)*XM(I,1)
  238. BB = (OMEG*J/2)*XASM(I,1)
  239. ZZ(NDDL*(1+(J-2))+I,NDDL*(1+(J-2))+I) = AA
  240. ZZ(NDDL*(1+(J-2))+I,NDDL*(1+(J-1))+I) = BB
  241. ZZ(NDDL*(1+(J-1))+I,NDDL*(1+(J-2))+I) = -BB
  242. ZZ(NDDL*(1+(J-1))+I,NDDL*(1+(J-1))+I) = AA
  243. ENDDO
  244. ENDDO
  245. c -cas conservatif (i.e. sans amortissement)
  246. ELSE
  247. DO J=2,2*NHBM,2
  248. DO I=1,NDDL
  249. AA = XK(I,1) - ((OMEG*J/2)**2)*XM(I,1)
  250. ZZ(NDDL*(1+(J-2))+I,NDDL*(1+(J-2))+I) = AA
  251. ZZ(NDDL*(1+(J-2))+I,NDDL*(1+(J-1))+I) = 0.D0
  252. ZZ(NDDL*(1+(J-1))+I,NDDL*(1+(J-2))+I) = 0.D0
  253. ZZ(NDDL*(1+(J-1))+I,NDDL*(1+(J-1))+I) = AA
  254. ENDDO
  255. ENDDO
  256. ENDIF
  257.  
  258. END
  259.  
  260.  
  261.  
  262.  
  263.  

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