Télécharger la plaquette Cast3M


fr en

Télécharger d2vlfa.eso

Retour à la liste

Numérotation des lignes : 

d2vlfa

  1. C D2VLFA    SOURCE    BP208322  20/09/18    21:15:17     10718          
  2. c
  3.       SUBROUTINE D2VLFA(Q1,Q2,FTOTA,NA1,IPALA,IPLIA,XPALA,XVALA,
  4.      &                  NLIAA,PDT,T,NPAS,IND,FINERT,IVINIT,FTEST,
  5.      &                  KTOTXA,KTOTVA,GETJAC)
  6.  
  7.       IMPLICIT INTEGER(I-N)
  8.       IMPLICIT REAL*8(A-H,O-Z)
  9.  
  10. *--------------------------------------------------------------------*
  11. *     Operateur DYN* : Calcul des forces de choc base A              *
  12. *--------------------------------------------------------------------*
  13. *                                                                    *
  14. *     Parametres:                                                    *
  15. *                                                                    *
  16. * e   Q1(.,.) Vecteur des deplacements generalises.                  *
  17. * e   Q2(.,.) Vecteur des vitesses generalisees.                     *
  18. * es  FTOTA   Forces exterieures totalisees sur la base A.           *
  19. * es  KTOTXA  Matrice tangente,X des efforts non-lineaires (base A)  *
  20. * es  KTOTVA  Matrice tangente,V des efforts non-lineaires (base A)  *
  21. * e   NA1     Nombre total d'inconnues en base A.                    *
  22. * e   IPALA   Renseigne sur le type de la liaison.                   *
  23. * e   IPLIA   Tableau contenant les numeros "DYNE" de la liaison.    *
  24. * e   XPALA   Tableau contenant les parametres de la liaison.        *
  25. * es  XVALA   Tableau contenant les variables internes des liaisons  *
  26. * e   NLIAA   Nombre de liaisons sur la base A.                      *
  27. * e   PDT     pas de temps                                           *
  28. * e   T       temps                                                  *
  29. * e   NPAS    Numero du pas de temps                                 *
  30. * e   IND     Indice du demi-pas de temps                            *
  31. *             = 2 si 1er demi-pas  = 1 si 2eme demi-pas              *
  32. * es  FINERT  Forces d'inertie                                       *
  33. * e   IVINIT  =1 si vitesses initiales, =0 sinon                     *
  34. * e   GETJAC  .TRUE. si on doit calculer la jacobienne KTOTXA,KTOTVA *
  35. *                                                                    *
  36. *   remarque:                                                        *
  37. *   =========                                                        *
  38. *   Si jeu negatif (cas particulier de la base A ou il n'y a pas de  *
  39. *      normale), on renverse les variables avec XNORM.               *
  40. *                                                                    *
  41. *--------------------------------------------------------------------*
  42. *
  43.       INTEGER IPALA(NLIAA,*),IPLIA(NLIAA,*)
  44.       REAL*8  XPALA(NLIAA,*),Q1(NA1,*),Q2(NA1,*),FTOTA(NA1,*)
  45.       REAL*8  XVALA(NLIAA,4,*),FINERT(NA1,*)
  46.       PARAMETER (XZERO = 0.D0, XONE=1.D0)
  47.       REAL*8  FTest(nA1,4)
  48. cbp, supprime le 2020-08-07          REAL*8  FTOTA0(NA1,4)
  49.       REAL*8  KTOTXA(NA1,*),KTOTVA(NA1,*) 
  50.       LOGICAL GETJAC
  51. *
  52.       XFIN = 0.D0
  53.       PDTS2 = PDT
  54.  
  55. *--------------------------------------------------------------------*
  56. *     BOUCLE SUR LES LIAISONS
  57. *--------------------------------------------------------------------*
  58.  
  59.       DO 10 I = 1,NLIAA
  60.  
  61.          ITYP = IPALA(I,1)
  62.          icond= IPALA(I,2)
  63.          iannul= 0
  64.  
  65.          IF (ICOND .NE. 1 ) GOTO 199
  66. *        CAS DES LIAISONS CONDITIONNELLES : 
  67. *>>>>>>> BOUCLE SUR LES LIAISONS "TESTS"     <<<<<<<<<<<<<<<<<<<<<<<<<<
  68.          DO 101 j = 4,20
  69.  
  70.            jliai = ipala(i,j)
  71.            jpliai = abs ( jliai)
  72.            if (jliai.EQ.0) goto 101
  73.  
  74.            jtyp = ipala(jpliai,1)
  75.            do 102 jk = 1,4
  76.            do 103 ik = 1,nA1
  77.                    ftest(ik,jk) = 0d0
  78. cbp, supprime le 2020-08-07    ftota0 (ik,jk) = ftota(ik,jk)
  79.  103       continue
  80.  102       continue
  81.  
  82. *      >>> CALCUL DES FORCES DES LIAISONS TEST <<<
  83. *
  84. * ------ choc elementaire POINT_PLAN sans amortissement
  85. *
  86.          IF (JTYP.EQ.1) THEN
  87.             XRAID = XPALA(jpliai,1)
  88.             XJEU  = XPALA(jpliai,2)
  89.             ETA = 0.D0
  90.             XNORM = XONE
  91.             IF (XJEU.LT.0D0) THEN
  92.                 XNORM = -XONE
  93.                 XJEU = -XJEU
  94.             ENDIF
  95.             INA1 = IPLIA(jpliai,1)
  96.             XDEP = XNORM*Q1(INA1,IND)
  97.             CALL DYCHEL(XDEP,XRAID,XJEU,ETA,XFL,DFDX,iannul)
  98.             XVALA(jpliai,IND,1) = XNORM*XFL
  99.             XVALA(jpliai,IND,4) = XNORM*XDEP
  100.             FTest(INA1,IND) =Ftest(INA1,IND) + XNORM*XFL
  101. *
  102. * ------ choc elementaire POINT_PLAN avec amortissement
  103. *
  104.          ELSE IF (JTYP.EQ.2) THEN
  105.             XRAID = XPALA(jpliai,1)
  106.             XJEU = XPALA(jpliai,2)
  107.             XAMO = XPALA(jpliai,3)
  108.             ETA = 0.D0
  109.             XNORM = XONE
  110.             IPERM = 0
  111.             IF (XJEU.LT.XZERO) THEN
  112.                 XNORM = -XONE
  113.                 XJEU = -XJEU
  114.             ENDIF
  115.             INA1 = IPLIA(jpliai,1)
  116.             XDEP = XNORM * Q1(INA1,IND)
  117. cbp,2020-09            IF (((NPAS.EQ.1).AND.(IND.EQ.3))) THEN
  118.                XVIT = XNORM * Q2(INA1,IND)
  119. cbp,2020-09            ELSE
  120. cbp,2020-09               IND2 = IND + 1
  121. cbp,2020-09               XDEPM1 =  XNORM * Q1(INA1,IND2)
  122. cbp,2020-09               XVIT = (XDEP - XDEPM1) / PDTS2
  123. cbp,2020-09            ENDIF
  124.             XVALA(jpliai,IND,3) = XNORM*XVIT
  125.             CALL DYCHAM(XDEP,XVIT,XRAID,XJEU,XAMO,ETA,
  126.      &                  XFL,DFDX,DFDV,IPERM,iannul)
  127. *           bp,2020 : pertinence des 2 lignes ci-dessous (cf. dycham) ? 
  128. *           -> a verifier + tard ...
  129. c             IF(ETA.EQ.0.)THEN
  130. c             IF (XDEP.GT.0.D0 .AND. XFL.GT.0.D0) XFL = 0.D0
  131. c             IF (XDEP.LT.0.D0 .AND. XFL.LT.0.D0) XFL = 0.D0
  132. c             ENDIF
  133.             XVALA(jpliai,IND,1) = XNORM*XFL
  134.             XVALA(jpliai,IND,4) = XNORM*XDEP
  135.             FTest(INA1,IND) =Ftest(INA1,IND) + XNORM*XFL
  136. *
  137. * ------ choc elementaire POINT_PLAN_FLUIDE
  138. *
  139.          ELSE IF (JTYP.EQ.3) THEN
  140.             XINER = XPALA(jpliai,1)
  141.             XCONV = XPALA(jpliai,2)
  142.             XVISC = XPALA(jpliai,3)
  143.             XPCEL = XPALA(jpliai,4)
  144.             XPCRA = XPALA(jpliai,5)
  145.             XJEU = XPALA(jpliai,6)
  146.             INA1 = IPLIA(jpliai,1)
  147.             XDEP = Q1(INA1,IND)
  148. cbp,2020-09            IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  149.                XVIT = Q2(INA1,IND)
  150. cbp,2020-09            ELSE
  151. cbp,2020-09               IND2 = IND + 1
  152. cbp,2020-09               XDEPM1 = Q1(INA1,IND2)
  153. cbp,2020-09               XVIT = (XDEP - XDEPM1) / PDTS2
  154. cbp,2020-09            ENDIF
  155.             IF (XJEU.GT.0D0) THEN
  156.                XDH= XJEU - XDEP
  157.                XNORM = XONE
  158.             ELSE
  159.                XDH= XDEP - XJEU
  160.                XNORM = -XONE
  161.             ENDIF
  162. *   Calcul de la masse ajoutee
  163.             XXIN = -XINER / XDH
  164.             FINERT(INA1,IND) = FINERT(INA1,IND) + XXIN
  165. *   Calcul de la force de convection
  166.             CALL DYFCON(XDH,XDEP,XVIT,XJEU,XCONV,XFCO,iannul)
  167. *   Calcul de la force de viscosite
  168.             CALL DYFVIS(XDH,XDEP,XVIT,XJEU,XVISC,XFVI,iannul)
  169. *   Calcul de la force de perte de charge
  170.             CALL DYFPDC(XDH,XDEP,XVIT,XJEU,XPCEL,XPCRA,XFPE,iannul)
  171.             XFL = (XFCO* XNORM )  + XFVI + XFPE
  172.             XVALA(jpliai,IND,1) = XDEP
  173.             XVALA(jpliai,IND,2) = XVIT
  174.             XVALA(jpliai,IND,3) = XXIN
  175.             XVALA(jpliai,IND,4) = XFCO*XNORM
  176.             XVALA(jpliai,IND,5) = XFVI
  177.             XVALA(jpliai,IND,6) = XFPE
  178.             FTest(INA1,IND) =Ftest(iNA1,IND) + XFL
  179. *
  180. * ------ force elementaire de COUPLAGE EN VITESSE
  181. *
  182.          ELSE IF (JTYP.EQ.4) THEN
  183.             INA1 = IPLIA(jpliai,1)
  184.             INA2 = IPLIA(jpliai,2)
  185.             XDEP = Q1(INA2,IND)
  186. cbp,2020-09            IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  187.                XVIT = Q2(INA2,IND)
  188. cbp,2020-09            ELSE
  189. cbp,2020-09               IND2 = IND + 1
  190. cbp,2020-09               XDEPM1 = Q1(INA2,IND2)
  191. cbp,2020-09               XVIT = (XDEP - XDEPM1) / PDTS2
  192. cbp,2020-09            ENDIF
  193.             XCPLGE = XPALA(jpliai,1)
  194.             CALL DYCPLV(XVIT,XCPLGE,XFL,iannul)
  195.             XVALA(jpliai,IND,3) = XVIT
  196.             XVALA(jpliai,IND,1) = XFL
  197.             XVALA(jpliai,IND,4) = XDEP
  198.             FTest(INA1,IND) =Ftest(INA1,IND) + XFL
  199. *
  200. * ------ force elementaire de COUPLAGE EN DEPLACEMENT
  201. *
  202.          ELSE IF (JTYP.EQ.5) THEN
  203.             INA1 = IPLIA(jpliai,1)
  204.             INA2 = IPLIA(jpliai,2)
  205.             XDEP = Q1(INA2,IND)
  206.             XCPLGE = XPALA(jpliai,1)
  207.             jfonct = ipala(jpliai,3)
  208.  
  209. cbp       - calcul eventuel d'un produit de convolution :
  210. c           XDEP=\int_0^T h(\tau)*Qj(t-\tau) d\tau
  211.             if(jfonct.eq.100) then
  212.                IP1=IPALA(jpliai,4)
  213. c              t_{n+1} ou t_{n}
  214. c                if (IND.eq.1.or.IND.eq.3) then 
  215. c              t_{n+1} ou t_{0} pour diff centreee ?
  216. c                if (IND.eq.1.or.IND.eq.2) then 
  217.                  IP2=IPALA(jpliai,5)
  218. c c              t_{n+1/2} ou t_{n-1/2} (seulement pour devogelaere)
  219. c                else
  220. c                  IP2=IPALA(jpliai,6)
  221. c                endif
  222.                CALL DYCPL1(IP1,IP2,XDEP,NPAS,PDT,XCONV)
  223.                XFL=XCPLGE*XCONV
  224.  
  225. cbp       - calcul d'un produit de convolution GRANGER_PAIDOUSSIS
  226.             elseif(jfonct.eq.101) then
  227.                IP1=IPALA(jpliai,4)
  228.                IP2=IPALA(jpliai,5)
  229.                IP3=IPALA(jpliai,6)
  230.                VSD=XPALA(jpliai,2)
  231.                XA0=XPALA(jpliai,3)
  232.                CALL DYCPL2(IP1,IP2,IP3,VSD,XA0,XDEP,NPAS,PDT,XCONV)
  233.                XFL=XCPLGE*XCONV               
  234.  
  235. cbp       - calcul du produit par une fonction trigo cos(q)
  236. c           rem : on pourrait aussi prevoir cos(coef2*q) avec coef2=XPALA(jpliai,2)
  237.             elseif(jfonct.eq.1) then
  238.                 XFL = XCPLGE * COS(XDEP)
  239.             elseif(jfonct.eq.2) then
  240.                 XFL = XCPLGE * SIN(XDEP)  
  241.  
  242. cbp       - calcul du produit par une fonction temporelle cos(wt)*q
  243.             elseif(jfonct.ge.11.and.jfonct.le.12) then
  244.                 XFREQ  = XPALA(jpliai,2)
  245.                 XTIME = DBLE(NPAS-1)*PDT
  246.                 if(IND.EQ.2) XTIME=XTIME-PDTS2
  247.                 if(jfonct.eq.1) XFONCT = COS(XFREQ*XTIME)
  248.                 if(jfonct.eq.2) XFONCT = SIN(XFREQ*XTIME)
  249.                 XFL = XCPLGE * XDEP * XFONCT
  250.  
  251. c         - simple raideur (ou raideur en puissance)
  252.             else
  253.                   Xexpo = XPALA(jpliai,2)
  254.                   if(Xexpo.eq.XONE) then
  255.                     XFL = XCPLGE * XDEP
  256.                   elseif(Xexpo.eq.3.D0) then
  257.                     XFL = XCPLGE * (XDEP**3)
  258.                     IF(GETJAC) DFDX = Xexpo*XCPLGE*(XDEP**2) 
  259.                   else
  260.                     XFL = XCPLGE * (XDEP**Xexpo)
  261.                   endif
  262.             endif
  263. cbp, ici    prise en compte de la possibilite de liaison conditionnelle 
  264.             IF(iannul.ne.0) XFL=0.D0           
  265.             XVALA(jpliai,IND,1) = XFL
  266.             XVALA(jpliai,IND,4) = XDEP
  267.             FTest(INA1,IND) =Ftest(INA1,IND) + XFL
  268. *
  269. * ------ force elementaire de type POLYNOMIALE
  270. *
  271.          ELSE IF (JTYP.EQ.6) THEN
  272. *           nombre de modes "origine"
  273.             NMOD = IPALA(I,2)
  274.             CALL D2POL1(Q1,Q2,NA1,IPLIA,XPALA,XVALA,NLIAA,IND,PDT,
  275.      &                  NPAS,jpliai,NMOD,FTest,IVINIT,iannul)
  276. *
  277. * ------ choc elementaire ...
  278. *
  279. *        ELSE IF (JTYP.EQ.  ) THEN
  280. *              .......
  281. *              .......
  282. *
  283.          ENDIF
  284.  
  285.  
  286. *       >>> TEST SUR j POUR VOIR SI ELLE ANNULE LA LIAISON I <<<
  287.            xff = 0d0
  288.            do 104 ik = 1,na1
  289.               do 105 jk = 1,4
  290.                    xff = xff + ( ftest(ik,jk) ** 2)
  291.  105          continue
  292.  104       continue
  293.            xff = xff ** .5
  294.            if (   ((xff .le. 1e-20 ) .and. ( jliai .gt. 0) )
  295.      &       .OR. ((xff .gt. 1e-20 ) .and. ( jliai .lt. 0) )  )
  296.      &        then
  297.               iannul = 1
  298.            endif
  299.  
  300.  
  301.  101     CONTINUE
  302. *>>>>>>> FIN DE BOUCLE SUR LES LIAISONS "TESTS"     <<<<<<<<<<<<<<<<<<<<
  303.  
  304.  199     CONTINUE
  305. *>>>>>>> CALCUL DES FORCES DE LIAISONS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  306.  
  307. * ------ choc elementaire POINT_PLAN sans amortissement
  308. *
  309.          IF (ITYP.EQ.1) THEN
  310.             XRAID = XPALA(I,1)
  311.             XJEU = XPALA(I,2)
  312.             ETA = 0.D0
  313.             XNORM = XONE
  314.             IF (XJEU.LT.0D0) THEN
  315.                 XNORM = -XONE
  316.                 XJEU = -XJEU
  317.             ENDIF
  318.             INA1 = IPLIA(I,1)
  319.             XDEP = XNORM*Q1(INA1,IND)
  320.             CALL DYCHEL(XDEP,XRAID,XJEU,ETA,XFL,DFDX,iannul)
  321.             XVALA(I,IND,1) = XNORM*XFL
  322.             XVALA(I,IND,4) = XNORM*XDEP
  323.             FTOTA(INA1,IND) = FTOTA(INA1,IND) + XNORM*XFL
  324.             IF (GETJAC) THEN
  325.               KTOTXA(INA1,INA1)  = KTOTXA(INA1,INA1) + DFDX
  326.             ENDIF
  327.  
  328. *
  329. * ------ choc elementaire POINT_PLAN avec amortissement
  330. *
  331.          ELSE IF (ITYP.EQ.2) THEN
  332.             XRAID = XPALA(I,1)
  333.             XJEU = XPALA(I,2)
  334.             XAMO = XPALA(I,3)
  335.             ETA = 0.D0
  336.             XNORM = XONE
  337.             IPERM = 0
  338.             IF (XJEU.LT.XZERO) THEN
  339.                 XNORM = -XONE
  340.                 XJEU = -XJEU
  341.             ENDIF
  342.             INA1 = IPLIA(I,1)
  343.             XDEP = XNORM * Q1(INA1,IND)
  344. cbp,2020-09            IF (((NPAS.EQ.1).AND.(IND.EQ.3))) THEN
  345.                XVIT = XNORM * Q2(INA1,IND)
  346. cbp,2020-09            ELSE
  347. cbp,2020-09               IND2 = IND + 1
  348. cbp,2020-09               XDEPM1 =  XNORM * Q1(INA1,IND2)
  349. cbp,2020-09               XVIT = (XDEP - XDEPM1) / PDTS2
  350. cbp,2020-09            ENDIF
  351.             XVALA(I,IND,3) = XNORM*XVIT
  352.             CALL DYCHAM(XDEP,XVIT,XRAID,XJEU,XAMO,ETA,
  353.      &                  XFL,DFDX,DFDV,IPERM,iannul)
  354. *           bp,2020 : pertinence des 2 lignes ci-dessous (cf. dycham) ? 
  355. *           -> on commente... a verifier + tard
  356. *            IF (XDEP.GT.0.D0 .AND. XFL.GT.0.D0) XFL = 0.D0
  357. *            IF (XDEP.LT.0.D0 .AND. XFL.LT.0.D0) XFL = 0.D0
  358.             XVALA(I,IND,1) = XNORM*XFL
  359.             XVALA(I,IND,4) = XNORM*XDEP
  360.             FTOTA(INA1,IND) = FTOTA(INA1,IND) + XNORM*XFL
  361.             IF (GETJAC) THEN
  362.               KTOTXA(INA1,INA1)  = KTOTXA(INA1,INA1) + DFDX
  363.               KTOTVA(INA1,INA1)  = KTOTVA(INA1,INA1) + DFDV
  364.             ENDIF
  365. *
  366. * ------ choc elementaire POINT_PLAN_FLUIDE
  367. *
  368.          ELSE IF (ITYP.EQ.3) THEN
  369.             XINER = XPALA(I,1)
  370.             XCONV = XPALA(I,2)
  371.             XVISC = XPALA(I,3)
  372.             XPCEL = XPALA(I,4)
  373.             XPCRA = XPALA(I,5)
  374.             XJEU = XPALA(I,6)
  375.             INA1 = IPLIA(I,1)
  376.             XDEP = Q1(INA1,IND)
  377. *
  378. cbp,2020-09            IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  379.                XVIT = Q2(INA1,IND)
  380. cbp,2020-09            ELSE
  381. cbp,2020-09               IND2 = IND + 1
  382. cbp,2020-09               XDEPM1 = Q1(INA1,IND2)
  383. cbp,2020-09               XVIT = (XDEP - XDEPM1) / PDTS2
  384. cbp,2020-09            ENDIF
  385.             IF (XJEU.GT.0D0) THEN
  386.                XDH= XJEU - XDEP
  387.                XNORM = XONE
  388.             ELSE
  389.                XDH= XDEP - XJEU
  390.                XNORM = -XONE
  391.             ENDIF
  392. *   Calcul de la masse ajoutee
  393.             XXIN = -XINER / XDH
  394.             FINERT(INA1,IND) = FINERT(INA1,IND) + XXIN
  395. *   Calcul de la force de convection
  396.             CALL DYFCON(XDH,XDEP,XVIT,XJEU,XCONV,XFCO,iannul)
  397. *   Calcul de la force de viscosite
  398.             CALL DYFVIS(XDH,XDEP,XVIT,XJEU,XVISC,XFVI,iannul)
  399. *   Calcul de la force de perte de charge
  400.             CALL DYFPDC(XDH,XDEP,XVIT,XJEU,XPCEL,XPCRA,XFPE,iannul)
  401.             XFL = (XFCO* XNORM )  + XFVI + XFPE
  402.             XVALA(I,IND,1) = XDEP
  403.             XVALA(I,IND,2) = XVIT
  404.             XVALA(I,IND,3) = XXIN
  405.             XVALA(I,IND,4) = XFCO*XNORM
  406.             XVALA(I,IND,5) = XFVI
  407.             XVALA(I,IND,6) = XFPE
  408.             FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL
  409. * TODO :
  410. *           IF (GETJAC) THEN
  411. *             KTOTXA(INA1,INA1)=KTOTXA(INA1,INA1)+DFDXc*XNORM+DFDXv+DFDXp
  412. *             KTOTVA(INA1,INA1)=KTOTVA(INA1,INA1)+DFDVc*XNORM+DFDVv+DFDVp
  413. *           ENDIF
  414. *            
  415. *
  416. * ------ force elementaire de COUPLAGE EN VITESSE
  417. *
  418.          ELSE IF (ITYP.EQ.4) THEN
  419.             INA1 = IPLIA(I,1)
  420.             INA2 = IPLIA(I,2)
  421.             XDEP = Q1(INA2,IND)
  422. cbp,2020-09            IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  423.                XVIT = Q2(INA2,IND)
  424. cbp,2020-09            ELSE
  425. cbp,2020-09               IND2 = IND + 1
  426. cbp,2020-09               XDEPM1 = Q1(INA2,IND2)
  427. cbp,2020-09               XVIT = (XDEP - XDEPM1) / PDTS2
  428. cbp,2020-09            ENDIF
  429.             XCPLGE = XPALA(I,1)
  430.             CALL DYCPLV(XVIT,XCPLGE,XFL,iannul)
  431. c             WRITE(*,*) 't=',(NPAS*PDT),'dx/dt=',XVIT,'F2=',XFL
  432.             XVALA(I,IND,3) = XVIT
  433.             XVALA(I,IND,1) = XFL
  434.             XVALA(I,IND,4) = XDEP
  435.             FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL
  436.             IF (GETJAC) THEN
  437.             IF(IANNUL.EQ.0) KTOTVA(INA1,INA2)=KTOTVA(INA1,INA2)+XCPLGE
  438.             ENDIF
  439. *
  440. * ------ force elementaire de COUPLAGE EN DEPLACEMENT
  441. *
  442.          ELSE IF (ITYP.EQ.5) THEN
  443.             INA1 = IPLIA(I,1)
  444.             INA2 = IPLIA(I,2)
  445.             XDEP = Q1(INA2,IND)
  446.             XCPLGE = XPALA(I,1)
  447.             jfonct = ipala(I,3)
  448. c           avant d'avoir defini DFDX pour toutes les liaisons, mise a 0
  449.             DFDX=0.d0
  450.  
  451. cbp         calcul eventuel d'un produit de convolution :
  452. c           XDEP=\int_0^T h(\tau)*Qj(t-\tau) d\tau
  453.             if(jfonct.eq.100) then
  454.                IP1=IPALA(I,4)
  455. c              t_{n+1} ou t_{n}
  456. c                if (IND.eq.1.or.IND.eq.3) then 
  457. c              t_{n+1} ou t_{0} pour diff centreee ?
  458. c                if (IND.eq.1.or.IND.eq.2) then 
  459.                  IP2=IPALA(I,5)
  460. c c              t_{n+1/2} ou t_{n-1/2} (seulement pour devogelaere)
  461. c                else
  462. c                  IP2=IPALA(I,6)
  463. c                endif
  464.                CALL DYCPL1(IP1,IP2,XDEP,NPAS,PDT,XCONV)
  465.                XFL=XCPLGE*XCONV
  466.  
  467. cbp       - calcul d'un produit de convolution GRANGER_PAIDOUSSIS
  468. *TODO           (cas particulier HBM : la force fluide-elastique est 
  469. *            prise en compte dans Z directement)
  470.             elseif(jfonct.eq.101) then
  471.                IP1=IPALA(I,4)
  472.                IP2=IPALA(I,5)
  473.                IP3=IPALA(I,6)
  474.                VSD=XPALA(I,2)
  475.                XA0=XPALA(I,3)
  476.                CALL DYCPL2(IP1,IP2,IP3,VSD,XA0,XDEP,NPAS,PDT,XCONV)
  477.                XFL=XCPLGE*XCONV              
  478.  
  479. cbp       - calcul du produit par une fonction trigo cos(q)
  480. c           rem : on pourrait aussi prevoir cos(coef2*q) avec coef2=XPALA(I,2)
  481.  
  482.             elseif(jfonct.eq.1) then
  483.                 XFL = XCPLGE * COS(XDEP)
  484.                 IF(GETJAC) DFDX= -XCPLGE * SIN(XDEP)
  485.             elseif(jfonct.eq.2) then
  486.                 XFL = XCPLGE * SIN(XDEP)
  487.                 IF(GETJAC) DFDX= XCPLGE * COS(XDEP)
  488.  
  489. cbp       - calcul du produit par une fonction temporelle cos(wt)*q
  490.             elseif(jfonct.ge.11.and.jfonct.le.12) then
  491.                 XFREQ  = XPALA(I,2)
  492.                 XTIME = DBLE(NPAS-1)*PDT
  493.                 if(IND.EQ.2) XTIME=XTIME-PDTS2
  494.                 if(jfonct.eq.1) XFONCT = COS(XFREQ*XTIME)
  495.                 if(jfonct.eq.2) XFONCT = SIN(XFREQ*XTIME)
  496.                 XFL = XCPLGE * XDEP * XFONCT
  497.  
  498. c         - simple raideur (ou raideur en puissance)
  499.             else
  500.                   Xexpo = XPALA(I,2)
  501.                   if(Xexpo.eq.XONE) then
  502.                     XFL = XCPLGE * XDEP
  503.                     IF(GETJAC) DFDX = XCPLGE
  504.                   elseif(Xexpo.eq.3.D0) then
  505.                     XFL = XCPLGE * (XDEP**3)
  506.                     IF(GETJAC) DFDX = Xexpo*XCPLGE*(XDEP**2) 
  507.                   else
  508.                     XFL = XCPLGE * (XDEP**Xexpo)
  509.                     IF(GETJAC) DFDX = Xexpo*XCPLGE*(XDEP**(Xexpo-1.d0)) 
  510.                   endif
  511.             endif
  512. cbp, ici    prise en compte de la possibilite de liaison conditionnelle 
  513.             IF(iannul.ne.0) XFL=0.D0
  514.             XVALA(I,IND,1) = XFL
  515.             XVALA(I,IND,4) = XDEP
  516.             FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL
  517.             IF(GETJAC) THEN
  518.               KTOTXA(INA1,INA1)  = KTOTXA(INA1,INA1) + DFDX
  519.             ENDIF
  520.  
  521. *
  522. * ------ force elementaire de type POLYNOMIALE
  523. *
  524.          ELSE IF (ITYP.EQ.6) THEN
  525. *           nombre de modes "origine"
  526.             NMOD = IPALA(I,2)
  527.             CALL D2POL1(Q1,Q2,NA1,IPLIA,XPALA,XVALA,NLIAA,IND,PDT,
  528.      &                  NPAS,I,NMOD,FTOTA,IVINIT,iannul)
  529. *
  530. * ------ choc elementaire ...
  531. *
  532. *        ELSE IF (ITYP.EQ.  ) THEN
  533. *              .......
  534. *              .......
  535. *
  536.          ENDIF
  537.  
  538. *
  539. * la suite n'est plus utile car on passe iannul aux s_p de calcul des
  540. * forces de liaisons.
  541.  
  542. *       si la liaison etait annulee on l'annule
  543. *         if ( ( icond.eq. 1 ) .and. ( iannul.eq.1)) then
  544. *              on annule l'increment de ftotb
  545. *                do 112 ik = 1,na1
  546. *                  do 113 jk = 1,4
  547. *                        ftota (ik,jk) = ftota0(ik,jk)
  548. * 113              continue
  549. * 112            continue
  550.  
  551. *         end if
  552.  
  553.   10  CONTINUE
  554. *--------------------------------------------------------------------*
  555. *     FIN DE BOUCLE SUR LES LIAISONS
  556. *--------------------------------------------------------------------*
  557. *
  558.       END
  559.  
  560.  
  561.  
  562.  
  563.  
  564.  
  565.  
  566.  
  567.  
  568.  
  569.  
  570.  
  571.  
  572.  
  573.  

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