wing Translation Listing



 FORTRAN CALCULUS (VERSION 2.0) ------------------------- 04/12/2008   07:52:34
 ####  LABEL ---+---------+---------+---------+---------+---------+---------+--
~
    1        GLOBAL ALL                      ! Wing Design Optimization
    2        PROBLEM WING (40000,5000)
    3          DYNAMIC XS,Y,ALPHA,ALPHAL,ALPHAU,DELALP,EIS,EIL,EIU,DELEI
    4          @NODIAG(1006)
    5          CALL SETUP
    6          FIND EIS,ALPHA,DYF; ! Flexural rigidity vector, length fraction
           ~    IN FLEX; BY THOR(TCON);
           ~    WITH BOUNDS DELEI,DELALP,DYB;
           ~    WITH LOWERS EIL,ALPHAL,DYL;
           ~    WITH UPPERS EIU,ALPHAU,DYU;
           ~    HOLDING WEIGHT; MATCHING ALTOTL;
           ~    TO MAXIMIZE WINGLIFT            ! Lift/Weight
    7          @PAGESUM(’Final Wing Curvature’)
    8          PRINT *,’ SOLUTION TO ODE-S’
    9          TABULATE XS,Y
   10        END
~
      
   11        PROCEDURE SETUP
   12          PI=3.141592654
   18>         O=0: WL=1: A=1: EPSLN=0.1: DYF=0.175: DYLDX=0 ! 10 deg, 0 deg
   21>         H=0.05: NEI=5: WTMAX=7
   22*         NPTS=(WL-O)/H+1
   23*         ALLOT XS(NPTS),Y(NPTS)
   24          ALLOT ALPHA(NEI),ALPHAL(NEI),ALPHAU(NEI),ALPHAU(NEI),DELALP(NEI)
   25          ALLOT EIS(NEI),EIL(NEI),EIU(NEI),DELEI(NEI)
   28>         <ALPHA>=(0.2): <ALPHAL>=(0.1): <ALPHAU>=(1.0)
   30>         <DELALP>=(0.0333): <DELEI>=(0.333)
   32>         <EIS>=DATA(6,5,4,3,2): <EIL>=DATA(3,2,2.5,1,1)
   33          <EIU>=DATA(10,7,5,4,4)
   34          DYL=0
   35          DYU=.2
   36          DYB=.025
   37          DO 10 I=1,NPTS
   38            XS(I)=O+(I-1)*H
   39            Y(I)=EPSLN*SIN(PI*XS(I)/2)
   40    10    CONTINUE
   41*         @STREAM(’SUMMARY’)
   42>         @PAGESUM(’Constant constraints: EIL,EIU,ALPHAL and ALPHAU’)
   43          TABULATE EIL,EIU,ALPHAL
   44        END
~
      
      
   45        MODEL FLEX
   46          DY0=DYF
   47          FIND DY0; IN BEAMIVP;  ! Beam equation boundary value problem
           ~    BY AJAX(KNOB); TO MATCH YPA
   48          EIAVG=0
   49          DO 10 K=1,NEI
   50    10    EIAVG=EIAVG+ALPHA(K)*EIS(K) ! Avg. flexural rigidity
   51          LFLG=1
   52*         WINGLIFT=$INTEGGL(FORCE,O,WL,4)/EIAVG ! Gauss quadrature integra
   53          LFLG=0
   54*         WEIGHT=WTMAX-EIAVG       ! Weight (inequality) constraint
   55          ALTOTL=ARRAYSUM(ALPHA)-1 ! Total length (equality) constraint
   56          @PAGESUM(’Intermediate Optimization Results’)
   57          ROWPRINT EIAVG,WINGLIFT,WEIGHT,ALTOTL
   58          TABULATE EIS,ALPHA
   59        END
~
      
   60        MODEL BEAMIVP               ! Beam equation initial value problem
   64>         X=O: Y0=0: DY0DX=DY0: XF=H
   65          Y(1)=Y0
   66          INITIATE ATHENA(CNTRL); FOR BEAMODE; ! Integrate beam equation
           ~    EQUATIONS D2Y0DX2/DY0DX,DY0DX/Y0; OF X; STEP H; TO XF
   67          DO 20 I=2,NPTSC
   68            ALTOTL=-1
   69            DO 10 LL=1,NEI
   70              ALTOTL=ALTOTL+ALPHA(LL)
   71     10     CONTINUE
   72*           INTEGRATE BEAMODE; BY ATHENA
   73            Y(I)=Y0
   74            XF=XF+H
   75     20   CONTINUE
   76*         TERMINATE BEAMODE
   77          YPA=DY0DX-DYLDX  ! Boundary constraint - must be zero
   78        END
~
      
   79        MODEL BEAMODE  ! Cantilever beam differential equation
   80          SUM=0
   81          DO 10 J=1,NEI
   82            IF(X.GE.SUM*WL.AND.X.LT.(ALPHA(J)+SUM)*WL) THEN
   83              EI=EIS(J)
   84              GOTO 20
   85            ENDIF
   86            SUM=SUM+ALPHA(J)
   87    10    CONTINUE
   88*   20    D2Y0DX2=-$INTEGGL(FORCE,O,X,4)/EI*(1+DY0DX**2)**1.5
   89        END
~
      
   90        FMODEL FORCE(R) ! Quadrature integrand - lift or moment
   91          IF(LFLG.GT.0) THEN  ! Calculate lift integrand
   92            FORCE=A*COS(PI*R/(1.95*WL))-EPSLN*TERP(R)
   93          ELSE  ! Calculate moment integrand
   94*           FORCE=R*(A*COS(PI*R/(1.95*WL))-EPSLN*TERP(R))
   95          ENDIF
   96        END
~
      
   97        FMODEL TERP(Z) ! Lagrange interpolation of wing formula
   98          IF(Z.EQ.XS(1)) THEN
   99            TERP=Y(1)
  100          ELSEIF(Z.EQ.XS(NPTS)) THEN
  101            TERP=Y(NPTS)
  102          ELSEIF(Z.LT.XS(2)) THEN
  103            TERP=POLY(XS(1),XS(2),XS(3),Y(1),Y(2),Y(3),Z)
  104          ELSE
  105*           DO 10 II=2,NPTS-1
  106              IF(Z.GE.XS(II)) GOTO 20
  107    10      CONTINUE
  108*   20      TERP=POLY(XS(II-1),XS(II),XS(II+1),Y(II-1),Y(II),Y(II+1),Z)
  109          ENDIF
  110        END
~
      
  111        FMODEL POLY(XI,XJ,XK,YI,YJ,YK,ZZ)
  112          GI=(ZZ-XJ)*(ZZ-XK)/((XI-XJ)*(XI-XK))
  113          GJ=(ZZ-XI)*(ZZ-XK)/((XJ-XI)*(XJ-XK))
  114          GK=(ZZ-XI)*(ZZ-XJ)/((XK-XI)*(XK-XJ))
  115          POLY=GI*YI+GJ*YJ+GK*YK
  116        END
~
      
  117        CONTROLLER TCON(THOR)
  118          REMAX=12
  119          SUMOUT=1
  120          DETAIL=1
  121          DETOUT=0
  122        END
~
      
  123        CONTROLLER KNOB(AJAX)
  124          REMAX=10
  125          DAMP=0
  126          ZERO=H**3
  127          DETAIL=1
  128          DETOUT=0
  129          SUMOUT=1
  130        END
~
      
  131        CONTROLLER CNTRL(ATHENA)
  132          ORDER=2
  133        END