- #1

- 9

- 0

REAL :: l !LENGTH OF THE BAR

REAL :: temp !TEMPERATURE IN KELVIN

REAL :: d_0 !CONSTANT DIFFUSION PARAMETER

REAL :: q !ACTIVATION ENERGY

REAL :: c_0 !INITIAL CONCENTRATION AT X=0

REAL :: r=8.31 !IDEAL GAS CONSTANT IN J/molK

REAL :: c_r !DESIRED END CONCENTRATION

INTEGER :: imax=30 !MAXIMUM ITERATIONS

REAL :: eps=1E-4 !MAXIUMUM ERROR

REAL :: a, b !INTERVAL OF EVALUATION [A,B]

REAL :: t

REAL :: exp, df !, erf, c, c_avg

END MODULE space_data

PROGRAM carburization

USE space_data

IMPLICIT NONE

CHARACTER :: choice

REAL :: f, root

INTEGER :: ercode

DO

WRITE(*,*) 'a) Enter new data'

WRITE(*,*) 'b) Calculate time until desired concentration'

WRITE(*,*) 'q) Quit'

WRITE(*,'(A,\)') 'Enter choice: '

READ(*,*) choice

SELECT CASE(choice)

CASE('A','a')

CALL call_data()

CASE('B','b')

CALL bisection(f, root, ercode)

CALL results()

CASE('Q','q')

STOP

END SELECT

END DO

CONTAINS

SUBROUTINE call_data()

USE space_data

IMPLICIT NONE

WRITE(*,'(A,\)') 'Enter initial end concentration (%): '

READ(*,*) c_0

WRITE(*,'(A,\)') 'Enter diffusion parameter (m^3/sec): '

READ(*,*) d_0

WRITE(*,'(A,\)') 'Enter activation energy (J): '

READ(*,*) q

WRITE(*,'(A,\)') 'Enter temperature (K): '

READ(*,*) temp

WRITE(*,'(A,\)') 'Enter length of bar (m): '

READ(*,*) l

WRITE(*,'(A,\)') 'Enter desired average final concentration (%): '

READ(*,*) c_r

WRITE(*,*) ''

WRITE(*,*) ' The bisection method needs an initial search interval [A,B].'

WRITE(*,*) ''

DO

WRITE(*,'(A,\)') 'Enter a A (Years):'

READ(*,*) a

WRITE(*,'(A,\)') 'Enter a B (Years):'

READ(*,*) b

IF(a < b) THEN

EXIT

ELSE

WRITE(*,*) 'Error- left end point needs to be smaller than right end point.'

WRITE(*,*) 'Please enter again:'

END IF

END DO

END SUBROUTINE call_data

END PROGRAM carburization

SUBROUTINE bisection(f, root, ercode)

USE space_data

IMPLICIT NONE

INTEGER::I

REAL,EXTERNAL::F

REAL,INTENT(OUT)::ROOT !THE ROOT WE ARE TRYING TO FIND

REAL::X1,X2,X3,F1,F2,F3,D,D01 !X VALUES, F(X) VALUES AND INTERVAL LENGTHS

INTEGER,INTENT(OUT)::ERCODE !ERROR CODE WILL RETURN AN INT VALUE BASED ON OUTCOME

!-1 MEANS ITERATION CNT EXCEED

!-2 MEANS NO ROOT IN INTERVAL A-B

!-3 UNKNOWN ERROR

! 0 ROOT FOUND, BISECTION METHOD SUCCESSFUL

X1=A

X3=B

F1 = F(X1)

F3 = F(X3)

IF(F1*F3>0) THEN

ERCODE = -2 !-2 MEANS NO ROOT IN INTERVAL A-B

ELSE

D=1

I=0

D01=A-B

DO

X2 = (X1+X3)/2.

F2 = F(X2)

IF(D<EPS) THEN

ROOT = X2

ERCODE = 0 !SUCCESS

EXIT

ELSE IF(I>IMAX) THEN

ERCODE = -1 !-1 MEANS ITERATION CNT EXCEED

EXIT

END IF

IF(F1*F2<0) THEN !LEFT

D = (X2-X1)/D01

F3 = F2

X3 = X2

ELSE IF(F2*F3<0)THEN!RIGHT

D = (X3-X2)/D01

F1 = F2

X1 = X2

ELSE IF(F2 == 0) THEN

ROOT = X2

ERCODE = 0 !SUCCESS

EXIT

ELSE

ERCODE = -3 !-3 UNKNOWN ERROR

EXIT

END IF

I = I+1

END DO

END IF

END SUBROUTINE

SUBROUTINE results()

USE space_data

IMPLICIT NONE

WRITE(*, '(T5,A)') 'RESULTS:'

WRITE(*, '(T5, 50("-"))')

WRITE(*,'(T5, "|", 1X, A,T40, "|", F12.3, 1X, "|")') 'Initial end concentration (%)',c_0

WRITE(*,'(T5, "|", 1X, A,T40, "|", ES12.3, 1X, "|")') 'Diffusion parameter (m^3/sec):',d_0

WRITE(*,'(T5, "|", 1X, A,T40, "|", ES12.3, 1X, "|")') 'Activationenergy (J):',q

WRITE(*,'(T5, "|", 1X, A,T40, "|", F12.2, 1X, "|")') 'Temperature (K):',temp

WRITE(*,'(T5, "|", 1X, A,T40, "|", F12.2, 1X, "|")') 'Length of bar (m):',l

WRITE(*,'(T5, "|", 1X, A,T40, "|", F12.3, 1X, "|")') 'Final concentration (%):',c_r

WRITE(*, '(T5, 50("-"))')

WRITE(*,'(T5, "|", 1X, A,T40, "|", F12.3, 1X, "|")') 'Computed diffusion time (days):',t

WRITE(*, '(T5, 50("-"))')

WRITE(*,*)

END SUBROUTINE results

REAL FUNCTION f()

USE space_data

IMPLICIT NONE

f = c_avg - c_r

END FUNCTION f

REAL FUNCTION c ()

USE space_data

IMPLICIT NONE

REAL :: z, x

df = d_0*(exp**(-q/(r*temp)))

x=l

z = x/(2*df*t)

c = (c_0/2.)*(1 - ERF(z))

END FUNCTION c

REAL FUNCTION c_avg ()

USE space_data

IMPLICIT NONE

REAL :: delta_x, idx

INTEGER :: i

delta_x = .2

DO i = 0, 5

idx = i*delta_x

c_avg = c_avg + (1/6.)*( c(idx,t) )

END DO

END FUNCTION c_avg