Télécharger hbmfn2.eso

Retour à la liste

Numérotation des lignes :

hbmfn2
  1. C HBMFN2 SOURCE CB215821 23/01/25 21:15:20 11573
  2.  
  3. SUBROUTINE HBMFN2(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,KTPHI,KTLIAA,
  4. & KTLIAB,KTPAS,KTFEX,KOCLFA,KOCLB1,FNL)
  5.  
  6. *=======================================================================
  7. * Calcul de FNL(X) et dFNLdX par AFT (Alternating Frequency Time)
  8. *=======================================================================
  9.  
  10. *----- Declarations ----------------------------------------------------
  11.  
  12. IMPLICIT INTEGER(I-N)
  13. IMPLICIT REAL*8(A-H,O-Z)
  14. INTEGER NP,NDDL,NT,NFFT
  15. segment mwork
  16. REAL*8 XT(NDDL,NFFT),FT(NDDL,NFFT),VT(NDDL,NFFT)
  17. REAL*8 XAUX(NDDL,4),VAUX(NDDL,4),VCX(NFFT),VCV(NFFT)
  18. REAL*8 MATJX(NDDL,NDDL,NFFT), MATJV(NDDL,NDDL,NFFT)
  19. REAL*8 auxDL(2*NHBM+1,NFFT), JFR(2*NHBM+1,2*NHBM+1)
  20. REAL*8 auxVDL(2*NHBM+1,NFFT), JVFR(2*NHBM+1,2*NHBM+1)
  21. REAL*8 JACV(NT,NT),FNL2(NT)
  22. endsegment
  23. c segment to call fftpack5.1
  24. segment wfft51
  25. real*8 work(lenwrk)
  26. real*8 wwsave(lensav)
  27. real*8 XX51(NCOU)
  28. real*8 XV51(NCOU)
  29. endsegment
  30. REAL*8 PDT,FNL(NT),V1(NT)
  31. LOGICAL RIGIDE
  32. *
  33. *-INC TMDYNC.INC
  34. ************************** debut TMDYNC.INC ****************************
  35.  
  36. * TMDYNC : FUTUR INCLUDE POUR LES SEGMENTS DE L'OPERATEUR DYNC
  37. * TODO : a extraire dans un include des que stabilise
  38. *
  39. * Segment des variables generalisees:
  40. * -----------------------------------
  41. SEGMENT MTQ
  42. REAL*8 Q1(NT1)
  43. REAL*8 OMEG,XPARA
  44. REAL*8 JAC(NT1,NT1),ZZ(NT1,NT1),RX(NT1,NT1)
  45. REAL*8 dX(NT1), dw, dv
  46. ENDSEGMENT
  47. * Q1 : vecteur des inconnues frequentielles de dimension (2h+1)*n
  48. * Q1 = {q_0 q_c1 q_s1 ... q_sh}
  49. * avec q_i vecteur de dimension n ou n=nombre de modes
  50. * OMEG : frequence fondamentale de l'approximation
  51. * XPARA: parametre de continuation (par defaut la frequence)
  52. * \in [PARINI,PARFIN]
  53. * RX : matrice jacobienne = ZZ + dFnl/dX
  54. * JAC : jacobienne des efforts non-lineaires = dFnl/dX
  55. * ZZ : matrice dynamique associee aux matrices modales K, M et C
  56. * lineaires et constantes
  57. * {dX,dw,(dv)} : vecteur tangent utilise pour la prediction
  58. *
  59. *
  60. * Segment contenant les matrices XK, XASM et XM:
  61. * ---------------------------------------------
  62. SEGMENT MTKAM
  63. REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  64. REAL*8 GAM(NPC1,nl1),IGAM(nl1,NPC1),DL(nl1)
  65. * REAL*8 GAMFIN(NPC2,nl1)
  66. ENDSEGMENT
  67. * XK,XASM et XM : matrices de raideur, amortissement et masse
  68. * GAM et IGAM : matrices pour la FFT et son inverse
  69. * GAMFIN :
  70. *
  71. * Segment des deformees modales:
  72. * ------------------------------
  73. * (idem DYNE)
  74. SEGMENT MTPHI
  75. INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  76. INTEGER IAROTA(NSB)
  77. REAL*8 XPHILB(NSB,NPLSB,NA2,IDIMB)
  78. ENDSEGMENT
  79. *
  80. * Segment descriptif des liaisons en base A:
  81. * ------------------------------------------
  82. * (idem DYNE)
  83. SEGMENT MTLIAA
  84. INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
  85. REAL*8 XPALA(NLIAA,NXPALA)
  86. ENDSEGMENT
  87. *
  88. * Segment descriptif des liaisons en base B:
  89. * ------------------------------------------
  90. * (idem DYNE)
  91. SEGMENT MTLIAB
  92. INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  93. REAL*8 XPALB(NLIAB,NXPALB)
  94. REAL*8 XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  95. ENDSEGMENT
  96. *
  97. * Segment representant les chargements exterieurs:
  98. * -----------------------------------------------
  99. SEGMENT MTFEX
  100. REAL*8 FEXA(NT1)
  101. REAL*8 FEXPSM(NPLB,NPC1,2,IDIMB)
  102. INTEGER BAL
  103. ENDSEGMENT
  104. * FEXA : Vecteur des efforts ext. sous la forme de coefficients de
  105. * Fourier et exprimes en base A
  106. * FEXPSM: chargement/deplacement statique lie aux modes negliges
  107. * (neglige aussi les Fnl). Dans DYNC toujours =0, cree pour
  108. * compatibilite avec calcul des Fnl.
  109. * BAL : indique s'il s'agit d'un chargement de type balourd
  110. * (cad proportionnel a OMEG**2)
  111. *
  112. * Segment "local" pour DEVLFA:
  113. * ----------------------------
  114. SEGMENT LOCLFA
  115. REAL*8 FTEST(NA1,4)
  116. ENDSEGMENT
  117. *
  118. * Segment "local" pour DEVLB1:
  119. * ----------------------------
  120. SEGMENT LOCLB1
  121. REAL*8 FTEST2(NPLB,6)
  122. ENDSEGMENT
  123. *
  124. * Segment contenant les variables au cours d un pas de temps:
  125. * ----------------------------------------------------------
  126. SEGMENT MTPAS
  127. REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
  128. REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4)
  129. REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
  130. REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
  131. REAL*8 KTOTXA(NA1,NA1),KTOTVA(NA1,NA1)
  132. REAL*8 KTOTXB(NPLB,IDIMB,IDIMB), KTOTVB(NPLB,IDIMB,IDIMB)
  133. ENDSEGMENT
  134. * FTOTA/B/BA : forces sur base A, B et B projetees sur A
  135. * XPTB : deplacement du point d'une liaison en base B
  136. * XVALA/B : grandeurs de la liaison en base A/B a stocker
  137. * FEXB : forces exterieures en base B (a priori uniquement
  138. * pour les moments appliques aux rotations rigides ?)
  139. * XCHPFB : forces de contact en base B (lorsqu'on considere un
  140. * maillage de contact dans certaines liaisons)
  141. * KTOTXA/XB/VA/VB : Jacobienne par rapport au deplacement/vitesse en
  142. * base A/B (= contributions a dFnl/dX)
  143. *
  144. *
  145. * Segment des points de reference des modes (base A):
  146. * --------------------------------------------------
  147. SEGMENT MPREF
  148. INTEGER IPOREF(NPREF)
  149. ENDSEGMENT
  150. *
  151. * Segment des points en base B:
  152. * -----------------------------
  153. SEGMENT NCPR(NBPTS)
  154. * NCRP(#global) = #local dans XPTB (1er indice)
  155. *
  156. * Segment des parametres numeriques pour la continuation:
  157. * ------------------------------------------------------
  158. SEGMENT PARNUM
  159. CHARACTER*4 TYPS
  160. REAL*8 DS,DSMAX,DSMIN,ANGMIN,ANGMAX,ITERMOY,ISENS,TOLMIN
  161. REAL*8 PARINI,PARFIN
  162. INTEGER ITERMAX,NBPAS
  163. LOGICAL JANAL
  164. ENDSEGMENT
  165. *
  166. * Segment des resultats:
  167. * ---------------------
  168. SEGMENT PSORT
  169. REAL*8 QSAVE(NT1,NPAS),WSAVE(NPAS),LSAVE(2,2*NA1,NPAS)
  170. REAL*8 VSAVE(NPAS)
  171. LOGICAL ZSAVE(NPAS)
  172. CHARACTER*2 TYPBIF(NBIFU)
  173. REAL*8 QBIFU(NT1,NBIFU),WBIFU(NBIFU),WBIF2(NBIFU)
  174. REAL*8 QPSIR(NT1,NBIFU),QPSII(NT1,NBIFU)
  175. INTEGER CBIF
  176. ENDSEGMENT
  177. * QSAVE(i,j) = Q harmonique i au pas j
  178. * VSAVE(j) = parametre de continuation (si non w) au j-eme pas
  179. * ZSAVE(j) = stabilite au j-eme pas
  180. * LSAVE(1,j) : partie reelle de l'exposant de Floquet
  181. * LSAVE(2,j) : partie imaginaire de l'exposant de Floquet
  182. * TYPBIF = {LimitPoint, BranchPoint, NeimarkSacker, PeriodDoubling}
  183. * QBIFU,WBIFU : vecteur Q et w au point de bifurcation
  184. * WBIF2 : partie imaginaire de l'exposant de Floquet
  185. * QPSIR,QPSII : vecteur propre au point de bifurcation
  186.  
  187. * Segment des tableaux de travail:
  188. * -------------------------------
  189. SEGMENT MTEMP
  190. REAL*8 RW(NT1),A,T0(NT1+1),TP(NT1+1),AMPX,AUX
  191. REAL*8 T02(NT1+2), TP2(NT1+2)
  192. INTEGER IPIV(NT1),IPIV2(NT1+1),IPIV3(NT1+2)
  193. REAL*8 res
  194. REAL*8 RHS(NT1+1),Ja(NT1+1,NT1+1)
  195. REAL*8 QOLD(NT1),OMEGOLD
  196. REAL*8 MATJA(NT1+1,NT1+1),Rw2(NT1)
  197. REAL*8 Jaa(NT1+2,NT1+2),RHS2(NT1+2),Ra(NT1),VV,VVOLD
  198. ENDSEGMENT
  199. * Jacobiennes augmentees
  200. * Ja : [ RX Rw ; dX dw]
  201. * Jaa: [ RX Rw Ra; gx 0 0; dX dw da]
  202.  
  203. * SEGMENT NNNN
  204. * REAL*8 IGAM2(nl1,NPC2),DL2(nl1)
  205. * ENDSEGMENT
  206.  
  207. *************************** fin TMDYNC.INC *****************************
  208.  
  209. -INC PPARAM
  210. -INC CCOPTIO
  211.  
  212. REAL*8 ZERO,ONE,TWO
  213. PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0)
  214.  
  215. segini,MWORK
  216. * Variables generalisees: coefficients de Fourier et frequence
  217. MTQ=KTQ
  218. * Matrices XK,KASM,XM,GAM,IGAM,DL
  219. MTKAM=KTKAM
  220. * Deformees modales
  221. MTPHI = KTPHI
  222. NSB = XPHILB(/1)
  223. NPLSB = XPHILB(/2)
  224. NA2 = XPHILB(/3)
  225. IDIMB = XPHILB(/4)
  226. * Liaisons sur base A
  227. MTLIAA = KTLIAA
  228. NLIAA = IPALA(/1)
  229. * Liaisons sur base B
  230. MTLIAB = KTLIAB
  231. NLIAB = IPALB(/1)
  232. NIP = XABSCI(/2)
  233. NPLB = JPLIB(/1)
  234. * Forces externes
  235. MTFEX=KTFEX
  236. * Truc local base A
  237. LOCLFA = KOCLFA
  238. * Truc local base B
  239. LOCLB1 = KOCLB1
  240. * Inconnues sur un pas de temps (AFT)
  241. MTPAS=KTPAS
  242. RIGIDE =.FALSE.
  243. c Initialisation FFT
  244. IERRD=0
  245. NCOU=NFFT
  246. lenwrk = 2*NFFT
  247. lensav = NFFT + INT(LOG(REAL(NFFT))/LOG(TWO)) + 4
  248. segini,wfft51
  249.  
  250.  
  251. *------ Time-domain displacements and velocities: XT, VT ---------------
  252.  
  253. CALL HBMDVEC(NT,NHBM,NDDL,Q1,OMEG,V1)
  254. CALL IFT(Q1,NT,NDDL,NHBM,NFFT,GAM,XT)
  255. CALL IFT(V1,NT,NDDL,NHBM,NFFT,GAM,VT)
  256.  
  257.  
  258. *------ Time-domain nonlinear forces, tangent matrix: FT, KTOTBA -------
  259.  
  260. c >>> Loop over time steps >>>
  261. PDT = ONE/NFFT
  262. DO I = 1,NFFT
  263.  
  264. c initialisation des vecteurs du pas de temps
  265. DO L=1,NDDL
  266. XAUX(L,1) = XT(L,I)
  267. VAUX(L,1) = VT(L,I)
  268. FTOTA(L,1)=ZERO
  269. ENDDO
  270. DO JJ=1,NDDL
  271. DO IJ=1,NDDL
  272. KTOTXA(IJ,JJ) = ZERO
  273. KTOTVA(IJ,JJ) = ZERO
  274. ENDDO
  275. ENDDO
  276.  
  277. c Fnl(x,v) en base A (ddls modaux)
  278. IF (NLIAA.NE.0) THEN
  279. CALL D2VLFA(XAUX,VAUX,FTOTA,NDDL,IPALA,IPLIA,XPALA,XVALA,
  280. & NLIAA,PDT,ZERO,I,1,FINERT,1,FTEST,
  281. & KTOTXA,KTOTVA,.TRUE.)
  282. ENDIF
  283.  
  284. c Fnl(x,-) en base B (ddls physique)
  285. IF (NLIAB.NE.0) THEN
  286. CALL D2VLFB(XAUX,VAUX,FTOTA,NDDL,IPALB,IPLIB,XPALB,XVALB,NLIAB,
  287. & XPHILB,JPLIB,NPLB,IDIMB,FTOTB,FTOTBA,XPTB,PDT,ZERO,I,IBASB,IPLSB,
  288. & INMSB,IORSB,NSB,NPLSB,NA2,1,FEXPSM,NPC1,IERRD,FTEST2,XABSCI,
  289. & XORDON,NIP,FEXB,RIGIDE,IAROTA,XCHPFB,
  290. & KTOTXA,KTOTVA,KTOTXB,KTOTVB,.TRUE.)
  291. ENDIF
  292.  
  293. c stockage de la force + JAC au pas de temps I
  294. DO J=1,NDDL
  295. FT(J,I) = FTOTA(J,1)
  296. DO K = 1,NDDL
  297. MATJX(J,K,I) = KTOTXA(J,K)
  298. MATJV(J,K,I) = KTOTVA(J,K)
  299. ENDDO
  300. ENDDO
  301.  
  302. ENDDO
  303. c <<< end of Loop over time steps <<<
  304.  
  305.  
  306. c------ Fourier coefficients ------------------
  307. c Initialize the WWSAVE array.
  308. CALL rfft1i (NFFT,wwsave,lensav,ier)
  309. IF (IERR.ne.0) RETURN
  310.  
  311. c------ Fourier coefficients of nonlinear forces: FNL ------------------
  312. * CALL DFT(FT,NDDL,NHBM,NFFT,IGAM,FNL)
  313. DO IM = 1,NDDL
  314. DO KI = 1,NFFT
  315. XX51(KI) = FT(IM,KI)
  316. ENDDO
  317. CALL rfft1f(NFFT,1,XX51,NFFT,wwsave,lensav,work,lenwrk,ier)
  318. IA = 0
  319. DO KI = IM,NT-(NDDL-IM),NDDL
  320. IA = IA+1
  321. FNL(KI) = XX51(IA)
  322. ENDDO
  323. ENDDO
  324.  
  325. c------ Fourier coefficients of nonlinear tangent matrix: dFNLdX -------
  326.  
  327. c >>> Loop over modes I,J >>>
  328. DO I = 1,NDDL
  329. DO J = 1,NDDL
  330.  
  331. c Extract row vectors for the (i,j)-th derivative term
  332. DO IK = 1,NFFT
  333. VCX(IK) = MATJX(I,J,IK)
  334. VCV(IK) = MATJV(I,J,IK)
  335. ENDDO
  336. c First series of FFTs on the columns of diagonal dfdx matrix
  337. DO II = 1,NFFT
  338. DO KI = 1,NFFT
  339. XX51(KI) = ZERO
  340. XV51(KI) = ZERO
  341. ENDDO
  342. XX51(II) = VCX(II)
  343. XV51(II) = VCV(II)
  344. CALL rfft1f(NFFT,1,XX51,NFFT,wwsave,lensav,work,lenwrk,ier)
  345. CALL rfft1f(NFFT,1,XV51,NFFT,wwsave,lensav,work,lenwrk,ier)
  346. DO JK = 1,2*NHBM+1
  347. auxDL(JK,II) = NFFT*XX51(JK)
  348. auxVDL(JK,II) = NFFT*XV51(JK)
  349. ENDDO
  350. ENDDO
  351. c Second series of FFTs on the rows of auxDL
  352. DO II = 1,2*NHBM+1
  353. DO KI = 1,NFFT
  354. XX51(KI) = auxDL(II,KI)
  355. XV51(KI) = auxVDL(II,KI)
  356. ENDDO
  357. CALL rfft1f(NFFT,1,XX51,NFFT,wwsave,lensav,work,lenwrk,ier)
  358. CALL rfft1f(NFFT,1,XV51,NFFT,wwsave,lensav,work,lenwrk,ier)
  359. JFR(II,1) = XX51(1)
  360. JVFR(II,1) = XV51(1)
  361. DO JK = 2,2*NHBM+1
  362. JFR(II,JK) = XX51(JK)/2.
  363. JVFR(II,JK) = XV51(JK)/2.
  364. ENDDO
  365. ENDDO
  366.  
  367. c Assemble the total matrices JAC, JACV
  368. IL = 0
  369. DO NI = I,NDDL*(2*NHBM)+I,NDDL
  370. IL = IL + 1
  371. IC = 0
  372. DO NJ = J,NDDL*(2*NHBM)+J,NDDL
  373. IC = IC + 1
  374. JAC(NI,NJ) = JFR(IL,IC)
  375. JACV(NI,NJ)= JVFR(IL,IC)
  376. ENDDO
  377. ENDDO
  378.  
  379. ENDDO
  380. ENDDO
  381. c <<< fin de la boucle sur les modes I,J <<<
  382.  
  383. * Derivative of the velocity term
  384. CALL HBMDMAT(NT,NDDL,OMEG,JACV)
  385.  
  386. * Total nonlinear tangent stiffness: dFNLdX = JAC + w*JACV*Nabla
  387. DO I = 1,NT
  388. DO J = 1,NT
  389. JAC(I,J) = JAC(I,J) + JACV(I,J)
  390. ENDDO
  391. ENDDO
  392.  
  393.  
  394. c------ Menage ------------------
  395. segsup,MWORK
  396.  
  397. END
  398.  
  399.  
  400.  
  401.  

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