CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Dr. Wolfgang Send C DLR-Institute of Aeroelasticity, Goettingen, Germany C http://www.dlr.de/~wolfgang.send C C The programme FLT2DQ may be regarded as being an open source, C published on the occasion of the annual "German Aeronautical and C Astronautical Congress 2003", Munich 17-20 November 2003. C C >>> Disclaimer: C The author reserves the right not to be responsible for the C correctness, completeness or quality of the information provided. C Liability claims regarding damage caused by the use of any C information provided, including any kind of information which is C incomplete or incorrect, will therefore be rejected. C C The author appeals to the potential users for being honest and C leaving the source of their gained knowledge revealed. C C >>> Essential: C The programme requires the IMSL library being linked. C International Mathematical and Statistical Library C http://www.vni.com/index.html C C >>> Information: C The user programmable section defining the paths for the output C ist marked "OUTPUT PATHS". C The user programmable section directing the various output files C ist marked "DIRECTING FILES". C C >>> Included subroutines: C CNM3DK Calculates the thin plate 2D unsteady coefficients C for lift and moment in the calling programme C Incompressible flow using Theodorsen's function C Applied for Mach number <= 0.1 C C POSSIOF Calculates the thin plate 2D unsteady coefficients C in subsonic flow solving Possio's integral equation C Effective for Mach number > 0.1 C C FLT2DK_P Calculates the parameters for two types of output files C - Software Techplot (not Tecplot) C - Software Tecplot C C SCALE Scales the Tecplot axes according to respective data C C GETARG Reads the executable's command line parameters C The routine offers the option to include the respective C system call of the user's compiler C C >>> Input and output files: See below C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C *** FLT2DQ *** Created: 12-Mar-97 C Last revision: 26-May-03 C All comments in English C *** Updates *** C 08-Jul-04 The case fs=0 was not treated properly after the last C revision on 26-May-03 due to one misprint of 'rho' C C 24-Jan-05 Major improvement: Solution of Possio's integral C equation included for subsonic forces of the plunging C and pitching plate. Subroutine possiof provided by C Dr. Volker Carstens (DLR-Institute of Aeroelasticity) C C Reference: V. Carstens, Berechnung der instationaeren Druck- C verteilung an harmonisch schwingenden Gittern in ebener C Unterschallstroemung, Teil II, DFVLR IB 253-75J02 C DFVLR-AVA, Goettingen 1975 C C *----------------------------------------------------------------* C | M * X.tt + D * X.t + ( K + P ) * X = Q(X,X.t) + R(t) | C *----------------------------------------------------------------* C C Q(X,X.t) = -fs * ( QR * X + QS * X.t ) C C (qr11 qr12) (qs11/(2pi*f01) qs12/(2pi*f02) ) C QR = ( ) ; QS = ( ) C (qr21 qr22) (qs21/(2pi*f01) qs22/(2pi*f02) ) C C ( m1 s ) ( d1 0 ) ( k1 km) (q11 q12) C M = fm * ( ), D = ( ), K = ( ), Q = ( ) C ( s m2) ( 0 d2) ( km k2) (q21 q22) C C C ( +1 -1) (1) (X1) C P = kap * ( ), R = a0*cos(2pi*f0*t+ph0) * ( ), X = ( ) C ( -1 +1) (0) (X2) C C Important: Where do the QR and QS values come from? C ========== The input data CQij = (Rij,Sij) are assumed to be C complex fluid forces from an outside flow solution for C the input value omred=0, or they are inline provided C from the analytical 2D thin plate solution for omred>0. C Thus, the supporting force Q in the differential C equation is Q = -Q.fluid, because the Q.fluid is C compiled from the fluid's reaction to the system's C motion. Using the frequency domain values in a time C integration, the contribution of the time derivative C in the imaginary parts has to be eliminated: C QRij = Rij, QSi1 = S1j*uom01, QS2i = Si2*uom02 C where uom0i=1/(2pi*f0i), f0i the natural frequencies C from the eigenvalue equation det(-om**2*M+K) = 0. C If "strong coupling" is detected, the oms0i both are C taken from the coupling frequency for the second run. C C >>> Some analysis programmes, e.g. ZAERO, do not follow the C above explained sign convention, but take the fluid C force positive. In those cases, the coupling constant C fs (which is basically the stagnation pressure) has to C enter with a negative sign to meet the programme's C logic. C C >>> "Strong coupling" is defined by the distance of the C maxima in the power spectra of the two variables X1 and C X2. The power spectra are resolved according to the C frequency resolution of the fourier analysis. If the C frequencies of the maxima are more distant than KDIST C neighbouring frequency values, no strong coupling is C assumed. C C Notice: ph0 input is taken from "phase" [deg] converted to [rad] C C Understanding the Fourier analysis: C X1 may be identified as being the plunge mode C X2 may be identified as being the pitch mode C C M : Mass matrix - Derived from the 2 DOFs model C D : Damping matrix - No cross over influence assumed C K : Stiffness matrix - Elastic properties C C P : Coupling by a spring (alternative to km) C C Coupled system solved using IMSL-DIVPAG C 2. order system with 2 variables reduced to 1. order system with C four variables: C C *---------------------------------------------* C | X1.t = Y1, X2.t = Y2, X1 = Y3, X2 = Y4 | C *---------------------------------------------* C C Programme scheme: C C 1. Reading input data C 2. First loop C 2.1 Integrating the differential equation C 2.2 Fourier analysis of the solution C 3. Second loop if "strong coupling" is detected C 3.1 Integrating the differential equation C 3.2 Fourier analysis of the solution C 4. Output of the solution and the analysis data C 5. Parameters for graphical display using C - TechPlot (not Tecplot!) C - Tecplot C C Call: C {pgm} {sys} [{input}] C C {pgm}: Name of the executable C z.B.: flt2dq.exec, flt2dq.exe C {sys}: Respective operating system for which file is compiled C WIN/win - Windows 32 Bit environment C LNX/lnx - Linux C SGI/sgi - SGI Irix C {input}: Input file C [ ] optional, default is {input} = flt2dq.inp C C Output files: C C flt2dq.dat Result of the time integration and the fourier C analysis of the time history. C C flt2dq_fl.dat Flutter analysis in the frequency domain with C increasing stagnation pressure for constant C reduced frequency C (Output directed to the local directory) C C flt2dq.itp Techplot legend file (not Tecplot!) C flt2dq.lay Tecplot layout file C flt2dq.sty Tecplot style file C C Additional input files: C C flt2dq.lay0 Tecplot sample layout file C flt2dq.sty0 Tecplot sample style file C Both files are processed and adapted to the C respective intervals and output data returning C C Using Tecplot: The preset *.lay and *.sty files for Tecplot are C implemented in a macro file. The call is: C C >tecplot -p flt2dq.mcr C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Program FLT2DQ Implicit Logical(B),COMPLEX*16(C),Real*8(D-H,O-Y), Character(Z) Parameter (n=4,nc=n/2,nparam=50,istore=10001,maxstep=1000000) Parameter (input=35,ioutfl=36,idata=37) Dimension da(n,n),param(nparam),y(n),y0(n) Dimension cmata(nc,nc),cmatq(nc,nc),dmatm(nc,nc),dmatk(nc,nc), # cmatm_inv(nc,nc),dmatm_inv(nc,nc),cmatb(nc,nc), # ceval(nc),cevec(nc,nc) Dimension Tsto(istore),Y1sto(istore),Y2sto(istore), # Y3sto(istore),Y4sto(istore) Dimension FouF(istore),FouCY3(istore),FouPY3(istore), # FouCY4(istore),FouPY4(istore) Dimension Gen_f(2,istore),reg_p(2,istore),Sum_p(2,istore) Character z_pathw*40,z_pathl*40,z_paths*40,z_comment*40 Character z_zeile*30,z_sys*3,z_getarg*20,z_input*20 Real*8 m1,m2,k1,k2,kap,a0,km Real*8 ans,ans1,lmd1,lmd2,qh1,qh2 C ... For POSSIOF Dimension ip_ifvek(3),CQVAL(3,3) Character z_pout*20 Common /POSSIO/cqval,p_mach,p_redf,p_xr,p_xrud,ip_ifvek, * ip_idr,z_pout Common /FLT2DQPar/m1,m2,sm,d1,d2,k1,k2,km,kap,a0,f0,ph0,phase Common /FLT2DQEvl/fs,omred,xdr,dmach,uom01,uom02,f01,f02 Common /FLT2DQLay/t_ini,t_end,z_comment Common /QPar/RQ11,RQ12,RQ21,RQ22,SQ11,SQ12,SQ21,SQ22 Common /Gen_f/y1_ff,y2_ff,y1_rp,y2_rp External DIVPAG,DSET,UMACH,DEVCCG,DLINRG,ZGEMM External DFCN,DFCNJ Data pi2 /6.2831 85307 17958d0/,dlim/1d-10/ C C >>> OUTPUT PATHS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Notice: lpath{w/l/s} and lengths of the character strings C have to be equal. C xxx Data lpathw,z_pathw /21,'D:\TPWIN\USER\DGLR03\'/ C xxx Data lpathl,z_pathl /31,'/home/wsend/send/public/flt2dq/'/ C xxx Data lpaths,z_paths /31,'/usr/people/send/public/flt2dq/'/ Data lpathw,z_pathw /1,' '/ Data lpathl,z_pathl /1,' '/ Data lpaths,z_paths /1,' '/ C C CCC End of programmable section CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C --- Reading the parameters of the command line ----------------------- C In case that no equivalent routine may be found, the GETARG call C has to be commented out. The variable below has to be activated. call GETARG (1,z_getarg) If (z_getarg .eq. ' ') then z_getarg(1:3) = 'win' z_input = 'flt2dq.inp' BWIN = .true. Else z_sys = z_getarg(1:3) BWIN = z_sys .eq. 'WIN' .or. z_sys .eq. 'win' BSGI = z_sys .eq. 'SGI' .or. z_sys .eq. 'sgi' BLNX = z_sys .eq. 'LNX' .or. z_sys .eq. 'lnx' BUNIX = BSGI .OR. BLNX C Reading name for input file or data from flt2dq.inp z_getarg = ' ' call GETARG (2,z_getarg) C xxx z_getarg = 'DesignatedInputFile' If (z_getarg .ne. ' ') then z_input = z_getarg Else z_input = 'flt2dq.inp' End if End if c --- ------------------------------------------------------------------ C *** Constants, initial values, control data WRAD = PI2/3.6d2 URAD = 1D0/WRAD PI = PI2/2d0 C ... Set fixed options for subroutine POSSIOF (part 1) C See subroutine for legend of options (no aileron motion) ip_ifvek(1) = 1 ip_ifvek(2) = 1 ip_ifvek(3) = 0 ip_idr = 0 p_xrud = 0d0 z_pout = 'possio_sub.out' C ... Constant kdist defines the "strong coupling" criterion kdist = 5 Open (input,FILE=z_input) Read (input,'(F10.0,20X,A40)') t_ini,z_comment Read (input,'(F10.0)') t_end Read (input,'( I10)') npoints Read (input,'(F10.0)') y0(1) Read (input,'(F10.0)') y0(2) Read (input,'(F10.0)') y0(3) Read (input,'(F10.0)') y0(4) Read (input,'(F10.0)') m1 Read (input,'(F10.0)') m2 C ... Mass matrix may be scaled by fm Read (input,'( A30)') z_zeile If (z_zeile(25:25) .eq. '<') then Read (z_zeile,'(F10.0,4X,F10.0)') sm,fm If (fm .lt. 1d-10) fm = 1d0 else fm = 1d0 Read (z_zeile,'(F10.0,4X,F10.0)') sm End if m1 = m1*fm m2 = m2*fm sm = sm*fm Read (input,'(F10.0)') d1 Read (input,'(F10.0)') d2 C ... Stiffness matrix optional off diagonal term km Read (input,'( A30)') z_zeile If (z_zeile(25:25) .eq. '<') then Read (z_zeile,'(F10.0,4X,F10.0)') k1,km else km = 0d0 Read (z_zeile,'(F10.0,4X,F10.0)') k1 End if Read (input,'(F10.0)') k2 Read (input,'(F10.0)') kap Read (input,'(F10.0)') a0 Read (input,'(F10.0)') f0 Read (input,'(F10.0)') phase Read (input,'(F10.0)') tol Read (input,'( I10)') method Read (input,'( I10)') iout06 Read (input,'(A1)',end=5) z_end C ... Parameters for aerodynamic loads 5 If (z_end .eq. '.') then Read (input,'(F10.0,4(2X,F10.0))') fs,omred,xdr,dmach,scale Read (input,'(F10.0,2X,F10.0)') RQ11,SQ11 Read (input,'(F10.0,2X,F10.0)') RQ12,SQ12 Read (input,'(F10.0,2X,F10.0)') RQ21,SQ21 Read (input,'(F10.0,2X,F10.0)') RQ22,SQ22 Read (input,'(A1)',end=6) z_end b_aero = .true. else if (z_end .eq. 'i' .or. z_end .eq. 'I' ) then Read (input,'(F10.0,4(2X,F10.0))') fs,omred,xdr,dmach,scale Read (input,'(A1)',end=6) z_end b_aero = .true. else fs = 0d0 b_aero = .false. end if C ... Parameters for flutter calculation 6 If (z_end .eq. '.') then Read (input,'(F10.0)') x_ref Read (input,'(F10.0)') rho Read (input,'(F10.0)') q0_ini Read (input,'( I10)') iflint Read (input,'( I10)') iflplus Read (input,'( I10)') kdist b_flut = .true. else x_ref = 1d0 rho = 1.2d0 q0_ini = 0d0 iflint = 20 iflplus = 5 b_flut = .false. end if 7 Close(input) Write (06,'(/''*** Input data read -''/ # '' Computation has started. ''/)') C ... Phase in radians ph0 = phase * wrad C ... Too many points If (npoints .gt. istore) then Write (06,'(/''!!! Too many nodes chosen: '',I6/ # '' Maximum number allowed: '',I6/)') # npoints,istore stop '*** FLT2DQ terminated due to npoints > istore' end if C *** Aerodynamic loads if (abs(omred) .gt. dlim) then C ... 2D thin plate approximation If (abs(scale) .lt. dlim) scale = 1d0 C ... Effect of mach number included if (abs(dmach) .le. 1d-1 .or. dmach .lt. 0) then dmach = abs (dmach) fakma = 1d0/sqrt(1d0-dmach**2) C >>> Notice: Moments calculated clockwise (math. negative) CALL CNM3DK (omred,xdr,CNeta,CMeta,CNxsi,CMxsi) CIOMS = cmplx(0d0,omred) CNh = CIOMS*CNeta*pi2*fakma*scale CMh = CIOMS*CMeta*pi2*fakma*scale CNalf = (CNxsi+CNeta)*pi2*fakma*scale CMalf = (CMxsi+CMeta)*pi2*fakma*scale RQ11 = dreal(CNh) RQ12 = dreal(CNalf) RQ21 = -dreal(CMh) RQ22 = -dreal(CMalf) SQ11 = dimag(CNh) SQ12 = dimag(CNalf) SQ21 = -dimag(CMh) SQ22 = -dimag(CMalf) Write (06,'(/''>>> Aerodynamic loads from '', # ''2D thin plate theory. ''/)') else C Set variable options for subroutine POSSIOF (part 2) p_mach = dmach p_redf = omred p_xr = (xdr-5d-1)*2d0 call POSSIOF RQ11 = dreal(CQVAL(1,1))*scale RQ12 = dreal(CQVAL(1,2))*scale RQ21 = dreal(CQVAL(2,1))*scale RQ22 = dreal(CQVAL(2,2))*scale SQ11 = dimag(CQVAL(1,1))*scale SQ12 = dimag(CQVAL(1,2))*scale SQ21 = dimag(CQVAL(2,1))*scale SQ22 = dimag(CQVAL(2,2))*scale Write (06,'(/''>>> Aerodynamic loads from '', # ''2D Possio integral equation ''/)') end if If (fs .lt. 0d0) then Write (06,'(/''!!! Sign of input fs changed '', # ''from minus to plus. ''/)') fs = abs (fs) end if else Write (06,'(/''>>> Aerodynamic loads from '', # ''external coefficients. ''/)') end if BQIJ = .TRUE. C ... Input assigned to the mass matrix do 10 j=1,n do 10 i=1,n 10 da(i,j) = 0d0 do 20 i=1,n 20 da(i,i) = 1d0 da(1,1) = m1 da(2,2) = m2 da(2,1) = sm da(1,2) = sm C *** Eigenvalues and eigenvectors (compiled without damping) C Output from IMSL analysis C C ... Inverse of the mass matrix dmatm do 30 j=1,nc do 30 i=1,nc 30 dmatm(i,j) = da(i,j) call DLINRG(nc,dmatm,nc,dmatm_inv,nc) dmatk(1,1) = k1 dmatk(2,2) = k2 dmatk(2,1) = km dmatk(1,2) = km C ... Option for estimate of flutter boundary if (iflint .gt. 0) then bfldata = .true. open(ioutfl,file='flt2dq_fl.dat') write (ioutfl,'(/8x,''****** Flutter Limits *******''// # '' q0 [N/m**2]'','' u0 [m/s]'', # '' f1 [Hz]'','' d1 [1/s]'', # '' g1 [-]'', # '' f2 [Hz]'','' d2 [1/s]'', # '' g2 [-]'', # '' R-ev21 [-]'','' I-ev21 [-]'', # '' B-ev21 [-]'','' P-ev21 [-]'', # '' R-ev22 [-]'','' I-ev22 [-]'', # '' B-ev22 [-]'','' P-ev22 [-]'', # '' q0 [N/m**2]'')') else bfldata = .false. end if fs_dif = fs/dfloat(iflint) fs_ini = abs(q0_ini)*dsign(1d0,fs) fs_i = fs_ini do 50 ifl = 1,iflint+iflplus+1 fs_i = fs_i+fs_dif u0 = sqrt(2d0*abs(fs_i)/rho) cmatq(1,1) = fs_i*dcmplx(RQ11,SQ11) cmatq(2,2) = fs_i*dcmplx(RQ22,SQ22) cmatq(2,1) = fs_i*dcmplx(RQ21,SQ21) cmatq(1,2) = fs_i*dcmplx(RQ12,SQ12) do 40 j=1,nc do 40 i=1,nc cmata(i,j) = dcmplx(dmatk(i,j),0d0)+cmatq(i,j) 40 cmatm_inv(i,j) = dcmplx(dmatm_inv(i,j),0d0) call ZGEMM ('N','N',nc,nc,nc,(1d0,0d0),cmatm_inv,nc,cmata,nc, # (0d0,0d0),cmatb,nc) call DEVCCG (nc,cmatb,nc,ceval,cevec,nc) C ... Natural frequencies and eigenvectors cfrq1 = sqrt (ceval(1))/pi2 cfrq2 = sqrt (ceval(2))/pi2 if (abs (cevec(1,1)) .lt. dlim) cevec(1,1) = dcmplx(1d0,0d0) if (abs (cevec(1,2)) .lt. dlim) cevec(1,2) = dcmplx(1d0,0d0) cevec11 = cevec(1,1)/cevec(1,1) cevec21 = cevec(2,1)/cevec(1,1) cevec12 = cevec(1,2)/cevec(1,2) cevec22 = cevec(2,2)/cevec(1,2) revec21 = dreal(cevec21) sevec21 = dimag(cevec21) revec22 = dreal(cevec22) cuevl1 = (1d0,0d0)/ceval(1) cuevl2 = (1d0,0d0)/ceval(2) if (abs (u0) .lt. dlim) then omred1 = 0d0 omred2 = 0d0 else omred1 = dreal(ceval(1))*(x_ref/2d0)/u0 omred2 = dreal(ceval(2))*(x_ref/2d0)/u0 end if C gcoef1 = dimag(ceval(1))/dreal(ceval(1)) C gcoef2 = dimag(ceval(2))/dreal(ceval(2)) gcoef1 = dimag(cfrq1)/dreal(cfrq1)*pi2 gcoef2 = dimag(cfrq2)/dreal(cfrq2)*pi2 call vdmplr(revec21,sevec21,devec21,pevec21,-11) call vdmplr(revec22,sevec22,devec22,pevec22,-11) sevec22 = dimag(cevec22) If (bfldata) then write (ioutfl,'(F12.2,15F12.5,F12.2)') abs(fs_i), u0, # cfrq1,gcoef1,cfrq2,gcoef2, # cevec21,devec21,pevec21, # cevec22,devec22,pevec22,abs(fs_i) else write (06,'(/'' Eigenvalues and eigenvectors'', # '' (full system):''/)') write (06,'( '' frq1 = '',2F12.6)') cfrq1 write (06,'( '' vec11= '',2F12.6)') cevec11 write (06,'( '' vec21= '',2F12.6)') cevec21 write (06,'( '' frq2 = '',2F12.6)') cfrq2 write (06,'( '' vec12= '',2F12.6)') cevec12 write (06,'( '' vec22= '',2F12.6)') cevec22 end if 50 Continue if (bfldata) close (ioutfl) C *** The homogeneous system C Eigenvalues and eigenvectors (compiled without damping) C Output from the algebraic program REDUCE C ans1=(k1**2*m2**2-2.*k1*k2*m1*m2+4.*k1*k2*sm**2-4.* . k1*km*m2*sm+k2**2*m1**2-4.*k2*km*m1*sm+4.*km**2*m1* . m2)**0.5 ans=(0.5*ans1+0.5*k1*m2+0.5*k2*m1-km*sm)/(m1*m2-sm** . 2) LMD1 = ANS C ans1=(k1**2*m2**2-2.*k1*k2*m1*m2+4.*k1*k2*sm**2-4.* . k1*km*m2*sm+k2**2*m1**2-4.*k2*km*m1*sm+4.*km**2*m1* . m2)**0.5 ans=(-0.5*ans1+0.5*k1*m2+0.5*k2*m1-km*sm)/(m1*m2-sm . **2) LMD2 = ANS C ans1=(k1**2*m2**2-2.*k1*k2*m1*m2+4.*k1*k2*sm**2-4.* . k1*km*m2*sm+k2**2*m1**2-4.*k2*km*m1*sm+4.*km**2*m1* . m2)**0.5 det = (ans1* . m1-k1*m1*m2+2.0*k1*sm**2+k2*m1**2-2.0*km*m1*sm) If (abs (det) .gt. dlim) then ans=(-ans1*sm-k1*m2*sm-k2*m1*sm+2.0*km*m1*m2)/det QH1 = ANS else QH1 = 0d0 end if C ans1=(k1**2*m2**2-2.*k1*k2*m1*m2+4.*k1*k2*sm**2-4.* . k1*km*m2*sm+k2**2*m1**2-4.*k2*km*m1*sm+4.*km**2*m1* . m2)**0.5 det = (ans1* . m1-k1*m1*m2+2.0*k1*sm**2+k2*m1**2-2.0*km*m1*sm) If (abs (det) .gt. dlim) then ans=(-ans1*sm-k1*m2*sm-k2*m1*sm+2.0*km*m1*m2)/det QH2 = ANS else QH2 = 0d0 end if C C ... Allocating eigenvalues to natural frequencies f001 = sqrt (k1/m1)/pi2 f002 = sqrt (k2/m2)/pi2 ans1 = sqrt(LMD1)/pi2 ans2 = sqrt(LMD2)/pi2 f0a = ans1 f0b = ans2 if (f001 .gt. f002 .and. f0a .gt. f0b) then f01 = f0a f02 = f0b else f01 = f0b f02 = f0a end if write (06,'(/'' Eigenvalues and eigenvectors'', # '' (only structural part):''/)') write (06,'( '' frq1 = '', F12.6)') ans1 write (06,'( '' vec1 = '',2F12.6)') qh1,1d0 write (06,'( '' frq2 = '', F12.6)') ans2 write (06,'( '' vec2 = '',2F12.6)') qh2,1d0 write (06,'(/'' Time history of the solution:''/)') f00 = max(f01,f02) uom01 = 1d0/(pi2*f00) uom02 = 1d0/(pi2*f00) C C *** Numerical time integration C C ... Parameters C Second run for strong coupling BREPEAT = BQIJ 80 Continue Call DSET(nparam,0d0,param,1) ido = 1 C Method param(12) = method C Maximum number of steps allowed (def=500) param(4) = maxstep C Matrix used (constant) param(19) = 1 C ... Matrix type for da (if used) C 2: symmetric positive definite param(14) = 2 C Nonlinear solver method C 2: chord method (no user prov. Jacobian) param(13) = 2 C C *** Integration intval = npoints-1 integ = 0 t_dif = (t_end-t_ini)/float(max(1,intval)) t_int = t_ini t = t_ini C do 90 i=1,n 90 y(i) = y0(i) C *** Loop 100 Continue integ = integ+1 t_int = t_int+t_dif Call DIVPAG (ido,n,DFCN,DFCNJ,da,t,t_int,tol,param,y) Tsto (integ) = t_int Y1sto(integ) = y(1) Y2sto(integ) = y(2) Y3sto(integ) = y(3) Y4sto(integ) = y(4) gen_f(1,integ) = -y1_ff gen_f(2,integ) = -y2_ff reg_p(1,integ) = y1_rp reg_p(2,integ) = y2_rp if (integ .eq. 1) then sum_p(1,integ) = 0d0 sum_p(2,integ) = 0d0 else sum_p(1,integ) = sum_p(1,integ-1)+ # (-fs)*y1_ff*t_dif*y1sto(integ) sum_p(2,integ) = sum_p(2,integ-1)+ # (-fs)*y2_ff*t_dif*y2sto(integ) end if if (integ .lt. intval) then goto 100 else if (integ .eq. intval) then ido = 3 goto 100 end if C C *** Fourier analysis of the solution Call DFFTRF(npoints,Y3sto,FouCY3) Call DFFTRF(npoints,Y4sto,FouCY4) upkt = 1d0/float(npoints) u2pkt = 2d0*upkt limP = (npoints+1)/2 Periode = abs(t_end-t_ini) UT = 1d0/Periode Pfak = Periode/Float(npoints)**2 Pfak2 = Pfak+Pfak FouPY3(1) = Pfak*FouCY3(1)**2 FouPY4(1) = Pfak*FouCY4(1)**2 SumPY3 = FouPY3(1) SumPY4 = FouPY4(1) do 110 k = 2,limP FouF (k) = float(k-1)*UT PY3 = Pfak2*(FouCY3(2*k-2)**2+FouCY3(2*k-1)**2) PY4 = Pfak2*(FouCY4(2*k-2)**2+FouCY4(2*k-1)**2) FouPY3(k) = PY3 FouPY4(k) = PY4 SumPY3 = SumPY3+PY3 SumPY4 = SumPY4+PY4 110 Continue C Power spectrum normalized UsumPY3 = 1d0/SumPY3 UsumPY4 = 1d0/SumPY4 FouPY3(1) = UsumPY3*FouPY3(1) FouPY4(1) = UsumPY4*FouPY4(1) FOUPY3max = 0d0 FOUPY4max = 0d0 do 120 k = 2,limP FouPY3(k) = UsumPY3*FouPY3(k) FouPY4(k) = UsumPY4*FouPY4(k) FOUPY3k = FOUPY3(k) FOUPY4k = FOUPY4(k) If (FOUPY3k .gt. FOUPY3max) then FOUPY3max = FOUPY3k kPY3 = k End if If (FOUPY4k .gt. FOUPY4max) then FOUPY4max = FOUPY4k kPY4 = k End if 120 Continue If (Abs(kPY3-kPY4) .gt. kdist) then write (06,'('' No strong coupling detected. '')') Else write (06,'(''*** Strong coupling detected: '')') Dif_freq = Abs(FouF(kPY3)-FouF(kPY4)) write (06,'('' Frequency margin '',F5.3, '' Hz'')') # Dif_freq kmax = kPY3 C Evaluating phase shift and amplitude ratio FouFmax = FouF(kmax) d1_cos = FOUCY3(2*kmax-2)*u2pkt d1_sin = FOUCY3(2*kmax-1)*u2pkt call vdmplr(d1_cos,d1_sin,d1_abs,d1_phas,11) d2_cos = FOUCY4(2*kmax-2)*u2pkt d2_sin = FOUCY4(2*kmax-1)*u2pkt call vdmplr(d2_cos,d2_sin,d2_abs,d2_phas,11) C ... ... Advancing of Y1 (plunge) ahead of Y2 (pitch) C Sign convention of pitch is already taken into account C (If Y1 is expected to start "at the bottom" starten, the C initial value has to be set as Y1 -> -Y1. dif_phas = d1_phas-d2_phas if (dif_phas .lt. 0d0) dif_phas = dif_phas+360d0 dkap = dif_phas dlmb = d1_abs/d2_abs Write (06,'(/'' f [Hz]'', # '' lambda'','' kappa'', # '' Y1_0'','' Y2_0'')') write (06,'(F08.3,F08.4,F08.2,2F08.2/)') # FouFmax,dlmb,dkap,d1_abs,d2_abs C ... ... Second run with new om01, om02 in the case of strong coupling IF (BREPEAT) THEN BREPEAT = .false. uom01 = 1d0/(pi2*FouFmax) uom02 = 1d0/(pi2*FouFmax) Write (06,'(/''Second run with the '', # ''coupling frequency. ''/)') Goto 80 End if End if C Analysis in the lower half of the frequency range do 130 k = limp+1,npoints FouF (k) = float(k-1)*UT FouPY3(k) = 0d0 FouPY4(k) = 0d0 130 Continue u0 = sqrt(2d0*abs(fs)/rho) if (abs(u0) .lt. dlim) then omred_new = 0d0 else omred_new = pi2*FouFmax*(x_ref/2d0)/u0 end if Write (06,'(/'' f [Hz] '',''rho[kg/m3]'','' x_ref[m]'', # '' u0[m/s]'','' Ma[-]'','' omred'', # '' q0 '',''[N/m2]'')') write (06,'(F08.3,2F10.2,F08.2,2F06.3,F08.0/)') # FouFmax,rho,x_ref,u0,dmach,omred_new,abs(fs) C If (abs (fs) .gt. 1d-10) then C ... Printing the Q Matrix Write (06,'(/'' *** QR,QS matrices as eventually applied'')') Write (06,'( '' '','' QR11 '','' QS11 '', # '' QR12 '','' QS12 '')') Write (06,'(5X,4(1PG12.3))') # RQ11, SQ11/uom01, RQ12, SQ12/uom02 Write (06,'( '' '','' QR21 '','' QS21 '', # '' QR22 '','' QS22 '')') Write (06,'(5X,4(1PG12.3)/)') # RQ21, SQ21/uom01, RQ22, SQ22/uom02 end if C C *** Saving and printing results C 150 Continue C C >>> DIRECTING FILES CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C if (BWIN)then Open (idata,FILE=z_pathw(1:lpathw)//'flt2dq.dat') OPEN (16,FILE='flt2dq.lay0') OPEN (17,FILE='flt2dq.sty0') OPEN (25,FILE=z_pathw(1:lpathw)//'flt2dq.itp') OPEN (26,FILE=z_pathw(1:lpathw)//'flt2dq.lay') OPEN (27,FILE=z_pathw(1:lpathw)//'flt2dq.sty') else if (BLNX) then Open (idata,FILE='flt2dq.dat') OPEN (16,FILE=z_pathl(1:lpathl)//'flt2dq.lay0') OPEN (17,FILE=z_pathl(1:lpathl)//'flt2dq.sty0') OPEN (26,FILE='flt2dq.lay') OPEN (27,FILE='flt2dq.sty') else if (BSGI) then Open (idata,FILE='flt2dq.dat') OPEN (16,FILE=z_paths(1:lpaths)//'flt2dq.lay0') OPEN (17,FILE=z_paths(1:lpaths)//'flt2dq.sty0') OPEN (26,FILE='flt2dq.lay') OPEN (27,FILE='flt2dq.sty') end if C C CCC End of programmable setion CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C If (BLNX) then Else if (BSGI) then Else if (BWIN) then Else Write (06,'(/''!!! FLT2DQ - '', # ''No valid OS found for output file''/)') stop '!!! FLT2DQ stopped due do wrong option {sys}' end if Write (idata, # '(''Title = "FLT2DQ Data - Programmed by W. Send" '')') Write (idata, # '('' Variables = "t", "Y1", "Y2", "Y3", "Y4", '', # ''"FouF", "FouCY3", "FouPY3", "FouCY4", "FouPY4", '', # ''"GenF1", "GenF2", "RPart1", "RPart2", "SumP1", '', # ''"SumP2", "SumPAll" '')') Write (idata,'('' Zone I = '',I5,'', F=Point '')') npoints Write (idata,'(''# '',I09,4F12.5)') npoints,f01,0.01,f02,0.01 Write (idata,'(''# '',I09,4F12.5)') 0 ,f01,0.10,f02,0.10 Write (idata,'(17(1PG14.5))') t_ini,y0,FouF(1), # FouCY3(1),FouPY3(1), # FouCY4(1),FouPY4(1), # 0d0, 0d0, 0d0, 0d0, # 0d0, 0d0, 0d0 if (iout06 .gt. 0) then write (iout06,'(/8x,''****** Results *******''// # '' t '','' y(1) '', # '' y(2) '','' y(3) '', # '' y(4) '')') write (iout06,'(12F12.5)') t_ini,y0,FouF(1), # FouCY3(1),FouPY3(1), # FouCY4(1),FouPY4(1) write (iout06,'(1H )') end if C *** Output and extreme values y_max = 0d0 ydot_max = 0d0 do 200 i=1,intval if (i .eq. intval) write (06,'(1H )') write (idata,'(17(1PG14.5))') # Tsto(i) ,Y1sto(i),Y2sto(i), # Y3sto(i),Y4sto(i),FouF(i+1), # FouCY3(i+1),FouPY3(i+1), # FouCY4(i+1),FouPY4(i+1), # gen_f(1,i), gen_f(2,i), # reg_p(1,i), reg_p(2,i), # sum_p(1,i), sum_p(2,i), # sum_p(1,i)+ sum_p(2,i) y_max = max( y_max, abs(y3sto(i)),abs(y4sto(i)) ) ydot_max = max( ydot_max, abs(y1sto(i)),abs(y2sto(i)) ) if (iout06 .gt. 0) # write (iout06,'(12F12.5)') Tsto(i) ,Y1sto(i),Y2sto(i), # Y3sto(i),Y4sto(i),FouF(i+1), # FouCY3(i+1),FouPY3(i+1), # FouCY4(i+1),FouPY4(i+1) 200 continue close (idata) Call FLT2DQ_P (bunix,y_max,ydot_max,n,y0) close (25) close (16) close (26) close (17) close (27) C End C C *** Subroutines Subroutine DFCN(n,t,y,ys) Implicit Real*8(D-H,O-Z) Dimension y(n),ys(n) Real*8 m1,m2,k1,k2,kap,a0,km Common /FLT2DQPar/m1,m2,sm,d1,d2,k1,k2,km,kap,a0,f0,ph0,phase Common /FLT2DQEvl/fs,omred,xdr,dmach,uom01,uom02,f01,f02 Common /QPar/RQ11,RQ12,RQ21,RQ22,SQ11,SQ12,SQ21,SQ22 Common /Gen_f/y1_ff,y2_ff,y1_rp,y2_rp Data pi2/6.2831 85307 17958D0/ DATA ICNT /0/,T0 /0D0/ C C ... Counting calls ICNT = ICNT+1 DIFT = T-T0 T0 = T C ... external force extf = a0*cos(pi2*f0*t+ph0) C C ... Fluid forces taken negative y1_ff = SQ11*uom01*y(1) + SQ12*uom02*y(2) * + RQ11*y(3) + RQ12*y(4) y2_ff = SQ21*uom01*y(1) + SQ22*uom02*y(2) * + RQ21*y(3) + RQ22*y(4) C C ... Regular part of diff. eq. y1_rp = -d1*y(1)-(k1+kap)*y(3)-(km-kap)*y(4) y2_rp = -d2*y(2)-(km-kap)*y(3)-(k2+kap)*y(4) C C ... Right side of diff. eq. ys(1) = y1_rp - fs*y1_ff +extf ys(2) = y2_rp - fs*y2_ff ys(3) = y(1) ys(4) = y(2) Return End C Subroutine DFCNJ(n,t,y,dypdy) Implicit Real*8(D-H,O-Z) Dimension y(n),dypdy(n) Real*8 m1,m2,k1,k2,kap,a0,km Common /FLT2DQPar/m1,m2,sm,d1,d2,k1,k2,km,kap,a0,f0,ph0,phase Common /FLT2DQEvl/fs,omred,xdr,dmach,uom01,uom02,f01,f02 Common /QPAR/RQ11,RQ12,RQ21,RQ22,SQ11,SQ12,SQ21,SQ22 C ... Partial derivatives (in this case never called) Return End