管径正割法求解

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个参数给它。其中 none 是可选的。

但是当你称呼它的时候，你称呼它为

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) 的表达式尽可能简单地排列总是有用的。（一个技巧是在外部提取常数因子并预先计算它们。）

此外，出于调试目的，有时检查各种术语的物理尺寸和物理单位的一致性有时很有用。事实上，如果得到的解的量级太大或太小，可能会存在一些物理单位换算系数的问题，例如。

管径正割法求解

Secant method solving for pipe diameter

fortran

numerical-methods

fluid-dynamics