Option Explicit
Option Base 1
Sub Coefficients()
Dim Rswrl As Double
Dim Lswrl As Double
Dim Ror As Double
Dim Alpha As Double
Dim rho_liq As Double
Dim visc As Double
Dim rho_air As Double
Dim surfT_liquid As Double
Dim Omega As Double
Dim Re As Double
Dim Ain As Double
Dim V1 As Double
Dim Ri As Double
Dim R As Double
Dim psi As Double
Dim psipsi As Double
Dim a As Double
Dim b As Double
Dim h As Double
Dim sum As Double
Dim First As Double
Dim Last As Double
Dim A1 As Double
Dim A2 As Double
Dim B1 As Double
Dim B2 As Double
Dim C2 As Double
Dim A3 As Double
Dim B3 As Double
Dim B33 As Double
Dim C3 As Double
Dim B4 As Double
Dim A4 As Double
Dim C(1 To 2, 1 To 2) As Double
Dim D(1 To 2, 1 To 1) As Variant
Dim C_inv As Variant
Dim F As Variant
Dim Zswrl As Double
Dim dz1 As Double
Dim K1 As Double
Dim K2 As Double
Dim K3 As Double
Dim K4 As Double
Dim L1 As Double
Dim L2 As Double
Dim L3 As Double
Dim L4 As Double
Dim F1 As Double
Dim F2 As Double
Dim Z1 As Double
Dim Z2 As Double
Dim Z3 As Double
Dim Z4 As Double
Dim Del1 As Variant
Dim Del2 As Double
Dim Del3 As Double
Dim Del4 As Double
Dim W1 As Variant
Dim W2 As Double
Dim W3 As Double
Dim W4 As Double
Dim Del1St As Variant
Dim W1St As Variant
Dim Pi As Double
Dim npts As Long
Dim k As Long
Dim Nswrl As Long
Dim IZ As Long
' Input Values:
' Geometry:
Rswrl = ActiveSheet.Range("H18").Value / 2
Lswrl = ActiveSheet.Range("H14").Value
Ror = ActiveSheet.Range("H20").Value / 2
Alpha = ActiveSheet.Range("H21").Value
Ri = ActiveSheet.Range("H17").Value
' Properties:
rho_liq = ActiveSheet.Range("H8").Value
visc = ActiveSheet.Range("H9").Value
rho_air = ActiveSheet.Range("H10").Value
surfT_liquid = ActiveSheet.Range("H11").Value
' Preliminaries
Pi = 3.14159265
Ain = ActiveSheet.Range("H23").Value
V1 = ActiveSheet.Range("H24").Value
' Swirl Chamber
' Coefficients (constants - variable values must be added.
' Coefficient A1
npts = 20
b = 1#
a = 0#
h = (b - a) / npts
' First and Last Integration Points
First = (2 * a - a ^ 2) - (2 * a - a ^ 2) ^ 2
Last = (2 * b - b ^ 2) - (2 * b - b ^ 2) ^ 2
sum = First
' Odd terms
For k = 1 To npts - 1 Step 2
psi = 2 * k * h - (k * h) * (k * h)
sum = sum + 4 * (psi - psi * psi)
Next k
' Even terms
For k = 2 To npts - 1 Step 2
psi = 2 * k * h - (k * h) * (k * h)
sum = sum + 2 * (psi - psi * psi)
Next k
' Add final value
sum = sum + Last
A1 = h / 3 * sum
ActiveSheet.Range("R5") = A1
' Coefficient A2
npts = 20
b = 1#
a = 0#
h = (b - a) / npts
' First and Last Integration Points
First = (2 * a - a ^ 2)
Last = (2 * b - b ^ 2)
sum = First
' Odd terms
For k = 1 To npts - 1 Step 2
psi = 2 * k * h - (k * h) * (k * h)
sum = sum + 4 * psi
Next k
' Even terms
For k = 2 To npts - 1 Step 2
psi = 2 * k * h - (k * h) * (k * h)
sum = sum + 2 * psi
Next k
' Add final value
sum = sum + Last
A2 = h / 3 * sum
ActiveSheet.Range("R6") = A2
' Coefficient B1 and C2
b = 1#
a = 0#
B1 = (2 - 2 * b) - (2 - 2 * a)
ActiveSheet.Range("R7") = B1
C2 = B1
ActiveSheet.Range("R9") = C2
' Coefficient B2
npts = 20
b = 1#
a = 0#
h = (b - a) / npts
' First and Last Integration Points
First = (2 * a - a ^ 2) ^ 2
Last = (2 * b - b ^ 2) ^ 2
sum = First
' Odd terms
For k = 1 To npts - 1 Step 2
psipsi = (2 * k * h - (k * h) * (k * h)) * (2 * k * h - (k * h) * (k * h))
sum = sum + 4 * psipsi
Next k
' Even terms
For k = 2 To npts - 1 Step 2
psipsi = (2 * k * h - (k * h) * (k * h)) * (2 * k * h - (k * h) * (k * h))
sum = sum + 2 * psipsi
Next k
' Add final value
sum = sum + Last
B2 = h / 3 * sum
ActiveSheet.Range("R8") = B2
' Zone 1
' Circulation Constant
Omega = V1 * (Rswrl - Ri)
' Reynolds Number
Re = Omega / visc
' Initial value
W1 = 0.001
Del1 = 0.001
Z1 = 0
' Swirl Chamber Height
Zswrl = ActiveSheet.Range("H18").Value
' Number of Steps in Swirl Chamber Height
Nswrl = 20
' Delta-Z
dz1 = Zswrl / Nswrl
For IZ = 1 To Nswrl
' Swirl Chamber Coefficient Matrix
C(1, 1) = Del1 * W1
C(1, 2) = Del1 ^ 2
C(2, 1) = ((A2 - B2) * Del1 * W1 ^ 2)
C(2, 2) = ((A2 - 2 * B2 + 1) * Del1 ^ 2 * W1)
D(1, 1) = B1 / A1
D(2, 1) = C2 * W1
' Vector F Initial Conditions
F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D)
ActiveSheet.Range("Y5") = F(1, 1)
ActiveSheet.Range("Y6") = F(2, 1)
F1 = F(1, 1)
F2 = F(2, 1)
K1 = dz1 * F1
L1 = dz1 * F2
Z1 = Z1 + dz1 / 2
Del1 = Del1 + K1 / 2
W1 = W1 + L1 / 2
C(1, 1) = Del1 * W1
C(1, 2) = Del1 ^ 2
C(2, 1) = ((A2 - B2) * Del1 * W1 ^ 2)
C(2, 2) = ((A2 - 2 * B2 + 1) * Del1 ^ 2 * W1)
D(1, 1) = B1 / A1
D(2, 1) = C2 * W1
F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D)
F1 = F(1, 1)
F2 = F(2, 1)
K2 = dz1 * F1
L2 = dz1 * F2
Z2 = Z1 + dz1 / 2
Del2 = Del1 + K2 / 2
W2 = W1 + L2 / 2
C(1, 1) = Del2 * W2
C(1, 2) = Del2 ^ 2
C(2, 1) = ((A2 - B2) * Del2 * W2 ^ 2)
C(2, 2) = ((A2 - 2 * B2 + 1) * Del2 ^ 2 * W2)
D(1, 1) = B1 / A1
D(2, 1) = C2 * W2
F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D)
F1 = F(1, 1)
F2 = F(2, 1)
K3 = dz1 * F1
L3 = dz1 * F2
Z3 = Z1 + dz1 / 2
Del3 = Del1 + K2 / 2
W3 = W1 + L2 / 2
C(1, 1) = Del3 * W3
C(1, 2) = Del3 ^ 2
C(2, 1) = ((A2 - B2) * Del3 * W3 ^ 2)
C(2, 2) = ((A2 - 2 * B2 + 1) * Del3 ^ 2 * W3)
D(1, 1) = B1 / A1
D(2, 1) = C2 * W3
F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D)
F1 = F(1, 1)
F2 = F(2, 1)
K4 = dz1 * F1
L4 = dz1 * F2
Z3 = Z1 + dz1 / 2
Del4 = Del1 + K3
W4 = W1 + L3
C(1, 1) = Del4 * W4
C(1, 2) = Del4 ^ 2
C(2, 1) = ((A2 - B2) * Del4 * W4 ^ 2)
C(2, 2) = ((A2 - 2 * B2 + 1) * Del4 ^ 2 * W4)
D(1, 1) = B1 / A1
D(2, 1) = C2 * W4
Del1 = Del1 + (K1 + 2 * K2 + 2 * K3 + K4) / 6
W1 = W1 + (L1 + 2 * L2 + 2 * L3 + L4) / 6
Del1St(IZ) = Del1 'Type Mismatch
W1St(IZ) = W1 'Type Mismatch
Z1St(IZ) = Z1 'Type Mismatch
' Adjust Z1
Z1 = Z1 + dz1
Next IZ
' Spin Chamber
' Coefficients
A3 = A1
B3 = A2
B33 = B2
C3 = B1
B4 = C2
A4 = A1
ActiveSheet.Range("R10") = A3
ActiveSheet.Range("R11") = B3
ActiveSheet.Range("R12") = B33
ActiveSheet.Range("R13") = C3
ActiveSheet.Range("R14") = B4
ActiveSheet.Range("R15") = A4
' Spray Orifice
Omega = ActiveSheet.Range("H25")
Re = ActiveSheet.Range("H26")
R = Ror
' Solve Momentum Equations
End Sub