管道直径的割线法求解

Question

我正试图编写一个程序来解决我设计的泵系统的管径问题。我已经在纸上做过了，并且了解了方程的力学原理。如果有任何指导，我将不胜感激。

编辑：我已经用用户的一些建议更新了代码，仍然看到了快速的分歧。那里的猜测太高了。如果我弄清楚了，我会把它更新为工作状态。

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

Answer 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（第二个伪参数）获得任何参数，因为它会在第一个错误时停止。

Answer 2

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)的表达式总是很有用的。（一个技巧是提取外部的常量因子并预先计算它们。）

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