0

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 understand the mechanics of the equations. I would appreciate any guidance.

EDIT: I have updated the code with some suggestions from users, still seeing quick divergence. The guesses in there are way too high. If I figure this out I will update it to working.

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
omitted
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.90
D=11.0
phead = (pres/gamma)
pump = (nu*ppow)/(gamma*Q)
hf = phead + z + pump

maxit=10
I = 1

DO
xnew=xold-((fx(xold,L,Q,hf,rho,mu,rough)*(xold-xolder))/ &
      (fx(xold,L,Q,hf,rho,mu,rough)-fx(xolder,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)) <= 1.0d-10) THEN
EXIT
END IF

IF (I >= maxit) THEN
EXIT
END IF 
END DO

RETURN

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
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.14159265d0, g=9.81d0
REAL(DP), INTENT(IN)::L,Q,rough,rho,mu,hf,D
REAL(DP)::f, fric, reynold, coef

fric=(hf/((L/D)*(((4.0*Q)/(pi*D**2))/2*g)))

reynold=((rho*(4.0*Q/pi*D**2)*D)/mu)

coef=(rough/(3.7d0*D))

f=(1/SQRT(fric))+2.0d0*log10(coef+(2.51d0/(reynold*SQRT(fric))))

END FUNCTION
4

2 回答 2

1

您非常清楚地将接口（和实现）中的函数声明为

FUNCTION f(L,D,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)::fx
END FUNCTION

所以你需要向它传递 7 个参数。它们都不是可选的。

但是当你调用它时，你称它为

xnew=xold-fx(xold)*((xolder-xold)/(fx(xolder)-fx(xold))

为其提供一个参数。例如，当您尝试编译它时gfortran，编译器会抱怨没有为D（第二个虚拟参数）获取任何参数，因为它会因第一个错误而停止。

于 2015-12-14T00:13:44.470 回答
0

似乎 和 的初始值xold离xolder解太远了。如果我们将它们更改为

xold   = 3.0d-5
xolder = 9.0d-5

并将收敛阈值更严格地更改为

IF (ABS(fx(xnew,L,Q,hf,rho,mu,rough)) <= 1.0d-10) THEN

然后我们得到

...
Diameter =    7.8306011049894322E-005
Diameter =    7.4533171406818087E-005
Diameter =    7.2580746283970710E-005
Diameter =    7.2653611474296094E-005
Diameter =    7.2652684750264582E-005
Diameter =    7.2652684291155581E-005

在这里，我们注意到函数f(x)定义为

FUNCTION f(D,L,Q,hf,rho,mu,rough)
...
f = (1/(hf/((L/D)*((4*Q)/pi*D))))                                   !! (1)
    + 2.0 * log(  (rough/(3.7*D)) + (2.51/(((rho*((4*Q)/pi*D))/mu)  !! (2)
                    * (hf/((L/D)*((4*Q)/pi*D)))))                   !! (3)
               )
END FUNCTION

其中第 (1) 行和 (3) 行中的项都是常数，而第 (2) 行中的项是 上的一些常数D。所以，我们看到f(D) = c1 - 2.0 * log( D / c2 )，所以我们可以解析得到解D = c2 * exp(c1/2.0) = 7.26526809959e-5，这与上面的数值解非常吻合。要大致了解解决方案在哪里，将其绘制f(D)为 的函数很有用D，例如使用 Gnuplot。

但我担心表达式f(D)本身（在 Fortran 代码中给出）可能会由于许多括号而包含一些拼写错误。为避免此类问题，在编写程序之前首先将表达式安排得f(D)尽可能简单总是很有用的。（一个提示是提取外部常数因子并预先计算它们。）

此外，出于调试目的，有时检查物理尺寸和各种术语的物理单位的一致性很有用。事实上，如果得到的解的幅度太大或太小，例如，物理单位的转换因子可能会出现一些问题。

于 2015-12-14T04:24:45.227 回答