C DYCPL2 SOURCE BP208322 18/01/10 21:15:32 9684 C SUBROUTINE DYCPL2(IP1,IP3,IP4,VSURD,XA0,XDEP,NPAS,PDT,XCONV) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) *--------------------------------------------------------------------* * * * Operateur DYNE : * * ________________ * * * * Calcul de la force fluidelastique par convolution. * * (cas particulier du modele de Granger Paidoussis) * * * * Parametres : * * * * e IP1 LISTREEL de grandeurs precalculees du modele * * e//s IP3 LISTREEL y_i(t_n-1) // y_i(t_n) * * e//s IP4 LISTREEL S_i(t_n-1) // S_i(t_n) * * e VSURD FLOTTTANT V/D * * e XA0 FLOTTTANT alpha_0 = 1 - \sum alpha_i * * e XDEP FLOTTTANT q(t_n) * * e PDT FLOTTANT = \Delta t * * s XCONV FLOTTANT int_0^t_n h(\tau)*Qj(t-\tau) d\tau * * * * si devogelaere (et pas differences_centrees), on a plutot : * * e PDT FLOTTANT = \Delta t/2 * * e//s IP3 LISTREEL y_i(t_n-1/2) // y_i(t_n) * * e//s IP4 LISTREEL S_i(t_n-1/2) // S_i(t_n) * * * * Auteur : * * BP, 2017-12-19 * * * *--------------------------------------------------------------------* -INC SMLREEL POINTEUR MLREE4.MLREEL c LOGICAL LWRITE * RECUPERATION * des listreels deja actif en entree (ip3 et ip4 sont modifiables) MLREE1=IP1 c MLREE2=IP2 c --> optimisation par precalcul dans 1 seul listreel MLREE3=IP3 MLREE4=IP4 c JG1=MLREE1.PROG(/1) c JG=JG1/3 JG=MLREE4.PROG(/1) c CALCUL DE L'INTEGRALE DE CONVOLUTION XCONV = 0.D0 c boucle sur i (degré d'approximation du modele (souvent 2)) DO 1 I=1,JG cc alpha_i et delta_i c Ai = MLREE1.PROG(I) c Di = MLREE2.PROG(I) c AiDi = A1 * Di c AUXi = Di*VSURD*PDT c --> optimisation par precalcul c alpha_i * delta_i AiDi = MLREE1.PROG(I) c YSHIFT=EXP(+delta_i*V/D*PDT) si DIFF_CENT (ou PDT/2 si DEVOGE) YSHIFT=MLREE1.PROG(I+JG) c XSHIFT=EXP(-delta_i*V/D*PDT) si DIFF_CENT (ou PDT/2 si DEVOGE) XSHIFT=MLREE1.PROG(I+2*JG) * recup des valeurs du pas precedent Yim1 = MLREE3.PROG(I) Sim1 = MLREE4.PROG(I) * decalage temporel * (pour eviter les exp(grand*t)=infini et exp(-grand*t)=0) c XSHIFT = EXP(-1.D0*AUXi) cc profitons du fait que exp(0)=1 cc Yim1 = Yim1 * XSHIFT c --> optimisation par precalcul Sim1 = Sim1 * XSHIFT * calcul des y_i(t=T) c Yi = EXP(AUXi)*XDEP c --> optimisation par precalcul Yi = YSHIFT * XDEP * contribution a l'integrale du pas : dSi = PDT * 0.5D0 * (Yi + Yim1) Si = Sim1 + dSi c dernier produit et cumul XCONVi = AiDi * XSHIFT * Si c IF(LWRITE) write(*,*) '>dycpl2 i=',i,AUXi,Yim1,Sim1,Yi,Si XCONV = XCONV + XCONVi * stockage pour porchain pas de temps c MLREE3.PROG(I) = Yi MLREE3.PROG(I) = XDEP MLREE4.PROG(I) = Si 1 CONTINUE c AJOUT DU TERME DE PURE RAIDEUR XSTIFF = XA0 * XDEP / VSURD XCONV = XCONV + XSTIFF RETURN END