I am trying to write a program to solve for pipe diameter for a pump system I've designed. I've done this on paper and in excel and understand the mechanics of the equations. I would appreciate any guidance.
At this point I am compiling but getting negative answers on the first iteration, my guess now is that it has to do with my initial guesses, the variables being passed to the function, or the secant calculation code itself. The three checks for phead, pump, and hf are correct.
MODULE Sec
CONTAINS
SUBROUTINE Secant(fx,xold,xnew,xolder)
IMPLICIT NONE
INTEGER,PARAMETER::DP=selected_real_kind(15)
REAL(DP), PARAMETER:: gamma=62.4
REAL(DP)::z,phead,hf,L,Q,mu,rho,rough,eff,pump,nu,ppow,fric,pres,xnew,xold,xolder,D
INTEGER::I,maxit
INTERFACE
FUNCTION fx(D,L,Q,hf,rho,mu,rough)
Omitted
END FUNCTION
END INTERFACE
Q=0.0353196
Pres=-3600.0
z=-10.0
L=50.0
mu=0.0000273
rho=1.940
nu=0.5
rough=0.000005
ppow=412.50
xold=1.0
xolder=0.5
phead = (pres/gamma)
WRITE(*,*) phead
pump = (nu*ppow)/(gamma*Q)
WRITE(*,*)pump
hf = phead + z + pump
WRITE(*,*)hf
maxit=20
I = 1
DO
xnew=xold-fx(xold,L,Q,hf,rho,mu,rough)*((xolder-xold)/(fx(xolder,L,Q,hf,rho,mu,rough)-fx(xold,L,Q,hf,rho,mu,rough)))
xolder = xold
xold = xnew
I=I+1
WRITE(*,*) "Diameter = ", xnew
IF (ABS(fx(xnew,L,Q,hf,rho,mu,rough)) <= 0.05) THEN
EXIT
END IF
IF (I >= maxit) THEN
EXIT
END IF
END DO
END SUBROUTINE Secant
END MODULE Sec
PROGRAM Pipes
USE Sec
IMPLICIT NONE
INTEGER,PARAMETER::DP=selected_real_kind(15)
REAL(DP)::xold,xolder,xnew
INTERFACE
FUNCTION f(D,L,Q,hf,rho,mu,rough)
Omitted
END INTERFACE
CALL Secant(f,xold,xnew,xolder)
END PROGRAM Pipes
FUNCTION f(D,L,Q,hf,rho,mu,rough)
IMPLICIT NONE
INTEGER,PARAMETER::DP=selected_real_kind(15)
REAL(DP), PARAMETER::pi=3.14159265, g=9.81
REAL(DP), INTENT(IN)::L,Q,rough,rho,mu,hf,D
REAL(DP)::f
f=(1/(hf/((L/D)*((4*Q)/pi*D))))+2.0*log((rough/(3.7*D))+(2.51/(((rho*((4*Q)/pi*D))/mu)*(hf/((L/D)*((4*Q)/pi*D))))))
END FUNCTION
