Télécharger la plaquette Cast3M


fr en

Télécharger hbm26.eso

Retour à la liste

Numérotation des lignes : 

hbm26

  1. C HBM26     SOURCE    CB215821  24/04/12    21:16:13     11897          
  2. *
  3. C DYNE26    SOURCE    BP208322  19/04/03    21:15:10     10174          
  4.       SUBROUTINE HBM26(IBAS,KTKAM,KTLIAB,KTPHI,IA1,IB,ICOMP,RIGIDE,ITKM)
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8(A-H,O-Z)
  7. *--------------------------------------------------------------------*
  8. *                                                                    *
  9. *     Operateur DYNE : algorithme de Fu - de Vogelaere               *
  10. *     ________________________________________________               *
  11. *                                                                    *
  12. *     Transpose l'information des objets de Castem2000 dans des      *
  13. *     tableaux de travail.                                           *
  14. *                                                                    *
  15. *     Parametres:                                                    *
  16. *                                                                    *
  17. * e   IBAS    Table representant une base modale                     *
  18. * es  KTKAM   Segment contenant les matrices XK, XASM et XM          *
  19. * es  KTPHI   Segment des deformees modales                          *
  20. * e   KTLIAB  Segment des liaisons sur base B                        *
  21. * es  IA1     Compteur                                               *
  22. * e   IB      Compteur de la sous base                               *
  23. * es  RIGIDE  Vrai si l'on a un corps rigide, faux sinon             *
  24. * e   ITKM    >0 si table RAIDEUR_ET_MASSE fournie                   *
  25. *                                                                    *
  26. *     Auteur, date de creation:                                      *
  27. *                                                                    *
  28. *     Lionel VIVAN, le 24 octobre 1989.                              *
  29. *                                                                    *
  30. *--------------------------------------------------------------------*
  31. -INC PPARAM
  32. -INC CCOPTIO
  33. -INC CCREEL
  34. -INC SMCHAML
  35. -INC SMELEME
  36. -INC SMMODEL
  37. *
  38. *-INC TMDYNC.INC
  39. ************************** debut TMDYNC.INC ****************************
  40.  
  41. *    TMDYNC : FUTUR INCLUDE POUR LES SEGMENTS DE L'OPERATEUR DYNC
  42. *    TODO   : a extraire dans un include des que stabilise
  43. *
  44. *     Segment des variables generalisees:
  45. *     -----------------------------------
  46.       SEGMENT MTQ
  47.          REAL*8 Q1(NT1)
  48.          REAL*8 OMEG,XPARA
  49.          REAL*8 JAC(NT1,NT1),ZZ(NT1,NT1),RX(NT1,NT1)
  50.          REAL*8 dX(NT1), dw, dv
  51.       ENDSEGMENT
  52. *     Q1   : vecteur des inconnues frequentielles de dimension (2h+1)*n
  53. *            Q1 = {q_0 q_c1 q_s1 ... q_sh}
  54. *            avec q_i vecteur de dimension n ou n=nombre de modes
  55. *     OMEG : frequence fondamentale de l'approximation
  56. *     XPARA: parametre de continuation (par defaut la frequence) 
  57. *            \in [PARINI,PARFIN]
  58. *     RX   : matrice jacobienne = ZZ + dFnl/dX
  59. *     JAC  : jacobienne des efforts non-lineaires  = dFnl/dX
  60. *     ZZ   : matrice dynamique associee aux matrices modales K, M et C 
  61. *            lineaires et constantes
  62. *     {dX,dw,(dv)} : vecteur tangent utilise pour la prediction
  63. *
  64. *
  65. *     Segment contenant les matrices XK, XASM et XM:
  66. *     ---------------------------------------------
  67.       SEGMENT MTKAM
  68.          REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  69.          REAL*8 GAM(NPC1,nl1),IGAM(nl1,NPC1),DL(nl1)
  70. *         REAL*8 GAMFIN(NPC2,nl1)
  71.       ENDSEGMENT
  72. *     XK,XASM et XM : matrices de raideur, amortissement et masse
  73. *     GAM et IGAM : matrices pour la FFT et son inverse
  74. *     GAMFIN      : 
  75. *
  76. *     Segment des deformees modales:
  77. *     ------------------------------
  78. *     (idem DYNE)
  79.       SEGMENT MTPHI
  80.          INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  81.          INTEGER IAROTA(NSB)
  82.          REAL*8  XPHILB(NSB,NPLSB,NA2,IDIMB)
  83.       ENDSEGMENT
  84. *
  85. *     Segment descriptif des liaisons en base A:
  86. *     ------------------------------------------
  87. *     (idem DYNE)
  88.       SEGMENT MTLIAA
  89.          INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
  90.          REAL*8  XPALA(NLIAA,NXPALA)
  91.       ENDSEGMENT
  92. *
  93. *     Segment descriptif des liaisons en base B:
  94. *     ------------------------------------------
  95. *     (idem DYNE)
  96.       SEGMENT MTLIAB
  97.          INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  98.          REAL*8  XPALB(NLIAB,NXPALB)
  99.          REAL*8  XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  100.       ENDSEGMENT
  101. *
  102. *     Segment representant les chargements exterieurs:
  103. *     -----------------------------------------------
  104.       SEGMENT MTFEX
  105.          REAL*8  FEXA(NT1)
  106.          REAL*8  FEXPSM(NPLB,NPC1,2,IDIMB)
  107.          INTEGER BAL
  108.       ENDSEGMENT
  109. *     FEXA : Vecteur des efforts ext. sous la forme de coefficients de
  110. *            Fourier et exprimes en base A
  111. *     FEXPSM: chargement/deplacement statique lie aux modes negliges
  112. *             (neglige aussi les Fnl). Dans DYNC toujours =0, cree pour
  113. *             compatibilite avec calcul des Fnl.
  114. *     BAL  : indique s'il s'agit d'un chargement de type balourd
  115. *            (cad proportionnel a OMEG**2)
  116. *
  117. *     Segment "local" pour DEVLFA:
  118. *     ----------------------------
  119.       SEGMENT LOCLFA
  120.          REAL*8 FTEST(NA1,4)
  121.       ENDSEGMENT
  122. *      
  123. *     Segment "local" pour DEVLB1:
  124. *     ----------------------------
  125.       SEGMENT LOCLB1
  126.          REAL*8 FTEST2(NPLB,6)
  127.       ENDSEGMENT
  128. *      
  129. *     Segment contenant les variables au cours d un pas de temps:
  130. *     ----------------------------------------------------------
  131.       SEGMENT MTPAS
  132.          REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
  133.          REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4)
  134.          REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
  135.          REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
  136.          REAL*8 KTOTXA(NA1,NA1),KTOTVA(NA1,NA1)
  137.          REAL*8 KTOTXB(NPLB,IDIMB,IDIMB), KTOTVB(NPLB,IDIMB,IDIMB) 
  138.       ENDSEGMENT
  139. *     FTOTA/B/BA : forces sur base A, B et B projetees sur A
  140. *     XPTB       : deplacement du point d'une liaison en base B
  141. *     XVALA/B    : grandeurs de la liaison en base A/B a stocker
  142. *     FEXB       : forces exterieures en base B (a priori uniquement 
  143. *                  pour les moments appliques aux rotations rigides ?)
  144. *     XCHPFB     : forces de contact en base B (lorsqu'on  considere un
  145. *                  maillage de contact dans certaines liaisons)
  146. *     KTOTXA/XB/VA/VB : Jacobienne par rapport au deplacement/vitesse en
  147. *                       base A/B (= contributions a dFnl/dX)
  148. *
  149. *
  150. *     Segment des points de reference des modes (base A):
  151. *     --------------------------------------------------
  152.       SEGMENT MPREF
  153.          INTEGER IPOREF(NPREF)
  154.       ENDSEGMENT
  155. *
  156. *     Segment des points en base B:
  157. *     -----------------------------
  158.       SEGMENT NCPR(XCOOR(/1)/(IDIM+1))
  159. *     NCRP(#global) = #local dans XPTB (1er indice)      
  160. *
  161. *     Segment des parametres numeriques pour la continuation: 
  162. *     ------------------------------------------------------
  163.       SEGMENT PARNUM
  164.          CHARACTER*4 TYPS
  165.          REAL*8   DS,DSMAX,DSMIN,ANGMIN,ANGMAX,ITERMOY,ISENS,TOLMIN
  166.          REAL*8   PARINI,PARFIN
  167.          INTEGER  ITERMAX,NBPAS
  168.          LOGICAL  JANAL
  169.       ENDSEGMENT
  170. *
  171. *     Segment des resultats:
  172. *     ---------------------
  173.       SEGMENT  PSORT
  174.         REAL*8 QSAVE(NT1,NPAS),WSAVE(NPAS),LSAVE(2,2*NA1,NPAS)
  175.         REAL*8 VSAVE(NPAS)
  176.         LOGICAL ZSAVE(NPAS)
  177.         CHARACTER*2 TYPBIF(NBIFU)
  178.         REAL*8 QBIFU(NT1,NBIFU),WBIFU(NBIFU),WBIF2(NBIFU)
  179.         REAL*8 QPSIR(NT1,NBIFU),QPSII(NT1,NBIFU)
  180.         INTEGER CBIF
  181.       ENDSEGMENT
  182. *     QSAVE(i,j) = Q harmonique i au pas j
  183. *     VSAVE(j)   = parametre de continuation (si non w) au j-eme pas
  184. *     ZSAVE(j)   = stabilite au j-eme pas
  185. *     LSAVE(1,j) : partie reelle de l'exposant de Floquet
  186. *     LSAVE(2,j) : partie imaginaire de l'exposant de Floquet
  187. *     TYPBIF = {LimitPoint, BranchPoint, NeimarkSacker, PeriodDoubling}
  188. *     QBIFU,WBIFU : vecteur Q et w au point de bifurcation
  189. *           WBIF2 : partie imaginaire de l'exposant de Floquet
  190. *     QPSIR,QPSII : vecteur propre au point de bifurcation
  191.  
  192. *     Segment des tableaux de travail:
  193. *     -------------------------------
  194.       SEGMENT  MTEMP
  195.          REAL*8 RW(NT1),A,T0(NT1+1),TP(NT1+1),AMPX,AUX
  196.          REAL*8 T02(NT1+2), TP2(NT1+2)
  197.          INTEGER IPIV(NT1),IPIV2(NT1+1),IPIV3(NT1+2)
  198.          REAL*8 res
  199.          REAL*8 RHS(NT1+1),Ja(NT1+1,NT1+1)
  200.          REAL*8 QOLD(NT1),OMEGOLD
  201.          REAL*8 MATJA(NT1+1,NT1+1),Rw2(NT1)
  202.          REAL*8 Jaa(NT1+2,NT1+2),RHS2(NT1+2),Ra(NT1),VV,VVOLD
  203.       ENDSEGMENT
  204. *     Jacobiennes augmentees 
  205. *     Ja : [ RX  Rw ; dX dw]
  206. *     Jaa: [ RX  Rw Ra; gx 0 0; dX dw da]
  207.  
  208. *      SEGMENT  NNNN
  209. *         REAL*8 IGAM2(nl1,NPC2),DL2(nl1)
  210. *      ENDSEGMENT
  211.  
  212. *************************** fin TMDYNC.INC *****************************
  213.  
  214. *
  215.       LOGICAL L0,L1,RIGIDE
  216.       CHARACTER*4 NOMTRI(6),NOMAXI(6),NOMPLA(3)
  217.       CHARACTER*8 CMOT,TYPRET,MORIGI,CHARRE
  218.       REAL*8      XAXROT(3),XROTA(2,3)
  219. *
  220. *     si IFOMOD = -1  :  modele PLAN
  221. *     si IFOMOD =  0  :  modele AXIS
  222. *     si IFOMOD =  1  :  modele FOUR
  223. *     si IFOMOD =  2  :  modele TRID
  224. *
  225. *     Les noms de composante sont
  226. *       - en modele PLAN   :  UX, UY, RT
  227. *       - en modele AXIS   :  UX, UY, RZ
  228. *       - en modele FOUR 1 :  UR, UZ, UT, RT
  229. *       - en modele TRID   :  UX, UY, UZ, RX, RY, RZ
  230. *
  231.       DATA NOMTRI/'UX  ','UY  ','UZ  ','RX  ','RY  ','RZ  '/
  232.       DATA NOMAXI/'UR  ','UT  ','UZ  ','RR  ','RT  ','RZ  '/
  233.       DATA NOMPLA/'UX  ','UY  ','RZ  '/
  234. *
  235.       MTKAM  = KTKAM
  236.       MTPHI  = KTPHI
  237.       MTLIAB = KTLIAB
  238. *
  239.       NLIAB  = IPALB(/1)
  240.       NPLB   = JPLIB(/1)
  241.       NSB    = XPHILB(/1)
  242.       NPLSB  = XPHILB(/2)
  243.       NA2    = XPHILB(/3)
  244.       IDIMB  = XPHILB(/4)
  245.       DEUXPI = 2.D0 * XPI
  246. *
  247.       IORSB(IB) = IA1 + 1
  248.       IAROTA(IB) = 0
  249.       IROT = 0
  250.       IN = 0
  251.  
  252. ************************************************************************
  253. *     table BASE_MODALE
  254. ************************************************************************
  255.  
  256.  10   CONTINUE
  257.       IN = IN + 1
  258.       TYPRET = ' '
  259.       CALL ACCTAB(IBAS,'ENTIER',IN,X0,' ',L0,IP0,
  260.      &                   TYPRET,I1,X1,CHARRE,L1,IBAMOD)
  261.       IF (IERR.NE.0) RETURN
  262. *    -on a bien un objet de type table
  263.       IF (IBAMOD.NE.0) THEN
  264.        IF (TYPRET.EQ.'TABLE   ') THEN
  265.  
  266.          IA1 = IA1 + 1
  267.  
  268. *        remplissage de XM et XK diagonale depuis la table BASE_MODALE
  269. *        sauf si deja fait car on a une table RAIDEUR_ET_MASSE !
  270.          IF (ITKM.LE.0) THEN
  271.            CALL ACCTAB(IBAMOD,'MOT',I0,X0,'MASSE_GENERALISEE',L0,IP0,
  272.      &                        'FLOTTANT',I1,XMASSE,' ',L1,IP1)
  273.            IF (IERR.NE.0) RETURN
  274.            XM(IA1,1) = XMASSE
  275.            CALL ACCTAB(IBAMOD,'MOT',I0,X0,'FREQUENCE',L0,IP0,
  276.      &                        'FLOTTANT',I1,XFREQ,' ',L1,IP1)
  277.            IF (IERR.NE.0) RETURN
  278.            OMEGA = XFREQ * DEUXPI
  279.            XK(IA1,1) = XMASSE * OMEGA * OMEGA
  280.            IF (IIMPI.EQ.333) THEN
  281.               WRITE(IOIMP,*)'HBM26 :   XM(',IA1,') =',XMASSE
  282.               WRITE(IOIMP,*)'HBM26 :   XK(',IA1,') =',XK(IA1,1)
  283.            ENDIF
  284.          ENDIF
  285.  
  286. *    si liaison_B existe, remplissage de IPLSB, XPHILB, IAROTA, INMSB...
  287.          IF (NLIAB.NE.0) THEN
  288.             CALL ACCTAB(IBAMOD,'MOT',I0,X0,'DEFORMEE_MODALE',L0,IP0,
  289.      &                         'CHPOINT',I1,X1,' ',L1,ICDM)
  290.             IF (IERR.NE.0) RETURN
  291.             CALL ACTOBJ('CHPOINT',ICDM,1)
  292.  
  293.             DO 12 ID = 1,IDIMB
  294.                IF (IFOUR.EQ.0 .OR. IFOUR.EQ.1) THEN
  295.                   CMOT = NOMAXI(ID)
  296.                ELSE
  297.                   IF (IFOMOD.EQ.-1) THEN
  298.                       CMOT = NOMPLA(ID)
  299.                   ELSE
  300.                       CMOT = NOMTRI(ID)
  301.                   ENDIF
  302.                ENDIF
  303.                IF (IIMPI.EQ.333)
  304.      &         WRITE(IOIMP,*)'HBM26 :  composante a extraire :',CMOT
  305.                ICOMP  = 0
  306.                DO 14 IP = 1,NPLB
  307.                   IPOINT = JPLIB(IP)
  308. *      On extrait du chpoint ICDM au point IPOINT de composante CMOT
  309.                   CALL EXTRA9(ICDM,IPOINT,CMOT,0,.FALSE.,XVAL,IRET)
  310.                   ICOMP = ICOMP + 1
  311. *          on ajuste la taille si necessaire
  312.                   IF(ICOMP.GT.NPLSB) THEN
  313.                       NPLSB=ICOMP
  314.                       SEGADJ MTPHI
  315.                   ENDIF
  316.                   IPLSB(IP) = ICOMP
  317. * suite a la modif dans extra9, car on attribue une valeur meme
  318. * si le point n'existe pas dans le chpoint
  319.                   IF (XVAL.NE.0.) THEN
  320.                     IF ((IBASB(IP).NE.0).AND.(IBASB(IP).NE.IB)) THEN
  321.                       call erreur (783)
  322.                       RETURN
  323.                     ENDIF
  324.                     IBASB(IP) = IB
  325.                   ELSEIF ((IB.EQ.NSB).AND.(IBASB(IP).EQ.0)) THEN
  326.                     IBASB(IP) = IB
  327.                   ENDIF
  328.                   XPHILB(IB,ICOMP,IN,ID) = XVAL
  329.       IF (IIMPI.EQ.333) THEN
  330.         WRITE(IOIMP,*)'HBM26 :   IPLSB(',IP,') =',IPLSB(IP)
  331.         WRITE(IOIMP,*)'HBM26 :   IBASB(',IP,') =',IBASB(IP)
  332.         XVA2 = XPHILB(IB,ICOMP,IN,ID)
  333.         WRITE(IOIMP,*)'HBM26 :   XPHILB(',IB,ICOMP,IN,ID,') =',XVA2
  334.       ENDIF
  335.  14            CONTINUE
  336.  12         CONTINUE
  337.          ENDIF
  338.  
  339. c *        Prise en compte d'un mode de rotation de corps rigide
  340.           MORIGI = ' '
  341.           CALL ACCTAB(IBAMOD,'MOT',I0,X0,'CORPS_RIGIDE',L0,IP0,
  342.      &                      MORIGI,I1,X1,CMOT,L1,IP1)
  343.           IF (IERR.NE.0) RETURN
  344.           IF (MORIGI.EQ.'MOT') THEN
  345.              IF (CMOT(1:4).EQ.'VRAI') THEN
  346.                  CALL ACCTAB(IBAMOD,'MOT',I0,X0,'CENTRE_DE_GRAVITE',
  347.      &          L0,IP0,'POINT',I1,X1,'  ',L1,ICDG)
  348.                  IF (IERR.NE.0) RETURN
  349.                  IAROTA(IB)=IA1
  350.                  IROT = IN
  351.              ENDIF
  352.           ENDIF
  353.          GOTO 10
  354.        ELSE
  355.          CALL ERREUR(491)
  356.          RETURN
  357.        ENDIF
  358.       ENDIF
  359. *    -fin du cas ou on a bien un objet de type table
  360.       INMSB(IB) = IN - 1
  361. *       
  362. ************************************************************************
  363. *     Remplissage des fausses deformees modales de rotations
  364. ************************************************************************
  365. *
  366. *50   continue
  367.       IF (IAROTA(IB).NE.0) THEN
  368.          RIGIDE = .TRUE.
  369.          MERR = 0
  370.          NPLUS = IN + 1
  371.          IF (NPLUS.GT.NA2) THEN
  372. * On reajuste le dimension NA2 de XPHILB
  373.                  NA2 = NPLUS
  374.                  SEGADJ MTPHI
  375.          ENDIF
  376.          DO 18 IP=1,NPLB
  377.              IPOINT=JPLIB(IP)
  378.              IPOS=IPLSB(IP)
  379.              IBBAS= IBASB(IP)
  380.              IF (IBBAS.EQ.IB) THEN
  381.                 DO 20 ID=(IDIM+1),IDIMB
  382.                   XAXROT(ID-IDIM) = XPHILB(IB,IPOS,IROT,ID)
  383.  20             CONTINUE
  384. * En tridimensionnel l'axe de rotation est le vecteur propre de rotation
  385. * On norme l axe du plan de rotation
  386.                 CALL DYNE41(XAXROT,MERR,IDIM)
  387. * En bidimensionnel l'axe de rotation est fixe
  388. * Calcul des fausses deformees modales de rotation
  389.                 CALL DYNE42(XROTA,XAXROT,IPOINT,ICDG,IDIMB,MERR)
  390.              DO 22 ID =1,IDIMB
  391.                  XPHILB(IB,IPOS,IN,ID)  = XROTA(1,ID)
  392.                  XPHILB(IB,IPOS,IN+1,ID)= XROTA(2,ID)
  393.  22          CONTINUE
  394.              ENDIF
  395.  18       CONTINUE
  396.       ENDIF
  397.  
  398.       IF (IIMPI.EQ.333) THEN
  399.         WRITE(IOIMP,*)'HBM26 :   INMSB(',IB,') =',INMSB(IB)
  400.         WRITE(IOIMP,*)'HBM26 :   IORSB(',IB,') =',IORSB(IB)
  401.         WRITE(IOIMP,*)'HBM26 :   IAROTA(',IB,') =',IAROTA(IB)
  402.       ENDIF      
  403.       END
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  

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