Iterative physics function

[Daniel89]

Hi, I am trying to make an iterative physics macro, but I have problems figuring out whats wrong. I have had success with a similar macro on another iterative macro earlier.

When I try to debug and run cursor to "P1" line, the formula goes blank and the debugger just stops...

Macro:

Function Navier_Stokes_compressibleflow_pressuredrop(SGavg, Dvisc, ID, T1, P2, T2, Lenght, PipeRoughness, scf_m, Zavg)
'March 2018 , function to calculate pressure drop for a compressible flow system with pure N2 gas
'INPUT:
'ID in mm
'T1 and T2 in celsius
'P2 in Bar
'Lenght in m
'Elevation in m
'Piperoughness in mm
'scf_m in standard cubic feet per minute
'Dvisc in Pa s
m = (scf_m / 35.3146667214886) / 60 ' m3/s
IDM = ID / 1000 ' ID in mm converted to m
T1K = T1 + 273.15 ' T1 in Celius converted to Kelvin
P2p = P2 * 10 ^ 5 'P2 in Bar converted to Pa
T2K = T2 + 273.15 ' T2 in Celius converted to Kelvin
R = 8.3145 'J/mol*K
Mw = 28.01348 'kg/kmol
Pi = 3.14159265358979
A = (IDM ^ 2 * Pi) / 4 'Area m2
Tavg = (T1K + T2K) * 0.5 ' average temperature Kelvin
Rho = SGavg * 1000
Relativeroughness = PipeRoughness / ID
P1new = 50 * 10 ^ 5 'initial guess
For i = 1 To 100
P1 = P1new
Pavg = (P1new + P2p) * 0.5
velocity = (m / A) / ((Pavg / (10 ^ 5)) / Zavg)
Nre = (Rho * velocity * IDM) / Dvisc ' Reynolds nr
Fd = 1.325 / (Log((Relativeroughness / 3.7) + (5.74 / (Nre ^ 0.9)))) ^ 2 'Darcy friction factor by Swamee Jain equation
P1 = (2 * ((R * Tavg) / Mw) * ((m / A) ^ 2) * ((Log(P1new / P2p)) * ((Fd * Lenght) / (2 * IDM))) - P2p ^ 2) ^ 0.5
Navier_Stokes_compressibleflow_pressuredrop = ((2 * ((R * Tavg) / Mw) * ((m / A) ^ 2) * ((Log(P1 / P2p)) * ((Fd * Lenght) / (2 * IDM))) - P2p ^ 2) ^ 0.5) / 10 ^ 5 'output in Bar
If Abs(P1new - P1) < 100 Then Exit Function
Next i
End
End
End Function

 

[Dryver14]

What exactly are you looping when you say i = 1 to 100, i is not mentioned anywhere in the loop that I can see.

In general (all the time?) a function returns one value.

Are you sure you are not running a sub routine. this line If Abs(P1new - P1) < 100 Then Exit Function should be placed earlier in the procedure ie: before the result of the function is given.

I have never seen a exit function statement before but I could be wrong on that point.

 

[Daniel89]

Hi you may be right,
I am looping to iterate P1, as the term is on both sides of the equation an interative procedure is needed. Pavg, velocity, Nre and Fd is also dependant on "new P1" for that loop, so they are included in the loop.

as I mentioned I had success with another macro that is working perfectly for me:

Function SRK_N2_Z(Bar, Celsius)
'Date: 08.02.2018
'Purpose: Calculate Z factor of N2 by Soavey Redlich Kwong Equation of state with iterations until error<10^-6
'Input:
'Pressure in Bar
'Temperature in Celsius
accFact = 0.04
Tc = 126.15 'K
Pc = 3397800 'Pa
R = 8.3145 'J/mol*K
TempK = Celsius + 273.15 'Kelvin
PressurePa = Bar * 10 ^ 5 'Pa
m = 28.01348 'kg/kmol
afactorN2 = 0.42748 * ((R ^ 2 * Tc ^ 2) / Pc)
bfactorN2 = 0.08664 * ((R * Tc) / (Pc))
Tr = TempK / Tc
Pr = PressurePa / Pc
alpha = (1 + (0.48508 + (1.55171 * accFact) - (0.15613 * accFact ^ 2)) * (1 - ((Tr) ^ (1 / 2)))) ^ 2
srkAfactor = (afactorN2 * alpha * PressurePa) / (R ^ 2 * TempK ^ 2)
srkBfactor = (bfactorN2 * PressurePa) / (R * TempK)
NewtonRhapsonZ = PressurePa / (R * TempK)
NewtonRhapsonA = (0.42748 * alpha * Pr) / (Tr ^ 2)
NewtonRhapsonB = (0.08664 * Pr) / Tr
Znew = 1 'initial guess
For i = 1 To 100
Z = Znew
F_Z = (Z ^ 3) - (Z ^ 2) + ((NewtonRhapsonA - NewtonRhapsonB - NewtonRhapsonB ^ 2) * Z) - (NewtonRhapsonA * NewtonRhapsonB)
D_Z = 3 * Z ^ 2 - 2 * Z + (NewtonRhapsonA - NewtonRhapsonB - NewtonRhapsonB ^ 2)
Znew = Z - (F_Z / D_Z)
SRK_N2_Z = Znew
If Abs(Znew - Z) < 10 ^ -6 Then Exit Function
Next i
End
End
End Function

 

[Dryver14]

Maybe someone else can see something in this I can't, but even in the function I cannon see the point of the loop.

If you comment out the loop do you get a different result?... The maths is beyond me by the way.

By the way in the first function, are you saying the formula goes blank when you say P1 = P1New

 

[Dryver14]

a follow up question is this line P2p = P2 * 10 ^ 5 'P2 in Bar converted to Pa referring to Range("P2") on the spreadsheet.

As you can see all the amounts for the equation are included in the function for your first function.. If you are obtaining numbers from cells then you need to reference them properly.

 

[offthelip]

I notice you have a variable called "lenght" which in incorrect spelling of Length but it is correctly used in both places
What I suggest is using the debug.print statement to print out all the variables of each equation just before you use them, this should give you values in the immediate window which hopefully will still be there when the debug stops.

 

[Daniel89]

Thx for showing interest

I dont understand your question: If you comment out the loop do you get a different result? can you elaborate?
In the first function the formula goes blank and goes out of the debugger mode when I come to line : P1 = (2 * ((R * Tavg) / Mw) * ((m / A) ^ 2) * ((Log(P1new / P2p)) * ((Fd * Lenght) / (2 * IDM))) - P2p ^ 2) ^ 0.5

P2 does not refer to a cell but "Pressure 2" P2 so I don't think the problem is there...

the equation is as follows:

as you see the P1(pressure1) which is the goal of the equation to find, is also on the right side of the equation. On these occations an iterative procedure is needed, every time you run the formula you will get closer to the correct answer. and you use the last P1 you got from the last iteration on the right side, until P1 does not change. for this procedure I have set the max allowed change to 100 Pascal, which is not much...

I may have to see if I can use another procedure, but I do not understand why it doesn't work. The last macro I posted ( Function SRK_N2_Z(Bar, Celsius)) works just fine, and from what I can see the procedure is the same...

 

[steve the fish]

Could you give some sensible values for the variables so i could test it? SGavg, Dvisc, ID, T1, P2, T2, Lenght, PipeRoughness, scf_m, Zavg

 

[Daniel89]

I don't know what happened but now the macro went through the first iteration, but then stepped out and did not repeat even though the absolute value of (P1new-P1) was more than 100.

the macro is this, and if someone wants to try it the input value cells is this:

SGavg: 0,057928051
Dvisc:1,97203*10^(-5)
ID: 49
T1: 20
P2: 1
T2: 20
Length: 10
Piperoughness: 0,1
Scf_m: 6000
Zavg: 1,00786436

Function Navier_Stokes_compressibleflow_pressuredrop(SGavg, Dvisc, ID, T1, P2, T2, Lenght, PipeRoughness, scf_m, Zavg)
'March 2018, function to calculate pressure drop for a compressible flow system with pure N2 gas
'INPUT:
'ID in mm
'T1 and T2 in celsius
'P2 in Bar
'Lenght in m
'Elevation in m
'Piperoughness in mm
'scf_m in standard cubic feet per minute
'Dvisc in Pa s
m = (scf_m / 35.3146667214886) / 60 ' m3/s
IDM = ID / 1000 ' ID in mm converted to m
T1K = T1 + 273.15 ' T1 in Celius converted to Kelvin
P2p = P2 * 10 ^ 5 'P2 in Bar converted to Pa
T2K = T2 + 273.15 ' T2 in Celius converted to Kelvin
R = 8.3145 'J/mol*K
Mw = 28.01348 'kg/kmol
Pi = 3.14159265358979
A = (IDM ^ 2 * Pi) / 4 'Area m2
Tavg = (T1K + T2K) * 0.5 ' average temperature Kelvin
Rho = SGavg * 1000
Relativeroughness = PipeRoughness / ID
P1new = 50 * 10 ^ 5 'initial guess
For i = 1 To 100
P1 = P1new
Pavg = (P1new + P2p) * 0.5
velocity = (m / A) / ((Pavg / (10 ^ 5)) / Zavg)
Nre = (Rho * velocity * IDM) / Dvisc
Fd = 1.325 / (Log((Relativeroughness / 3.7) + (5.74 / (Nre ^ 0.9)))) ^ 2 'Darcy friction factor by Swamee Jain equation
P1 = (2 * ((R * Tavg) / Mw) * ((m / A) ^ 2) * ((Log(P1new / P2p)) * ((Fd * Lenght) / (2 * IDM))) + P2p ^ 2) ^ 0.5 'to find new P1
Navier_Stokes_compressibleflow_pressuredrop = ((2 * ((R * Tavg) / Mw) * ((m / A) ^ 2) * ((Log(P1 / P2p)) * ((Fd * Lenght) / (2 * IDM))) + P2p ^ 2) ^ 0.5) / 10 ^ 5 'output in Bar
If Abs(P1new - P1) < 100 Then Exit Function
Next i
End
End Function

 

[wideboydixon]

Two things:

1. You never update the value for P1new so the iteration will never work and you'll hit the "End" statement which is bad
2. You have an "End" statement which is bad

I suggest you change this line:

P1 = (2 * ((R * Tavg) / Mw) * ((m / A) ^ 2) * ((Log(P1new / P2p)) * ((Fd * Lenght) / (2 * IDM))) + P2p ^ 2) ^ 0.5 'to find new P1

To this:

P1new = (2 * ((R * Tavg) / Mw) * ((m / A) ^ 2) * ((Log(P1new / P2p)) * ((Fd * Lenght) / (2 * IDM))) + P2p ^ 2) ^ 0.5 'to find new P1

At least then you should be able to get an answer. Using your parameters with an updated function I get an answer of:

1.00001632178924 (or 1,00001632178924)

WBD