Решение осевого метода для диаметра трубы

Question

Я пытаюсь написать программу для решения проблемы диаметра трубы для насосной системы, которую я разработал. Я сделал это на бумаге и понял механику уравнений. Буду признателен за любые рекомендации.

EDIT: Я обновил код с некоторыми предложениями от пользователей, все еще наблюдая быстрое расхождение. Догадки там слишком высокие. Если я это выясню, я обновлю его на работу.

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

Похоже, что начальные значения для 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) как можно более простым, прежде чем делать программу. (Один TIP должен извлекать постоянные факторы снаружи и предварительно вычислять их.)

Кроме того, для целей отладки иногда полезно проверять согласованность физических измерений и физических единиц различных терминов. Действительно, если величина полученного решения слишком велика или слишком мала, может быть, например, проблема коэффициентов пересчета для физических единиц.

Answer 2

Вы очень четко объявить функцию в интерфейсе (и реализации), как

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 (второй фиктивный аргумент), потому что он останавливается с первой ошибкой.