Télécharger la plaquette Cast3M


fr en

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