PROGRAM PROP C C 22-JAN-93 CDP 67 PLUS/VVPR N675- 1325/ T. Rose provided SSS C Alter updates C C PROGRAM TO CREATE DMIG CARDS TO ADD THE PROPELLER AERODYNAMIC TERMS C TO THE AERO SOLUTION C SOL 75 - USE DMIG WITH CASE CONTROL CALL OUT FOR EACH SET C C SETS ARE BASED ON VELOCITY,MACH NUMBER AND DENSITY C C ONE SUBCASE MUST EXIST FOR EACH SET - USE DMAP ALTER TO INCLUDE C THEM C C REAL MU,ETA,K,I1,I2,I3,J1,J2,J3,IP,M DIMENSION STIFF(4,4),DAMP(4,4) DIMENSION KGRID(2,1000),KDOF(2,1000),TERM(2,1000),NVAL(2) DIMENSION WASH(4,2000),IWDOF(2,2000) CHARACTER TITLE*80,ANS*1,INFILE*32,OUTFILE*32,DMIGFIL*32 INTEGER DNWASH COMMON OMEGARS,V,R,A,RHO,RHO0,R0(40),C0(40),NR0,I1,I2,I3,J1,J2,J3, + CONST1,CONST2,IP,MU,CR,M,BLADE,IPRINT,B1,B2,A0RAT, + TITLE,LAG,A0,AM,OMEGA,IGRID,IZDOF,IYDOF,ICASE,NVEL, + GEARAT,PROPIP,TURIP,DNWASH,DMIGFIL,IDOF1,IDOF2, + KGRID,KDOF,TERM,NVAL C C VERSION ID FOR PRINTOUT C VERSN = 3.0 C C VERSION 3.0 - ADDED DOWNWASH TERMS - 11/1988 C C ADDED TURBINE INERTIA AND GEAR RATIO C C CONSISTENT UNITS ASSUMED - POUND AND INCH UNLESS OTHERWISE NOTED C C INPUT DATA: C C OMEGA = PROPELLER ROTATION (RPM - INTERNALLY CHANGED TO RAD/SEC) C V = VELOCITY (input in ft/sec and converted internally) C SSOUND= SPEED OF SOUND(input in ft/sec and used to find m) C TURIP = WEIGHT MOMENT OF INERTIA OF TURBINE (used to calc IP) C GEARAT= GEAR RATIO " " C PROPIP= WEIGHT MOMENT OF INERTIA OF PROP " " C NOTE : GRAV = 386,0886 USED TO CONVERT WEIGHT TO MASS C IP = PROPELLER POLAR MASS MOMENT OF INERTIA - NO LONGER INPUT C CALCULATED AS IP=(PROPIP+GEARAT*TURIP)/GRAV C A0 = LIFT CURVE SLOPE - SCALING FACTOR FOR INTEGRALS - USES A0/ C AM = MAXIMUM LIFT CURVE SLOPE C TITLE = PROBLEM TITLE TO USE ON OUTPUT C A = ASPECT RATIO C LAG = LAG FLAG - 0 = CALCULATE LAG COEFFICIENTS C 1 = SET LAG COEFFICIENTS TO 0.0 (NO LAG) C RHO = DENSITY C RHO0 = REFERENCE DENSITY C NR0 = NUMBER OF SECTIONS USED TO DEFINE PROP C R0(I) = RADIUS(I) C C0(I) = CHORD(I) C BLADE = NUMBER OF BLADES ON THE PROP C CR = REFERENCE CHORD C M = MACH NUMBER C IGRID = GRID NUMBER FOR DMIG MATRICES C IYP = PITCH AXIS AT IGRID (NORMALLY 2, BUT IF COORD DON'T LINE U C THIS CAN BE USED TO ALIGN THE MATRICES C IZP = YAW AXIS AT GRID (NORMALLY 3, BUT IF COORD DON'T LINE UP, C THIS CAN BE USED TO ALIGN THE MATRICES C IPRINT= PRINT FLAG - 0 = ONLY INPUT C 1 = 0 + CALCULATED TERMS C 2 = 0 + 1 + MATRICES C 3 = 2 + INTERMEDIATE VALUES DURING CALCS C DNWASH - FLAG TO INCLUDE DOWNWASH TERMS FROM DMIG FILE V C 0 = NO DOWNWASH TERMS C >0 = INCLUDE DOWNWASH TERMS (NEXT ITEMS ARE REQUIRED) C DMIGFIL = FILE WITH DMIG CARDS FOR DOWNWASH TERMS (ONLY IF DNWASH> C IDOF1 = PROP DOF ASSOCIATED WITH FIRST TERM IN PARTN1 IN MSC/ C NASTRAN PART 1 RUN (0=NOT USED) C IDOF2 = PROP DOF ASSOCIATED WITH SECOND TERM IN PARTN1 IN MSC/ C NASTRAN PART 1 RUN (0=NOT USED) C C..................................................................... IPAGE = 1 ICASE = 1 ICONT = 100000 GRAV = 386.0886 c c set nvel = 0 as a flag for input of extra velocities from screen c nvel=0 C C CHECK ON TYPE OF INPUT DESIRED C WRITE(*,10) 10 FORMAT(' PROGRAM TO CREATE DMIG CARDS TO ADD THE ', ' 'PROPELLER AERODYNAMIC TERMS',/,' TO THE AERO SOLUTION', ' ' IN MSC/NASTRAN',//,' IS YOUR INPUT IN A FILE?(Y/N)?') READ(*,'(A1)')ANS IF(ANS.EQ.'Y')GOTO 20 IF(ANS.EQ.'y')GOTO 20 C C ENTER DATA FROM TERMINAL C CALL INPUT1 OPEN(UNIT=21,STATUS='NEW',FILE='AERODMIG.DAT') GO TO 25 20 CONTINUE C C DATA IS IN A INPUT FILE C write(*,21) 21 format(' Please enter the name of your input file') read(*,'(a32)')infile open(unit=5,status='old',file=infile) CALL INPUT2 25 CONTINUE IP = ( PROPIP + GEARAT * TURIP ) / GRAV IZP = IZDOF IYP = IYDOF IZPR = IZP + 3 IYPR = IYP + 3 write(*,26) 26 format (' Do you want your output to go to a file?') READ(*,'(A1)')ANS if(ans.ne.'Y'.and.ans.ne.'y')go to 28 write(*,27) 27 format(' Enter file name for output') read(*,'(a32)')outfile open(unit=6,status='new',file=outfile) 28 continue if(nvel.ne.0)go to 29 OPEN(UNIT=21,STATUS='NEW',FILE='AERODMIG.DAT') 29 TWOPI = 6.283185307 A0RAT = A0 / TWOPI OMEGARS = OMEGA * TWOPI / 60. C IF ASPECT RATIO INPUT, THEN DON'T CALCULATE IF(A.NE.0.)GO TO 56 C C CALCULATE ASPECT RATIO C ETAOLD = R0(1) / R0(NR0) CROLD = C0(1) / CR TOTAL = 0. IF(IPRINT.GE.3)WRITE(6,49) 49 FORMAT(' CALCULATING ASPECT RATIO OF PROP',/, + ' STEP C RATIO ETA DELTAETA CHANGE ', + ' TOTAL SO FAR',/) DO 55 IT = 2, NR0 CRNEW = C0(IT) / CR ETANEW = R0(IT) / R0(NR0) DELTA = ETANEW - ETAOLD CHANGE = (CRNEW + CROLD) * DELTA / 2. TOTAL = TOTAL + CHANGE CROLD = CRNEW ETAOLD = ETANEW IF(IPRINT.GE.3)WRITE(6,54)IT,CRNEW,ETANEW,DELTA,CHANGE,TOTAL 54 FORMAT(I5,5F10.6) 55 CONTINUE c correction for error 29655 - remove factor of 2 - changed to radius... c c A = 2.*R0(NR0) / (CR * TOTAL) c A = R0(NR0) / (CR * TOTAL) A = A * (1. - R0(1)/R0(NR0))**2 C C ASPECT RATIO CALCULATED C 56 CONTINUE C C CALCULATE MU AND CONSTANTS C 50 MU = V/(OMEGARS*R0(NR0)) CONST1 = MU*A/CR CONST2 = MU/CONST1 C C ALWAYS ECHO INPUT C WRITE(6,2000)VERSN,IPAGE,TITLE 2000 FORMAT(1H1,' PROGRAM TO GENERATE TERMS FOR USE IN PROPELLER', + ' WHIRL DYNAMICS IN MSC/NASTRAN',/,T32,'VERSION ',F4.2,T70, + 'PAGE',I5,//,1X,A80,//,' DATA AS INPUT:',/) IPAGE = IPAGE + 1 ILINE = 6 WRITE(6,2100)BLADE,A,TURIP,GEARAT,PROPIP,IP,CR,A0,AM,OMEGA, + OMEGARS,NR0,(I,R0(I),C0(I),I=1,NR0) 2100 FORMAT(' PROPELLER GEOMETRY DATA:',/, + ' NUMBER OF BLADES = ',F3.0, + ' ASPECT RATIO =',F9.3,/, + ' TURBINE WEIGHT MOMENT OF INERTIA = ',1PE13.6, + ' GEAR RATIO = ',0PF6.2,/, + ' PROPELLER WEIGHT MOMENT OF INERTIA = ',1PE13.6,/, + ' POLAR MASS MOMENT OF INERTIA = ',1PE13.6,/, + ' REFERENCE CHORD LENGTH = ',0PF6.3,/, + ' A0 = ',0PF10.7,' AM =',0PF10.6,/, + ' OMEGA = ',1PE13.6,' RPM = ',1PE13.6,' RADIANS/SEC',/, + ' PROPELLER GEOMETRY DEFINED BY',0PI3,' SECTIONS',/, + ' SECTION RADIUS CHORD',/,40(0PI5,F11.2,F10.4,/)) ILINE = ILINE + 6 + NR0 IF(ILINE.GT.45)CALL NEWPAGE(VERSN,ILINE,IPAGE,TITLE) if(lag.eq.0)write(6,2196) 2196 format(' Lag terms are included',/) if(lag.eq.1)write(6,2197) 2197 format(' Lag terms are not included',/) iline = iline + 2 WRITE(6,2200)ICASE,RHO,RHO0,V/12.,M 2200 FORMAT(' DATA FOR CASE',I4,/, + ' DENSITY = ',1PE13.6, +' REFERENCE DENSITY = ',1PE13.6,/, + ' VELOCITY = ',1PE13.6, +' ft/sec MACH NUMBER = ',1PE13.6,/) ILINE = ILINE + 3 C C END OF INPUT ECHO C C C SET LAG TERMS - for version not using Bessel Functions C B1 = 0. B2 = 0. IF(LAG.EQ.0)B1 = .165 IF (LAG.EQ.0)B2=.335 C INTEGRATE EQUATIONS C CALL INTGRT C C NOW HAVE I1-3, J1-3 C TEMP = OMEGARS*CR/V CZTHETA = -4.*TEMP*I1 CMTHETA = -2.*TEMP*J2 CYTHETA = -4.*TEMP*J1 CNTHETA = -2.*TEMP*I2 CZQ = -2.*CMTHETA CMQ = -2.*TEMP*I3 CYQ = 2.*CNTHETA CNQ = -2.*TEMP*J3 CZPSI = CYTHETA CMPSI = -CNTHETA CYPSI = -CZTHETA CNPSI = CMTHETA CZR = CYQ CMR = -CNQ CYR = -CZQ CNR = CMQ C ********************* C C Czq,Cnq,Cmr,Cyr set = 0.0 since sign convention not well defined C 1/18/1988 CZQ=0. CNQ=0. CMR=0. CYR=0. C ********************* C C ALL PROP AERO TERMS FOUND - CALC GYRO TERMS C CINERT = IP*OMEGARS C C PRINT INPUT AND ALSO RESULTANTS IF DESIRED C C C ADDITIONAL PRINTOUT IF DESIRED C IF(IPRINT.EQ.0)GO TO 4000 C C PRINT CALCULATED TERMS C IF(ILINE.GT.47)CALL NEWPAGE(VERSN,ILINE,IPAGE,TITLE) WRITE(6,2300)CZTHETA,CMTHETA,CYTHETA,CNTHETA 2300 FORMAT(/,' CALCULATED TERMS FOR AERO MATRICES',/, + ' CZTHETA = ',1PE13.6,' CMTHETA = ',1PE13.6,/, + ' CYTHETA = ',1PE13.6,' CNTHETA = ',1PE13.6) ILINE = ILINE + 3 IF(ILINE.GT.47)CALL NEWPAGE(VERSN,ILINE,IPAGE,TITLE) WRITE(6,2400)CZQ ,CMQ ,CYQ ,CNQ 2400 FORMAT( + ' CZQ = ',1PE13.6,' CMQ = ',1PE13.6,/, + ' CYQ = ',1PE13.6,' CNQ = ',1PE13.6) ILINE = ILINE + 2 IF(ILINE.GT.47)CALL NEWPAGE(VERSN,ILINE,IPAGE,TITLE) WRITE(6,2500)CZPSI ,CMPSI ,CYPSI ,CNPSI 2500 FORMAT( + ' CZPSI = ',1PE13.6,' CMPSI = ',1PE13.6,/, + ' CYPSI = ',1PE13.6,' CNPSI = ',1PE13.6) ILINE = ILINE + 2 IF(ILINE.GT.47)CALL NEWPAGE(VERSN,ILINE,IPAGE,TITLE) WRITE(6,2600)CZR ,CMR ,CYR ,CNR 2600 FORMAT( + ' CZR = ',1PE13.6,' CMR = ',1PE13.6,/, + ' CYR = ',1PE13.6,' CNR = ',1PE13.6,/) ILINE = ILINE + 2 IF(ILINE.GT.47)CALL NEWPAGE(VERSN,ILINE,IPAGE,TITLE) C C PRINT GYRO TERM C WRITE(6,2700)CINERT 2700 FORMAT(' GYRO DAMPING TERM = ',1PE13.6) ILINE = ILINE + 1 C C DONE PRINTING CALCULATED TERMS C 4000 CONTINUE C C CALC MATRIX TERMS C C C STIFFNESS TERMS C C CONST = RHO * V*V * R0(NR0)**3 * 3.14159 CONST1 = CONST/(2.*R0(NR0)) c c store values of constants for prop downwash c con=const con1=const1 c **************** C write(6,99501)con,con1 C99501 format(' Const =',e14.8,' const1 = 'e14.8) c DO 4010 I = 1,4 DO 4005 J = 1,4 DAMP(I,J) = 0.0 4005 STIFF(I,J) = 0.0 4010 CONTINUE STIFF(1,2) = -CZTHETA * CONST1 STIFF(2,2) = -CMTHETA * CONST STIFF(3,2) = -CYTHETA * CONST1 STIFF(4,2) = -CNTHETA * CONST STIFF(1,4) = -CZPSI * CONST1 STIFF(2,4) = -CMPSI * CONST STIFF(3,4) = -CYPSI * CONST1 STIFF(4,4) = -CNPSI * CONST C C DAMPING MATRIX C CONST1 = .5 * RHO * V * R0(NR0)**2 * 3.14159 C.... C CONST2 = CONST1 TIMES D/2 CONST2 = CONST1 * R0(NR0) DAMP(1,1) = -CZTHETA * CONST1 DAMP(1,2) = -CZQ * CONST2 C .... SIGN CHANGED 10/31/88 AS PER WPR...TLR C DAMP(1,3) = -CZPSI * CONST1 DAMP(1,3) = CZPSI * CONST1 DAMP(1,4) = -CZR * CONST2 DAMP(2,1) = -CMTHETA * CONST1 * 2. * R0(NR0) DAMP(2,2) = -CMQ * CONST2 * 2. * R0(NR0) C .... SIGN CHANGED 10/31/88 AS PER WPR...TLR C DAMP(2,3) = -CMPSI * CONST1 * 2. * R0(NR0) DAMP(2,3) = CMPSI * CONST1 * 2. * R0(NR0) DAMP(2,4) = -CMR * CONST2 * 2. * R0(NR0) + CINERT DAMP(3,1) = -CYTHETA * CONST1 DAMP(3,2) = -CYQ * CONST2 C .... SIGN CHANGED 10/31/88 AS PER WPR...TLR C DAMP(3,3) = -CYPSI * CONST1 DAMP(3,3) = CYPSI * CONST1 DAMP(3,4) = -CYR * CONST2 DAMP(4,1) = -CNTHETA * CONST1 * 2. * R0(NR0) DAMP(4,2) = -CNQ * CONST2 * 2. * R0(NR0) - CINERT C .... SIGN CHANGED 10/31/88 AS PER WPR...TLR C DAMP(4,3) = -CNPSI * CONST1 * 2. * R0(NR0) DAMP(4,3) = CNPSI * CONST1 * 2. * R0(NR0) DAMP(4,4) = -CNR * CONST2 * 2. * R0(NR0) C C PRINT MATRICES IF REQUESTED C IF(IPRINT.LT.2) GO TO 5000 IF(ILINE.GT.42) CALL NEWPAGE(VERSN,ILINE,IPAGE,TITLE) IZPR = IZP + 3 IYPR = IYP + 3 c WRITE(6,4200)IGRID,IZP,IZPR,IYP,IYPR WRITE(6,4200)IGRID,IZP,IYPR,IYP,IZPR 4200 FORMAT(' AERODYNAMIC STIFFNESS MATRIX WITHOUT DOWNWASH ', + 'FOR GRID ',I8,//, + ' DOF',3(I6,10X),I6,/) ILINE = ILINE + 2 WRITE(6,4300)IZP,(STIFF(1,J),J=1,4) 4300 FORMAT(I4,4(1PE13.6,3X)) WRITE(6,4300)IYPR,(STIFF(2,J),J=1,4) C WRITE(6,4300)IZPR,(STIFF(2,J),J=1,4) WRITE(6,4300)IYP,(STIFF(3,J),J=1,4) C WRITE(6,4300)IYPR,(STIFF(4,J),J=1,4) WRITE(6,4300)IZPR,(STIFF(4,J),J=1,4) ILINE = ILINE + 4 IF(ILINE.GT.42) CALL NEWPAGE(VERSN,ILINE,IPAGE,TITLE) C C WRITE DAMPING MATRIX C WRITE(6,4400)IGRID,IZP,IYPR,IYP,IZPR 4400 FORMAT(' AERODYNAMIC DAMPING MATRIX FOR GRID ',I8,//, + ' DOF',3(I6,10X),I6,/) ILINE = ILINE + 2 WRITE(6,4300)IZP,(DAMP(1,J),J=1,4) C WRITE(6,4300)IZPR,(DAMP(2,J),J=1,4) WRITE(6,4300)IYPR,(DAMP(2,J),J=1,4) WRITE(6,4300)IYP,(DAMP(3,J),J=1,4) C WRITE(6,4300)IYPR,(DAMP(4,J),J=1,4) WRITE(6,4300)IZPR,(DAMP(4,J),J=1,4) ILINE = ILINE + 4 C C MATRICES WRITTEN TO OUTPUT C 5000 CONTINUE C C CHECK IF DOWNWASH TERMS EXIST C IF(DNWASH.EQ.0)GO TO 4150 IF(IPRINT.GE.2) CALL NEWPAGE(VERSN,ILINE,IPAGE,TITLE) C C CALC DOWNWASH TERMS C CALL DOWNWASH(WASH,CON,CON1,KGRID,KDOF,TERM, + IZP,IYP,IGRID,IZPR,IYPR,IPRINT,IDOF1,IDOF2,NVAL,IWDOF, + NVFOUND,CZTHETA,CMTHETA,CYTHETA,CNTHETA,CZPSI,CMPSI, + CYPSI,CNPSI) 4150 CONTINUE C C CREATE DMIG CARDS IN FILE 'AERODMIG.DAT' - FILE 21 C WRITE(21,5030) 5030 FORMAT('$ RESULTS FROM PROGRAM PROP FOR THE FOLLOWING PROBLEM') WRITE(21,5050)TITLE 5050 FORMAT('$ ',A80) WRITE(21,5100)ICASE+100,ICASE+100,V/12.,M,RHO 5100 FORMAT('$ THE FOLLOWING MATRICES (AEROS',I3,' AND AEROB',I3, + ') ARE BASED ON ',/,'$ VELOCITY = ',1PE13.6,'(ft/sec)'/, +'$ MACH NO = ',1PE13.6,/, +'$ DENSITY = ',1PE13.6) IF(DNWASH.EQ.0) + CALL DMIG(STIFF,DAMP,IGRID,IZP,IZPR,IYP,IYPR,ICASE,ICONT) IF(DNWASH.GT.0)CALL DMIG1(STIFF,DAMP,IGRID,IZP,IZPR, + IYP,IYPR,ICASE,ICONT,WASH,IWDOF,IDOF1,IDOF2,NVAL) C C CHECK FOR MORE CASES C if(nvel.eq.0)go to 6000 if(icase.ge.nvel)go to 9999 call input1 go to 6050 6000 continue READ(5,*,END=9999,ERR=9999)V,SSOUND,RHO M = V / SSOUND v=v*12. SSOUND = SSOUND*12. 6050 ICASE = ICASE + 1 GO TO 50 9999 WRITE(6,10000) if(ans.ne.'Y'.and.ans.ne.'y')go to 10001 write(*,10000) 10000 FORMAT(' DONE') 10001 continue END C C SUBROUTINES .................................................. C SUBROUTINE INTGRT C C SUBROUTINE TO INTEGRATE I1-3,J1-3 USING TRAPEZOIDAL RULE C REAL MU,ETA,K,I1,I2,I3,J1,J2,J3,K041,K320,I1NEW,I1OLD,IP,M REAL I2NEW,I2OLD,I3NEW,I3OLD,J1NEW REAL J1OLD,J2NEW,J2OLD,J3NEW,J3OLD CHARACTER TITLE*80 CHARACTER DMIGFIL*32 INTEGER DNWASH DIMENSION KGRID(2,1000),KDOF(2,1000),TERM(2,1000),NVAL(2) COMMON OMEGARS,V,R,A,RHO,RHO0,R0(40),C0(40),NR0,I1,I2,I3,J1,J2,J3, + CONST1,CONST2,IP,MU,CR,M,BLADE,IPRINT,B1,B2,A0RAT, + TITLE,LAG,A0,AM,OMEGA,IGRID,IZDOF,IYDOF,ICASE,NVEL, + GEARAT,PROPIP,TURIP,DNWASH,DMIGFIL,IDOF1,IDOF2, + KGRID,KDOF,TERM,NVAL C CHECK = 1.-(A0/AM)**2 C C FIRST POINT - START INTEGRATION C I1 = 0. I2 = 0. I3 = 0. J1 = 0. J2 = 0. J3 = 0. I = 1 ETA = R0(I)/R0(NR0) SMMEE = SQRT( MU*MU + ETA*ETA ) DENOM = M * M * ( 1 + ETA*ETA / ( MU*MU )) IF(DENOM.GT.CHECK)DENOM=CHECK DENOM = 1 - DENOM CC **************************** C WRITE(6,1 9985)ETA,MU,M,DENOM C19985 FORMAT(' ETA = ',F10.6,' MU =',F10.6,' M = ',F10.6,' DENOM =', C + F13.8) CC **************************** DENOM = 2. + A * SQRT(DENOM) DENOM = SMMEE * DENOM K = C0(I) / (2. * R0(NR0) * SMMEE ) K041 = 1. + (.041 / K)**2 K320 = 1. + (.320 / K)**2 if(lag.ne.0)go to 8 c c use Theodorsen function to find F and G if lag included c CALL STPBS0(K,1,BJ0,BY0) CALL STPBS1(K,1,BJ1,BY1) DENOMb=(BJ1+BY0)**2 + (BY1-BJ0)**2 fk= (BJ1*(BJ1+BY0) + BY1*(BY1-BJ0) )/DENOMb gk=-(BY1*BY0 + BJ1*BJ0)/DENOMb C (CR + I*CI = THEODORSEN FUNCTION) go to 9 8 continue fk=1. gk=0. 9 continue c cc old lag term calculation - before Bessel function added c FK = 1. - B1/K041 - B2/K320 c GK = -B1 * (.041 / K) / K041 - B2 * (.320 / K) / K320 J1NEW = C0(I) * GK / DENOM I1NEW = C0(I) * FK / DENOM I2NEW = I1NEW * ETA * ETA J2NEW = J1NEW * ETA * ETA I3NEW = I2NEW * ETA * ETA J3NEW = J2NEW * ETA * ETA IF(IPRINT.GT.2)WRITE(6,8900)I IF(IPRINT.GT.2)WRITE(6,10)SMMEE,FK,GK IF(IPRINT.GT.2)WRITE(6,9000)DENOM,K,ETA,MU,M IF(IPRINT.GT.2)WRITE(6,9100)FK,GK,I1NEW, + J1NEW,I2NEW,J2NEW,I3NEW,J3NEW C C LOOP ON RADIUS C DO 1000 I = 2,NR0 I1OLD = I1NEW I2OLD = I2NEW I3OLD = I3NEW J1OLD = J1NEW J2OLD = J2NEW J3OLD = J3NEW ETAOLD = ETA ETA = R0(I)/R0(NR0) SMMEE = SQRT( MU*MU + ETA*ETA ) DENOM = M * M * ( 1 + ETA*ETA / ( MU*MU )) IF (DENOM.GT.CHECK)DENOM = CHECK DENOM = 1 - DENOM CC **************************** C WRITE(6,19986)ETA,MU,M,DENOM C19986 FORMAT(' SECOND SET',/, C + ' ETA = ',F10.6,' MU =',F10.6,' M = ',F10.6,' DENOM =', C + F13.8) CC **************************** DENOM = 2. + A * SQRT(DENOM) DENOM = SMMEE * DENOM K = C0(I) / (2. * R0(NR0) * SMMEE ) 10 FORMAT(' SMMEE =',1PE13.6,' F(K) = ',1PE13.6, + ' G(K) =',1PE13.6) IF(IPRINT.GT.2)WRITE(6,8900)I 8900 FORMAT(' INTEGRATION STEP ',I2,/) IF(IPRINT.GT.2)WRITE(6,9000)DENOM,K,ETA,MU,M, + I1OLD,J1OLD,I2OLD,J2OLD,I3OLD,J3OLD 9000 FORMAT(' VALUES IN SUBROUTINE INTGRT',/, + ' DENOM =',1PE13.6,' K =',1PE13.6,' ETA = ',1PE13.6,/, + ' MU = ',1PE13.6,' M = ',1PE13.6,/, + ' I1OLD = ',1PE13.6,' J1OLD = ',1PE13.6,/, + ' I2OLD = ',1PE13.6,' J2OLD = ',1PE13.6,/, + ' I3OLD = ',1PE13.6,' J3OLD = ',1PE13.6,/) K041 = 1. + (.041 / K)**2 K320 = 1. + (.320 / K)**2 if(lag.ne.0)go to 9014 c c use Theodorsen function to find F and G c CALL STPBS0(K,1,BJ0,BY0) CALL STPBS1(K,1,BJ1,BY1) DENOMb=(BJ1+BY0)**2 + (BY1-BJ0)**2 fk= (BJ1*(BJ1+BY0) + BY1*(BY1-BJ0) )/DENOMb gk=-(BY1*BY0 + BJ1*BJ0)/DENOMb C (CR + I*CI = THEODORSEN FUNCTION) go to 9016 9014 fk=1. gk=0. 9016 continue c cc old lag terms - before Bessel function added c FK = 1. - B1/K041 - B2/K320 c GK = -B1 * (.041 / K) / K041 - B2 * (.320 / K) / K320 IF(IPRINT.GT.2)WRITE(6,10)SMMEE,FK,GK J1NEW = C0(I) * GK / DENOM I1NEW = C0(I) * FK / DENOM I2NEW = I1NEW * ETA * ETA J2NEW = J1NEW * ETA * ETA I3NEW = I2NEW * ETA * ETA J3NEW = J2NEW * ETA * ETA IF(IPRINT.GT.2)WRITE(6,9100)FK,GK,I1NEW, + J1NEW,I2NEW,J2NEW,I3NEW,J3NEW 9100 FORMAT(' FK =',1PE13.6,' GK = ',1PE13.6,/, + ' I1NEW = ',1PE13.6,' J1NEW = ',1PE13.6,/, + ' I2NEW = ',1PE13.6,' J2NEW = ',1PE13.6,/, + ' I3NEW = ',1PE13.6,' J3NEW = ',1PE13.6,/) C C TRAPEZOIDAL INTEGRATION C DELTAE = ETA - ETAOLD I1 = I1 + (I1OLD + I1NEW) * DELTAE / 2. I2 = I2 + (I2OLD + I2NEW) * DELTAE / 2. I3 = I3 + (I3OLD + I3NEW) * DELTAE / 2. J1 = J1 + (J1OLD + J1NEW) * DELTAE / 2. J2 = J2 + (J2OLD + J2NEW) * DELTAE / 2. J3 = J3 + (J3OLD + J3NEW) * DELTAE / 2. IF(IPRINT.GT.2)WRITE(6,9200)DELTAE,I1,J1,I2,J2,I3,J3 9200 FORMAT(' DELTAE =',1PE13.6,/, + ' I1 = ',1PE13.6,' J1 = ',1PE13.6,/, + ' I2 = ',1PE13.6,' J2 = ',1PE13.6,/, + ' I3 = ',1PE13.6,' J3 = ',1PE13.6,/) C C END OF LOOP ON RADIUS C 1000 CONTINUE I1 = A0RAT * MU*MU * A * I1 / CR I2 = A0RAT * I2 * MU * A / CR I3 = A0RAT * I3 * A / CR J1 = A0RAT * J1 * MU*MU * A / CR J2 = A0RAT * J2 * MU * A / CR J3 = A0RAT * J3 * A / CR IF(IPRINT.GT.2)WRITE(6,9300)I1,J1,I2,J2,I3,J3 9300 FORMAT(' LOOP COMPLETED - VALUES BEFORE BLADE CORRECTION',/, + ' I1 = ',1PE13.6,' J1 = ',1PE13.6,/, + ' I2 = ',1PE13.6,' J2 = ',1PE13.6,/, + ' I3 = ',1PE13.6,' J3 = ',1PE13.6,/) IF(BLADE.EQ.4.) GO TO 2000 BRAT = BLADE / 4. I1 = I1 * BRAT I2 = I2 * BRAT I3 = I3 * BRAT J1 = J1 * BRAT J2 = J2 * BRAT J3 = J3 * BRAT IF(IPRINT.GT.2)WRITE(6,9400)I1,J1,I2,J2,I3,J3 9400 FORMAT(' VALUES AFTER BLADE CORRECTION',/, + ' I1 = ',1PE13.6,' J1 = ',1PE13.6,/, + ' I2 = ',1PE13.6,' J2 = ',1PE13.6,/, + ' I3 = ',1PE13.6,' J3 = ',1PE13.6,/) 2000 CONTINUE RETURN END C C END OF SUBROUTINE INTGRT C SUBROUTINE NEWPAGE(VERSN,ILINE,IPAGE,TITLE) CHARACTER TITLE*80 WRITE(6,100)VERSN,IPAGE,TITLE 100 FORMAT(1H1,' PROGRAM TO GENERATE TERMS FOR USE IN PROPELLER', + ' WHIRL DYNAMICS IN MSC/NASTRAN',/,T32,'VERSION ',F4.2,T70, + 'PAGE',I5,//,1X,A80,//,' DATA AS INPUT:',/) ILINE = 3 IPAGE = IPAGE + 1 RETURN END C C END OF SUBROUTINE NEWPAGE C SUBROUTINE DMIG(STIFF,DAMP,IGRID,IZP,IZPR,IYP,IYPR,ICASE,ICONT) c DIMENSION STIFF(4,4),DAMP(4,4),IDOF(4),IC(3) DIMENSION STIFF(4,4),DAMP(4,4),IDOF(4),IC(5) C C SUBROUTINE TO CREATE DMIG CARDS IN FILE 21 C C C STIFFNESS TERMS FIRST C ICAS=ICASE+100 WRITE(21,100)ICAS,0,1,1 100 FORMAT('DMIG AEROS',I3,3I8,/) IDOF(1) = IZP IDOF(2) = IYPR IDOF(3) = IYP IDOF(4) = IZPR IC(1) = ICONT - 4 IC(2) = ICONT - 3 IC(3) = ICONT - 2 IC(4) = ICONT - 1 IC(5) = ICONT C C Write stiffness matrix to file AERODMIG.DAT DO 1000 I = 1,4 DO 590 I99 = 1,5 590 IC(I99) = IC(I99) + 5 1000 WRITE(21,200)ICAS,IGRID,IDOF(I),IC(1), + (IC(J),IGRID,IDOF(J),STIFF(J,I),IC(J+1),J=1,4) 200 FORMAT('DMIG *AEROS',I3,2I16,T73,'*S',I6,/, + 4('*S',I6,2I16,1PE16.9,16X,'*S',I6/)) C C WRITE DAMPING matrix to file AERODMIG.DAT C WRITE(21,1100)ICAS,0,1,1 1100 FORMAT('DMIG AEROB',I3,3I8,/) IDOF(1) = IZP IDOF(2) = IYPR IDOF(3) = IYP IDOF(4) = IZPR C DO 2000 I = 1,4 DO 1590 I99 = 1,5 1590 IC(I99) = IC(I99) + 5 2000 WRITE(21,2200)ICAS,IGRID,IDOF(I),IC(1), + (IC(J),IGRID,IDOF(J),DAMP(J,I),IC(J+1),J=1,4) 2200 FORMAT('DMIG *AEROB',I3,2I16,T73,'*D',I6,/, + 4('*D',I6,2I16,1PE16.9,16X,'*D',I6,/)) ICONT = IC(5) + 1 RETURN END C SUBROUTINE INPUT1 C C SUBROUTINE TO READ INPUT FROM THE TERMINAL C COMMON OMEGARS,V,R,A,RHO,RHO0,R0(40),C0(40),NR0,I1,I2,I3,J1,J2,J3, + CONST1,CONST2,IP,MU,CR,M,BLADE,IPRINT,B1,B2,A0RAT, + TITLE,LAG,A0,AM,OMEGA,IGRID,IZDOF,IYDOF,ICASE,NVEL, + GEARAT,PROPIP,TURIP,DNWASH,DMIGFIL,IDOF1,IDOF2, + KGRID,KDOF,TERM,NVAL DIMENSION KGRID(2,1000),KDOF(2,1000),TERM(2,1000),NVAL(2) REAL IP,M CHARACTER DMIGFIL*32 INTEGER DNWASH CHARACTER TITLE*80,ANS*1,INFILE*32 c c if nvel is not 0, then only need to input changes for new case c if(nvel.ne.0)go to 2000 WRITE(*,100) 100 FORMAT(' INPUT FROM TERMINAL - PLEASE ENTER VALUES AS REQUESTED') WRITE(*,200) 200 FORMAT('(units are inch,pound,radian, and second unless noted)'/) WRITE(*,300) 300 FORMAT(' Enter problem title (80 characters or less)') read(*,'(A80)')title write(*,400) 400 format(/,' Enter propeller rate of revolution(RPM)') read(*,*)omega write(*,420) 420 format(' Enter the polar weight moment of inertia ', + 'of the turbine, the gear ratio, and the ',/,' ', + ' weight moment of inertia of the prop') read(*,*)turip,gearat,propip write(*,500) 500 format(/,' Enter lift curve slope for integrals ',/, + '(enter 0 to use default value of 2*pi)') read(*,*)a0 if(a0.eq.0.)a0=6.283185307 write(*,550) 550 format(' Enter maximum lift curve slope for integrals ',/, + '(enter 0 to use default value of 4*pi)') read(*,*)am if(am.eq.0.)am=6.283185307*2. write(*,600) 600 format(' Enter flag for consideration of lag effects', + ' (0=calculate lag terms, 1 = no lag terms)') read(*,*)lag if(lag.eq.0)write(*,700) 700 format(' Lag terms are to be calculated') if(lag.ne.0)write(*,800) 800 format(' Lag terms are not to be included') WRITE(*,900) 900 FORMAT(' Enter reference air density') read(*,*)rho0 write(*,1000) 1000 format(/,' Propeller data:',//, + ' How many blades does the prop have?') read(*,*) blade write(*,1100) 1100 format(' How many sections will you use to define the', + ' prop geometry?') read(*,*)nr0 1200 format(' Enter propeller geometry data starting at the hub') do 1400 i=1,nr0 write(*,1300)i 1300 format(' Enter Radius and Chord for section',I3) read(*,*)r0(i),c0(i) 1400 continue write(*,1450) 1450 format(' Enter reference chord and aspect ratio') write(*,1455) 1455 format(' (if input aspect ratio=0.0, then it will be ', + 'calculated)') read(*,*)cr,a write(*,1475) 1475 format(' Enter the NASTRAN GRIDid the prop matrices will', + ' attach to') read(*,*)igrid write(*,1480)igrid 1480 format(' Enter dof to be used for YAW terms for GRID ',I8, + /,' (Check manual to be sure of positive orientation!)') read(*,*)iydof write(*,1490)igrid 1490 format(' Enter dof to be used for PITCH terms for GRID ',I8, + /,' (Check manual to be sure of positive orientation!)') read(*,*)izdof C C VERSION 3.0 C write(*,1492) 1492 format(' Do you have a downwash matrix(1=yes,0=no)?') read(*,*)dnwash if(dnwash.eq.0)go to 1499 c c downwash input c write(*,1493) 1493 format(' Enter the name of the file with the DMIG downwash', + ' data in it') read(*,*)dmigfil write(*,1494) 1494 format(' Enter the dof of the prop grid associated with the ', + ' first column in the DMIG matrix (1-3)') read(*,*)idof1 write(*,1495) CALL READDMIG(DMIGFIL,KDOF,KGRID,TERM,1,NVAL,IPRINT) 1495 format(' Enter the dof of the prop grid associated with the ', + ' second column in the DMIG matrix (1-3)') read(*,*)idof2 IF(IDOF2.NE.0) +CALL READDMIG(DMIGFIL,KDOF,KGRID,TERM,2,NVAL,IPRINT) 1499 continue c c end of version 3.0 additions c write(*,1500) 1500 format(' How many different cases do you wish to consider?') read(*,*)nvel 2000 continue c c input values which change from case to case c write(*,2100)icase+1 2100 format(' Enter Velocity(ft/sec), Speed of sound(ft/sec),', + ' and air Density for case',i2) read(*,*)v,ssound,rho m = v/ssound v = v*12. SSOUND=SSOUND*12. c c temporily set iprint - can be set to 1,2, or 3 iprint = 2 c RETURN END SUBROUTINE INPUT2 COMMON OMEGARS,V,R,A,RHO,RHO0,R0(40),C0(40),NR0,I1,I2,I3,J1,J2,J3, + CONST1,CONST2,IP,MU,CR,M,BLADE,IPRINT,B1,B2,A0RAT, + TITLE,LAG,A0,AM,OMEGA,IGRID,IZDOF,IYDOF,ICASE,NVEL, + GEARAT,PROPIP,TURIP,DNWASH,DMIGFIL,IDOF1,IDOF2, + KGRID,KDOF,TERM,NVAL DIMENSION KGRID(2,1000),KDOF(2,1000),TERM(2,1000),NVAL(2) CHARACTER DMIGFIL*32 INTEGER DNWASH REAL IP,M CHARACTER TITLE*80,ANS*1,INFILE*32 c c subroutine to read input from a file c read(5,'(a80)')title read(5,*)omega,turip,gearat,propip,a0,am,lag,rho0,blade,cr,a,nr0 IF(A0.EQ.0.)A0=6.2831853 IF(AM.EQ.0.)AM=6.2831853*2. read(5,*)(r0(i),c0(i),i=1,nr0) read(5,*)igrid,iydof,izdof iprint = 3 C read(5,*)v,ssound,rho,IPRINT c c additions for version 3.0 c read(5,*)dnwash if(dnwash.eq.0)go to 500 read(5,'(a32)')dmigfil read(5,*)idof1 IPRINT = 2 CALL READDMIG(DMIGFIL,KDOF,KGRID,TERM,1,NVAL,IPRINT) read(5,*)idof2 if(idof2.ne.0) + CALL READDMIG(DMIGFIL,KDOF,KGRID,TERM,2,NVAL,IPRINT) c c end of version 3.0 additions c 500 IPRINT = 2 read(5,*)v,ssound,rho m = v/ssound V = V * 12. SSOUND = SSOUND*12. RETURN END SUBROUTINE STPBS0(X,NCODE,BJ0,BY0) C SUBROUTINE BES0. J AND Y BESSEL FUNCTIONS OF ORDER ZERO C E. ALBANO, ORGN 3721, EXT 1022, OCT. 1967 C COMPUTES J0(X) IF X IS GREATER THAN -3. C COMPUTES Y0(X) IF (X IS GREATER THAN E AND NCODE = 1 ), C WHERE INTEGER NAME(2) DATA NAME /4HSTPB,4HS0 / E=0.00001 C REF. US DEPT OF COMMERCE HANDBOOK (AMS 55) PG. 369 A=ABS(X) IF(A-3.) 10,10,100 10 Z=X*X/9. BJ0=1.+Z*(-2.2499997+Z*(1.2656208+Z*(-0.3163866+Z*(0.0444479 1 +Z*(-0.0039444+Z* 0.00021))))) IF(NCODE-1) 15,20,15 15 RETURN 20 IF(X-E) 200,25,25 25 BY0=0.63661977*BJ0*(ALOG(X)-.69314718)+.36746691+Z*(0.60559366+Z* 1 (-0.74350384+Z*(0.25300117+Z*(-0.04261214+Z*(0.00427916 2 -0.00024846*Z))))) RETURN 100 IF(X ) 250,250,110 110 U=1./SQRT(X) Z=3./X W=0.79788456+Z*(-0.00000077+Z*(-0.0055274+Z*(-0.00009512+Z* 1 (0.00137237+Z*(-0.00072805+0.00014476*Z))))) T=X-0.78539816+Z*(-0.04166397+Z*(-0.00003954+Z*(0.00262573+Z* 1 (-0.00054125+Z*(-0.00029333+0.00013558*Z))))) UW=U*W BJ0=UW*COS(T) IF(NCODE-1) 15,120,15 120 BY0=UW*SIN(T) 1000 RETURN 200 CONTINUE 250 CONTINUE GO TO 1000 END SUBROUTINE STPBS1(X,NCODE,BJ1,BY1) C SUBROUTINE BES1 J AND Y BESSEL FUNCTIONS OF FIRST ORDER C E. ALBANO, ORGN 3721, EXT 1022, OCT 1967 C COMPUTES J1(X) IF X IS GREATER THAN -3. C COMPUTES Y1(X) IF (X IS GREATER THAN E AND NCODE = 1 ), C WHERE INTEGER NAME(2) DATA NAME /4HSTPB,4HS1 / E=0.00001 C REF. US DEPT OF COMMERCE HANBOOK (AMS 58) PG. 370 A=ABS(X) IF(A-3.) 10,10,100 10 Z=X*X/9. BJ1=X*(0.5+Z*(-0.56249985+Z*(0.21093573+Z*(-0.03954289+Z* 1 (0.00443319+Z*(-0.00031761+0.00001109*Z)))))) IF(NCODE-1) 15,20,15 15 RETURN 20 IF(X-E) 200,25,25 25 BY1=0.63661977*BJ1*(ALOG(X)-.69314718)+(-0.6366198 +Z* 1 (0.2212091+Z*(2.1682709 +Z*(-1.3164827+Z*(0.3123951+Z* 2 (-0.0400976+0.0027873*Z))))))/X RETURN 100 IF(X) 250,250,110 110 U=1./SQRT(X) Z=3./X W=0.79788456+Z*(0.00000156+Z*(0.01659667+Z*(0.00017105+Z* 1 (-0.00249511+Z*(0.00113653-0.00020033*Z))))) T=X-2.35619449+Z*(0.12499612+Z*(0.00005650+Z*(-0.00637879+Z* 1 (0.00074348+Z*(0.00079824-0.00029166*Z))))) UW=U*W BJ1=UW*COS(T) IF(NCODE-1) 15,120,15 120 BY1=UW*SIN(T) 1000 RETURN 200 CONTINUE 250 CONTINUE GO TO 1000 END C C SUBROUTINE READDMIG(DMIGFIL,KDOF,KGRID,TERM,ICOL,NVAL,IPRINT) C C SUBROUTINE TO READ TERMS FROM MSC/NASTRAN DMIG FILE C C DMIGFIL = FILE WITH MSC/NASTRAN DMIG MATRIX C KGRID,KDOF = ASSOCIATED GRID AND DOF FOR TERMS C TERM = MATRIX TERMS C ICOL = COLUMN NUMBER TO USE IN MATRIX C DIMENSION KGRID(2,1000),KDOF(2,1000),TERM(2,1000),NVAL(2) CHARACTER DMIGFIL*32,CARD*16,MATRIX*8 C C OPEN FILE C OPEN(UNIT=11,STATUS='OLD',FILE=DMIGFIL) IERROR = 0 C C READ MATRIX NAME AND TYPE C READ(11,100)CARD,I1,I2,I3,I4 100 FORMAT(A16,4I8) MATRIX=CARD(9:16) WRITE(*,200)MATRIX 200 FORMAT(/,' MATRIX ',A8,' IS BEING READ FOR DOWNWASH TERMS',/) IF(I2.EQ.9.OR.I2.EQ.2)GO TO 400 C C WRONG MATRIX FORM C write(*,300)I2 300 FORMAT(' W A R N I N G : MATRIX IS NOT RECTANGULAR,', + ' FORM IS ',I2, ' PROGRAM IS CONTINUING ANYWAY') IERROR = 1 400 CONTINUE C C SET MATRIX TYPE: 1 = SINGLE PRECISION, 2 = DOUBLE PRECISION C ITYPE = 1 IF(I3.EQ.2.OR.I3.EQ.4)ITYPE = 2 C C READ INPUT DATA C I = 1 IC = 0 1000 IF(ITYPE.EQ.1)READ(11,1100,END=9999)CARD,IGRID,IDOF,VALUE IF(ITYPE.EQ.2)READ(11,1200,END=9999)CARD,IGRID,IDOF,VALUE 1100 FORMAT(A16,I8,I16,E16.9) 1200 FORMAT(A16,I8,I16,D16.9) IF(CARD(9:16).NE.MATRIX)GO TO 1300 C C NEW COLUMN C IC = IC + 1 IF(IC.GT.ICOL)GO TO 9999 if(iprint.ge.3)write(6,1250)icol c if(iprint.ge.3)write(*,1250)icol 1250 format(' Reading column ',i3) GO TO 1000 1300 CONTINUE C C CHECK IF CURRENT COLUMN IS THE RIGHT ONE C IF(ICOL.NE.IC)GO TO 1000 KGRID(IC,I) = IGRID KDOF(IC,I) = IDOF C CHANGE SIGN OF PROP WASH TERMS - 11/22/88 - AS PER WPR TERM(IC,I) = -VALUE NVAL(ICOL) = I if(iprint.ge.3)write(6,5000)igrid,idof,value c if(iprint.ge.3)write(*,5000)igrid,idof,value 5000 format(' GRID ',i8,' DOF ',i8,' has downwash term =',e13.6) I = I + 1 GO TO 1000 9000 WRITE(6,9100)ICOL,MATRIX 9100 FORMAT(' COLUMN ',I2,' OF MATRIX ',A8,' HAS BEEN READ') 9999 RETURN END C SUBROUTINE DOWNWASH(WASH,CON,CON1,KGRID,KDOF,TERM,IZP,IYP, + IGRID,IZPR,IYPR,IPRINT,IDOF1,IDOF2,NVAL,IWDOF,NVFOUND, + CZTHETA,CMTHETA,CYTHETA,CNTHETA,CZPSI,CMPSI,CYPSI,CNPSI) C C CALCULATE DOWNWASH TERMS AND PUT IN WASH C DIMENSION KGRID(2,1000),KDOF(2,1000),WASH(4,2000) DIMENSION TERM(2,1000),NVAL(2),IWDOF(2,2000) C WRITE(*,1) C 1 FORMAT(' IN SUBROUTINE DOWNWASH') C C WASH IS 4 PROP DOF X 1000 OTHER DOF C C FIRST DOF IN WASH IS IZP C SECOND DOF IN WASH IS IYPR C THIRD DOF IN WASH IS IYP C FOURTH DOF IN WASH IS IZPR C C IWDOF IS THE ASSOCIATED WASH GRIDS AND DOF C NOTE THERE IS THE POSSIBILITY OF REPEAT GRIDS AND DOF IN C THE PITCH AND YAW TERMS, SO THEY HAVE TO BE CHECKED C C ZERO WASH C c*********** ipold=iprint C iprint=3 c*********** NV=NVAL(1)+NVAL(2) DO 10 I = 1,4 DO 10 J = 1,NV 10 WASH(I,J)=0.0 C C CHECK IF FIRST DOF IS PITCH C C WRITE(*,14)IDOF1 C 14 FORMAT(' IDOF1 = ',I2) IF(IDOF1.EQ.IYP)GO TO 1000 C C FIRST DOF IS PITCH C IF(IPRINT.GE.3)WRITE(*,15) 15 FORMAT(' First DOF of downwash terms is pitch') DO 100 I = 1,NVAL(1) WASH(1,I)=WASH(1,I)-CON1*CZTHETA*TERM(1,I) WASH(2,I)=WASH(2,I)-CON*CMTHETA*TERM(1,I) WASH(3,I)=WASH(3,I)-CON1*CYTHETA*TERM(1,I) WASH(4,I)=WASH(4,I)-CON*CNTHETA*TERM(1,I) IWDOF(1,I)=KGRID(1,I) IWDOF(2,I)=KDOF(1,I) 100 CONTINUE NVFOUND=NVAL(1) GO TO 2000 1000 CONTINUE C C FIRST DOF IS YAW C DO 1100 I = 1,NVAL(1) WASH(1,I)=WASH(1,I)-CON1*CZPSI*TERM(1,I) WASH(2,I)=WASH(2,I)-CON*CMPSI*TERM(1,I) WASH(3,I)=WASH(3,I)-CON1*CYPSI*TERM(1,I) WASH(4,I)=WASH(4,I)-CON*CNPSI*TERM(1,I) IWDOF(1,I)=KGRID(1,I) IWDOF(2,I)=KDOF(1,I) 1100 CONTINUE NVFOUND=NVAL(1) 2000 CONTINUE C C CHECK SECOND DOF C IF(IDOF2.EQ.0)GO TO 9999 C C IDOF2 EXISTS C IF(IDOF2.EQ.IYP)GO TO 3000 C C IDOF2 IS PITCH C DO 2100 I = 1,NVAL(2) C C NEW DOF C NVFOUND=NVFOUND+1 I1=NVFOUND WASH(1,I1)=WASH(1,I1)-CON1*CZTHETA*TERM(2,I) WASH(2,I1)=WASH(2,I1)-CON*CMTHETA*TERM(2,I) WASH(3,I1)=WASH(3,I1)-CON1*CYTHETA*TERM(2,I) WASH(4,I1)=WASH(4,I1)-CON*CNTHETA*TERM(2,I) IWDOF(1,I1)=KGRID(2,I) IWDOF(2,I1)=KDOF(2,I) 2100 CONTINUE GO TO 9999 3000 CONTINUE C C IDOF2 IS YAW C DO 3100 I = 1,NVAL(2) C C NEW DOF C NVFOUND=NVFOUND+1 J1=NVFOUND WASH(1,J1)=WASH(1,J1)-CON1*CZPSI*TERM(2,I) WASH(2,J1)=WASH(2,J1)-CON*CMPSI*TERM(2,I) WASH(3,J1)=WASH(3,J1)-CON1*CYPSI*TERM(2,I) WASH(4,J1)=WASH(4,J1)-CON*CNPSI*TERM(2,I) IWDOF(1,J1)=KGRID(2,I) IWDOF(2,J1)=KDOF(2,I) 3100 CONTINUE C C DONE WITH IDOF2 C 9999 CONTINUE IF (IPRINT.GE.2)WRITE(6,9998) 9998 FORMAT(1H1,' DOWNWASH TERMS') if(iprint.ge.2) +WRITE(6,10000)(IWDOF(1,I),IWDOF(2,I),(WASH(J1,I),J1=1,4), + I=1,NVFOUND) 10000 FORMAT(' GRID ',I5,' DOF ',I2,' WASH=',4E14.6,/) iprint=ipold RETURN END C SUBROUTINE DMIG1(STIFF,DAMP,IGRID,IZP,IZPR, + IYP,IYPR,ICASE,ICONT,WASH,IWDOF,IDOF1,IDOF2,NVAL) c DIMENSION STIFF(4,4),DAMP(4,4),IDOF(4),IC(3),WASH(4,2000) DIMENSION STIFF(4,4),DAMP(4,4),IDOF(4),IC(5),WASH(4,2000) DIMENSION IWDOF(2,2000),NVAL(2) DIMENSION TEMP(4,2000),IGOT(4) INTEGER TDOF(2000),TGRID(2000) C C SUBROUTINE TO CREATE DMIG CARDS IN FILE 21 INCLUDING DOWNWASH TERMS C C C STIFFNESS TERMS FIRST C C WRITE TERMS INTO TEMPORARY MATRIX TEMP(I,J) C I=1--IZP;2--IYPR;3--IYP;4--IZPR C ICAS=ICASE+100 C C COMBINE WASG TERMS IN TEMP C IDOF(1) = IZP IDOF(2) = IYPR IDOF(3) = IYP IDOF(4) = IZPR NT=0 IT=0 IF(NVAL(1).EQ.0)GO TO 1100 DO 1000 I = 1,NVAL(1) TEMP(1,I) = WASH(1,I) TEMP(2,I) = WASH(2,I) TEMP(3,I) = WASH(3,I) TEMP(4,I) = WASH(4,I) TDOF(I) = IWDOF(2,I) 1000 TGRID(I) = IWDOF(1,I) NT = NVAL(1) C C FIRST SET OF TERMS PLACED IN TEMP C 1100 CONTINUE C C CHECK IF IDOF2 EXISTS C IF(IDOF2.EQ.0)GO TO 3000 C C IDOF2 EXISTS - PLACE TERMS IN TEMP C I1 = NVAL(1)+1 I2 = NVAL(1)+NVAL(2) DO 2000 I = I1,I2 C C CHECK FOR MATCH ON GRID AND DOF C ID1=0 1200 ID1=ID1+1 IF(ID1.GT.NT)GO TO 1500 IF(IWDOF(1,I).NE.TGRID(ID1))GO TO 1200 IF(IWDOF(2,I).NE.TDOF(ID1))GO TO 1200 C C MATCH ON GRID AND DOF C TEMP(1,ID1)=TEMP(1,ID1)+WASH(1,I) TEMP(2,ID1)=TEMP(2,ID1)+WASH(2,I) TEMP(3,ID1)=TEMP(3,ID1)+WASH(3,I) TEMP(4,ID1)=TEMP(3,ID1)+WASH(4,I) GO TO 1990 1500 CONTINUE C C NEW GRID AND DOF C NT = NT + 1 TEMP(1,NT)= WASH(1,I) TEMP(2,NT)= WASH(2,I) TEMP(3,NT)= WASH(3,I) TEMP(4,NT)= WASH(4,I) TDOF(NT) = IWDOF(2,I) TGRID(NT) = IWDOF(1,I) 1990 CONTINUE 2000 CONTINUE 3000 CONTINUE C C TEMP CONTAINS ALL INFO FROM WASH - NOW ADD STIFF C ID1 = 1 DO 4000 I = 1,NT IF(TGRID(I).NE.IGRID)GO TO 3990 IF(TDOF(I).EQ.IZP)IGOT(ID1)=1 IF(TDOF(I).EQ.IPR)IGOT(ID1)=2 IF(TDOF(I).EQ.IYP)IGOT(ID1)=3 IF(TDOF(I).EQ.IZPR)IGOT(ID1)=4 IF(IGOT(ID1).EQ.0) GO TO 3990 ID2=ID1 ID1=ID1+1 C C MATCH ON GRID AND DOF C TEMP(1,I) = TEMP(1,I) + STIFF(IGOT(ID2),1) TEMP(2,I) = TEMP(2,I) + STIFF(IGOT(ID2),2) TEMP(3,I) = TEMP(3,I) + STIFF(IGOT(ID2),3) TEMP(4,I) = TEMP(4,I) + STIFF(IGOT(ID2),4) C 3990 CONTINUE 4000 CONTINUE C C ID2 = NUMBER OF DOF MATCHED C IF(ID2.EQ.4)GO TO 5000 C C NEED TO APPEND STIFF TERMS ONTO TEMP C I1=4-ID2 DO 4500 I = 1,4 IDONT = 0 DO 4300 I3 = 1,I1 4300 IF(IGOT(I3).EQ.I)IDONT = I IF(IDONT.NE.0)GO TO 4400 C C APPEND TERMS C NT = NT + 1 TEMP(1,NT)=STIFF(I,1) TEMP(2,NT)=STIFF(I,2) TEMP(3,NT)=STIFF(I,3) TEMP(4,NT)=STIFF(I,4) TDOF(NT) = IDOF(I) TGRID(NT)=IGRID 4400 CONTINUE 4500 CONTINUE C C TEMP HAS ALL STIFFNESS TERMS C C C WRITE MATRIX TO FILE 21 AS DMIG CARDS C C C WRITE HEADER C WRITE(21,4600)ICAS,0,1,1 4600 FORMAT('DMIG AEROS',I3,3I8,/) DO 5000 I = 1,4 C C WRITE EACH SET OF TERMS C IC(1)=ICONT-1 IC(2)=ICONT WRITE(21,4650)ICAS,IGRID,IDOF(I),IC(2) 4650 FORMAT('DMIG *AEROS',I3,2I16,T73,'*S',I6) C C WRITE TERMS C DO 4700 I99=1,NT IC(1)=IC(2) IC(2)=IC(2)+1 4700 WRITE(21,4750)IC(1),TGRID(I99),TDOF(I99),TEMP(I,I99),IC(2) 4750 FORMAT('*S',I6,2I16,1PE16.9,16X,'*S',I6) ICONT = IC(2)+1 5000 CONTINUE C C STIFFNESS DMIG WRITTEN C C C WRITE DAMPING OUT C WRITE(21,5100)ICAS,0,1,1 5100 FORMAT('DMIG AEROB',I3,3I8,/) IDOF(1) = IZP IDOF(2) = IYPR IDOF(3) = IYP IDOF(4) = IZPR IC(1)=IC(1)+1 IC(2)=IC(1)+1 IC(3)=IC(1)+2 IC(4)=IC(1)+3 IC(5)=IC(1)+4 C DO 7000 I = 1,4 DO 5590 I99 = 1,5 5590 IC(I99) = IC(I99) + 5 7000 WRITE(21,7200)ICAS,IGRID,IDOF(I),IC(1), + (IC(J),IGRID,IDOF(J),DAMP(J,I),IC(J+1),J=1,4) 7200 FORMAT('DMIG *AEROB',I3,2I16,T73,'*D',I6,/, + 4('*D',I6,2I16,1PE16.9,16X,'*D',I6,/)) ICONT = IC(5) + 1 C C MATRICES WRITTEN C RETURN END