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

 

[steve the fish]

Using the supplied variables and equation P1 iterates down from 100739.9 to 100000. P1new is 500000 so Abs(P1new - P1) < 100 is never true. Is this a likely outcome?

 

[Daniel89]

Thx a lot wideboydixon The macro is ALIVE and working perfectly now, what a stupid mistake, I am novice in VBA programming...

Another question:

I have 2 macros, I want this macro1 to calculate a value within macro2:

macro1:

Function SRK_N2_Density(Bar, Celsius)
'Date: 08.02.2018, 'Purpose: Calculate SG of N2 with iterations until error<10^-6 by the Soavey Redlich Kwong Equation of State
'Purpose: Calculate SG of N2 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_Density = ((m * 10 ^ -3) / ((Znew * R * TempK) / (PressurePa / 1000))) * 1000 'output in kg/m3
If Abs(Znew - Z) < 10 ^ -6 Then Exit Function
Next i
End
End
End Function

macro2:

Function Navier_Stokes_compressibleflow_Pressuredrop(SGavg, Dvisc, ID, T1, P2, T2, Lenght, PipeRoughness, scf_m, Zavg)
'March 2018 made by Daniel Bjerkmo, 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
'Piperoughness in mm
'scf_m in standard cubic feet per minute
'Dvisc in Pa s
m = (scf_m / 35.3146667214886) / 60
Tavg = (T1K + T2K) * 0.5 ' average temperature Kelvin

m_kg_s = m * 1.14941818065257 ' kg/s ( HERE: I want to call on macro1 to calculate the number 1.14941818065257(standard value at 20°C) which is the density of 1sm3 of nitrogen gas, I want it to calculate on the basis of 1 bar and the Tavg elaborated in macro2 on the line above, to get precise results)

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
Rho = SGavg * 1000
Relativeroughness = PipeRoughness / ID
P1new = 50 * 10 ^ 5 'initial guess
For i = 1 To 1000
P1 = P1new
Pavg = (P1new + P2p) * 0.5 ' Average pressure
velocity = (m / A) / ((Pavg / (10 ^ 5)) / Zavg) 'average velocity in m/s for the system based on the average pressure and compressibility factor(Z)
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
P1new = (2 * ((R * Tavg) / Mw) * ((m_kg_s ^ 2) / (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_kg_s ^ 2) / (A ^ 2)) * ((Log(P1 / P2p)) * ((Fd * Lenght) / (2 * IDM))) + P2p ^ 2) ^ 0.5) - P2p) / 10 ^ 5 'output in Bar
If Abs(P1new - P1) < 10 Then Exit Function
Next i
End
End Function

I have tried several ways, but I cannot get it to work...

 

[Daniel89]

Thank you steve, the problem with the macro was explained by Wideboydixon. Working fine now Thx a lot for your interest

 

[wideboydixon]

Is it not just:

m_kg_s = m * SRK_N2_Density(1, Tavg) ' kg/s

?

WBD

 

[Daniel89]

Thank you very much Wideboydixon worked like a charm I am very new to VBA coding, I mostly use it for engineering purposes, it is not that easy to find help online for those purposes...

THANK YOU

 

[wideboydixon]

Glad I could help. Your

End

Statements still upset me though!

WBD

 

[Daniel89]

Hi again, I am struggling with something.

There is a certain range where the formula is valid.

this is my macro:

Function Navier_Stokes_compressibleflow_Pressuredrop_ID(PipeID_mm, T1, P2, T2, Pipelength_m, PipeRoughness_mm, Flowrate_scf_min)
'March 2018, function to calculate pressure drop for a compressible flow system with pure N2 gas, by the use of the Navier stokes momentum equation
'INPUT:
'PipeID-mm in milimeter
'T1 and T2 in celsius
'P2 in Bar
'Pipelength_m in m
'PipeRoughness_mm in mm
'Flowrate_scf_min in standard cubic feet per minute
M = (Flowrate_scf_min / 35.3146667214886) / 60 'm3/s
m_kg_s = M * ((SRK_N2_Density(1.01325, T1) + SRK_N2_Density(1.01325, T2)) * 0.5) ' kg/s
IDM = PipeID_mm / 1000 'PipeID_mm converted to m
Tavg = (T1 + T2) * 0.5
T1K = T1 + 273.15 ' T1 in Celius converted to Kelvin
T2K = T2 + 273.15 ' T2 in Celius converted to Kelvin
TavgK = (T1K + T2K) * 0.5 ' average temperature Kelvin
P2p = P2 * 10 ^ 5 'P2 in Bar converted to Pa
R = 8.3145 'J/mol*K
Mw = 0.02801348 'kg/mol
Pi = 3.14159265358979
g = 9.81 'm2/s2
Elevation = 300 ' m
A = (IDM ^ 2 * Pi) / 4 'Area m2
Relativeroughness = PipeRoughness_mm / PipeID_mm
P1new = 1 * 10 ^ 5 'initial guess
For i = 1 To 1000
P1 = P1new
Pavg = (P1new + P2p) * 0.5 ' Average pressure
D_visc_Pas = (N2_Viscosity_Dynamic((P1 * 10 ^ -5), T1) + N2_Viscosity_Dynamic(P2, T2)) * 0.5
SGavg = (SRK_N2_SG((P1 * 10 ^ -5), T1) + SRK_N2_SG(P2, T2)) * 0.5
velocity = (M / A) / ((Pavg / (101325)) / ((SRK_N2_Z((P1 * 10 ^ -5), T1) + SRK_N2_Z(P2, T2)) * 0.5)) 'average velocity in m/s for the system based on the average pressure and compressibility factor(Z)
Fd = Haaland_Equation_Darcy_frictionfactor(PipeRoughness_mm, PipeID_mm, SGavg, velocity, D_visc_Pas) 'Haaland Friction factor
P1new = (((2 * R * TavgK * m_kg_s ^ 2) / (Mw * A ^ 2)) * (Log(P1 / P2p) + ((Fd * Pipelength_m) / (2 * IDM))) + P2p ^ 2) ^ 0.5 'to find new P1
Navier_Stokes_compressibleflow_Pressuredrop_ID = (P1 / 10 ^ 5) - P2 'output in Bar
eq_14_1 = (Abs(SGavg * 1000 * g * Elevation) / ((Fd * Pipelength_m * m_kg_s ^ 2) / (2 * A ^ 2 * IDM * SGavg * 1000)))
eq_14_2 = (SGavg * 1000 * g) / ((Fd * m_kg_s ^ 2) / (2 * A ^ 2 * IDM * SGavg * 1000))
eq_15 = (((Fd) / (2 * IDM)) * (m_kg_s / A)) ^ 0.5
error_horisontal = -((Pi ^ 2 * IDM ^ 6 * Mw * (P1 ^ 2 - P2p ^ 2)) / ((Fd * Pipelength_m) ^ 2 * R * TavgK * (((2 * IDM) / (Fd * Pipelength_m)) + (1 / (Log(P1 / P2p))))))
If Abs(P1new - P1) < 10 Then Exit Function
Next i
End
End Function


between End and End function, I want this function: so that if the limits is exceeded the numbers in the cell where the functions is used turns red to show that the result is not valid.


If eq_14_1>eq_14_2 or 1<<eq_14_2 Then
eq_14_2 or eq_15<(Sgavg*1000) Then
ActiveCell.Select
With Selection.Font
.Color = -16776961 'To change the text within the cell where the function is used to red text
.TintAndShade = 0
End With
End if



Thx a lot to anyone that can help </eq_14_2>

 

[wideboydixon]

Please use code tags when posting code; it makes it so much easier. Here's a stab at a solution.

Function Navier_Stokes_compressibleflow_Pressuredrop_ID(PipeID_mm, T1, P2, T2, Pipelength_m, PipeRoughness_mm, Flowrate_scf_min)
'March 2018, function to calculate pressure drop for a compressible flow system with pure N2 gas, by the use of the Navier stokes momentum equation
'INPUT:
'PipeID-mm in milimeter
'T1 and T2 in celsius
'P2 in Bar
'Pipelength_m in m
'PipeRoughness_mm in mm
'Flowrate_scf_min in standard cubic feet per minute
M = (Flowrate_scf_min / 35.3146667214886) / 60 'm3/s
m_kg_s = M * ((SRK_N2_Density(1.01325, T1) + SRK_N2_Density(1.01325, T2)) * 0.5) ' kg/s
IDM = PipeID_mm / 1000 'PipeID_mm converted to m
Tavg = (T1 + T2) * 0.5
T1K = T1 + 273.15 ' T1 in Celius converted to Kelvin
T2K = T2 + 273.15 ' T2 in Celius converted to Kelvin
TavgK = (T1K + T2K) * 0.5 ' average temperature Kelvin
P2p = P2 * 10 ^ 5 'P2 in Bar converted to Pa
R = 8.3145 'J/mol*K
Mw = 0.02801348 'kg/mol
Pi = 3.14159265358979
g = 9.81 'm2/s2
Elevation = 300 ' m
A = (IDM ^ 2 * Pi) / 4 'Area m2
Relativeroughness = PipeRoughness_mm / PipeID_mm
P1new = 1 * 10 ^ 5 'initial guess
For i = 1 To 1000
    P1 = P1new
    Pavg = (P1new + P2p) * 0.5 ' Average pressure
    D_visc_Pas = (N2_Viscosity_Dynamic((P1 * 10 ^ -5), T1) + N2_Viscosity_Dynamic(P2, T2)) * 0.5
    Sgavg = (SRK_N2_SG((P1 * 10 ^ -5), T1) + SRK_N2_SG(P2, T2)) * 0.5
    velocity = (M / A) / ((Pavg / (101325)) / ((SRK_N2_Z((P1 * 10 ^ -5), T1) + SRK_N2_Z(P2, T2)) * 0.5)) 'average velocity in m/s for the system based on the average pressure and compressibility factor(Z)
    Fd = Haaland_Equation_Darcy_frictionfactor(PipeRoughness_mm, PipeID_mm, Sgavg, velocity, D_visc_Pas) 'Haaland Friction factor
    P1new = (((2 * R * TavgK * m_kg_s ^ 2) / (Mw * A ^ 2)) * (Log(P1 / P2p) + ((Fd * Pipelength_m) / (2 * IDM))) + P2p ^ 2) ^ 0.5 'to find new P1
    Navier_Stokes_compressibleflow_Pressuredrop_ID = (P1 / 10 ^ 5) - P2 'output in Bar
    eq_14_1 = (Abs(Sgavg * 1000 * g * Elevation) / ((Fd * Pipelength_m * m_kg_s ^ 2) / (2 * A ^ 2 * IDM * Sgavg * 1000)))
    eq_14_2 = (Sgavg * 1000 * g) / ((Fd * m_kg_s ^ 2) / (2 * A ^ 2 * IDM * Sgavg * 1000))
    eq_15 = (((Fd) / (2 * IDM)) * (m_kg_s / A)) ^ 0.5
    error_horisontal = -((Pi ^ 2 * IDM ^ 6 * Mw * (P1 ^ 2 - P2p ^ 2)) / ((Fd * Pipelength_m) ^ 2 * R * TavgK * (((2 * IDM) / (Fd * Pipelength_m)) + (1 / (Log(P1 / P2p))))))
    If Abs(P1new - P1) < 10 Then Exit For
Next i

If i > 1000 Then End

If eq_14_1 > eq_14_2 Or 1 < eq_14_2 Or eq_15 < (Sgavg * 1000) Then
    ActiveCell.Font.Color = -16776961
    ActiveCell.Font.TintAndShade = 0
End If

End Function

WBD

 

[Daniel89]

Thank you very much WBD

 

[Daniel89]

I have come a far way on my project thanks to the people in the forum, thanks.

I have another question:

I would like to state the variables at the top of the module so that all functions that need them in the rest of the module take values from the top. This way it is easy to change the purpose of the calculation tool to fit other gases and materials and so on

This is my code:

Option Explicit

Sub Declare_variables()
Dim Pi, g, Tc, Pc, R, M, Mw, accFact
Pi = 3.14159265358979
g = 9.81 'm/s2
Tc = 126.15 'K
Pc = 3397800 'Pa
R = 8.3145 'J/mol*K
M = 28.01348 'kg/kmol
Mw = 0.02801348 'kg/mol
accFact = 0.04 ' Accentric Factor for N2
End Sub

Function Kelvin(C)
'Adds 273.15 to the input Celsius value to convert to Kelvin temperature
Kelvin = C + 273.15
End Function
Function Pressure_Pa(P_Bar)
'Function to convert Bar to Pascal pressure
Pressure_Pa = P_Bar * 10 ^ 5
End Function

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
 TempK = Kelvin(Celsius) 'Kelvin
 PressurePa = Pressure_Pa(Bar) 'Pa
 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 For
    Next i
If i > 100 Then End
End Function

It isn't working, the SRK_N2_Z function is not taking the variables from the sub and I get the error:

Variable not defined!