Télécharger la plaquette Cast3M


fr en

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