Télécharger la plaquette Cast3M


fr en

Télécharger hbmfnl.eso

Retour à la liste

Numérotation des lignes : 

hbmfnl

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

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