Welcome to AE Resources
Converted document Converted document

FORTRAN 90+: DECLARATION OF VARIABLES DEMONSTRATION

In this demonstration, the use of DECLARATION OF VARIABLES is demonstrated.
In fluid mechanics or aerodynamics, it is possible to determine the local Mach number at a point when the pressure coefficient and free stream Mach numbers are known. This can be important when trying to determine if a point on a configuration is in the subsonic or supersonic Mach regimes. The equation to determine the pressure coefficient at point a, Cp, a is nonlinear:
(1) Cp, a = (2)/(γM2∞)⎡⎣⎛⎝(1 + (γ − 1) ⁄ 2M2∞)/(1 + (γ − 1) ⁄ 2M2a)⎞⎠γ ⁄ (γ − 1) − 1.0⎤⎦
M∞ and Cp, a are known, while Ma is not known. A bisection or Newton-Raphson method can be used to determine Ma through the process known as root finding. 
PROGRAM Cp_Mach
	IMPLICIT NONE
	REAL, PARAMETER :: pi=4.0*ATAN(1.0), gamma = 1.4
	INTEGER :: choice, nmax
	LOGICAL :: mystop = .TRUE.
	REAL :: M, Cp, Ma, ml, mr, tol, Cp_test, Cp_test2, slope
	
	WRITE(*,*) ’This program computes the Mach number at a certain point on an airfoil.’
	WRITE(*,*) ’The output should be within .01% accuracy.’
	WRITE(*,*) ’Required user input: Cp value and freestream Mach number.’
	WRITE(*,*) ’ ’
	
	!Input Mach number
	DO WHILE (mystop)
	     WRITE(*,*) ’Enter Mach (M-Infinity) number’
	     READ(*,*) M
	     ! Make sure the Mach number  is reasonable
	     IF ( M <= 0.0 .OR. M > 10.0) THEN
	        mystop = .TRUE.
	         WRITE(*,*) ’Too extreme of a Mach number.’
	         WRITE(*,*) ’Please choose a Mach number between [0.1, 10].’
	     ELSE
	         mystop = .FALSE.
	     END IF
	END DO
	
	!Input Cp value
	DO
	  WRITE(*,*) ’Enter coefficient of pressure (Cp) value’
	  READ(*,*) Cp
	  ! Make sure Cp is reasonable
	  IF (Cp <= 1.0 .AND. Cp >= -10.0) THEN
	     EXIT
	  ELSE
	     WRITE(*,*) ’ Cp values should be between [-10,1] ’
	  END IF	
	END DO
	!Set a few constants
	gam = gamma
	gm1 = gamma - 1.0
	
	!Input numerical method choice
	DO
	    WRITE(*,*) ’Enter <1> for Bisection Method or <2> for Newton-Raphson Method’
	    READ(*,*) choice
	    IF (choice == 1 .OR. choice == 2) THEN
	         EXIT
	    ELSE
	         WRITE(*,*) ’Please enter an appropriate choice (1 or 2)’
	   END IF
	END DO	
	! nmax = 0  !Make sure there are no infinite loops.
	SELECT CASE (choice)
	CASE (1)
	     WRITE(*,*) ’Bisection Method’
	     ! Set bounds
	     ml = 0.0
	     mr = 10.0
	     DO
	     ! initialize Ma at middle of bounds
	          Ma = (ml + mr) / 2.0
	          Cp_test = (2.0/(gam*M**2)) * (( (1. + (gm1/2.)*(M**2.)) &
	                / (1. + (gm1/2.)*(Ma**2.)) ) ** (gam/gm1) - 1.) - Cp
	           ! Within tolerance?
	            tol = ABS(Cp_test)
	            nmax = nmax + 1
	            IF (tol <= .0001) THEN !  Within tolerance of .1%
	                 EXIT
	            ELSE IF (nmax >= 500) THEN !Stop infinite loop
	                 WRITE(*,*) ’ Maximum iterations reached (500).’
	                 WRITE(*,*) ’Error Analysis:’
	                 WRITE(*,*) ’Either Mach or Cp may be too extreme, causing’
	                 WRITE(*,*) ’Ma to never be found.  Please check input and retry.’
	                 STOP
	            ELSE IF (Cp_test > 0.0) THEN
	                  ! Update the left bound
	                  ml = Ma
	            ELSE IF (Cp_test < 0.0) THEN
	                   ! Update the right bound
	                   mr = Ma
	             END IF
	       END DO
	CASE (2)
	       WRITE(*,*) ’Newton-Raphson Method’
	       Ma = .9*M ! initial guess
	       DO
	          Cp_test = (2.0/(gam*M**2)) * (( (1. + gm1/2.)*(M**2.)) &
	                  / (1. + (gm1/2.)*(Ma**2.)) ) ** (gam/gm1) - 1.) - Cp
	           tol = ABS(Cp_test)
	           nmax = nmax + 1
	           IF (tol <= .001) THEN !Within tolerance of .1%
	             EXIT
	           ELSE IF (nmax >= 500) THEN !Stop inifinite loop
	              WRITE(*,*) ’ Maximum iterations reached (500).’
	              WRITE(*,*) ’Error Analysis:’
	              WRITE(*,*) ’Either Mach or Cp may be too extreme, causing’
	              WRITE(*,*) ’Ma to never be found.  Please check input and retry.’
	              STOP
	           ELSE
	              ! Obtain local slope 
	              Ma = Ma + .00001
	              Cp_test2 = (2.0/(gam*M**2)) * (( (1. + (gm1/2.)*(M**2.)) &
	                      / (1. + (gm1/2.)*(Ma**2.)) ) ** (gam/gm1) - 1.) - Cp
	              slope = (Cp_test2 - Cp_test) / 0.00001
	              ! Apply formula M(i+1) = M - f(M) / f’(M)
	              Ma = (Ma - 0.000001) - (Cp_test) / slope
	            END IF
	          END DO
	END SELECT
WRITE(*,*) ’Ma =’, Ma	
END PROGRAM Cp_Mach