From 91c40a2ff5e2c1142ad185a0a3230374b0b50191 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=C3=85smund=20V=C3=A5ge=20Fannemel?=
 <34712686+asmfstatoil@users.noreply.github.com>
Date: Fri, 22 Nov 2024 13:08:06 +0100
Subject: [PATCH] refact: getPhases -> getPhase

---
 .../fluidmechanics/flownode/FlowNode.java     |  14 +-
 .../heatmasstransfercalc/FluidBoundary.java   |  32 ++--
 .../NonEquilibriumFluidBoundary.java          |  12 +-
 .../KrishnaStandartFilmModel.java             |  52 +++---
 .../ReactiveFluidBoundary.java                |   8 +-
 .../enhancementfactor/EnhancementFactor.java  |   2 +-
 .../EnhancementFactorAlg.java                 |  28 +--
 .../EnhancementFactorNumeric.java             |   4 +-
 .../InterphasePipeFlow.java                   |  13 +-
 .../InterphaseDropletFlow.java                |  43 +++--
 .../InterphaseStratifiedFlow.java             |  43 +++--
 .../InterphasePackedBed.java                  |   5 +-
 .../InterphaseStirredCellFlow.java            |  33 ++--
 .../stirredcellsolver/StirredCellSolver.java  |  47 +++--
 .../TwoPhaseFixedStaggeredGridSolver.java     | 162 ++++++++----------
 .../neqsim/thermo/ThermodynamicModelTest.java |   2 +-
 .../flashops/TPflash.java                     |   9 +-
 17 files changed, 243 insertions(+), 266 deletions(-)

diff --git a/src/main/java/neqsim/fluidmechanics/flownode/FlowNode.java b/src/main/java/neqsim/fluidmechanics/flownode/FlowNode.java
index f511f0c25..efb1c1720 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/FlowNode.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/FlowNode.java
@@ -477,9 +477,9 @@ public double calcNusseltNumber(double prandtlNumber, int phase) {
   /** {@inheritDoc} */
   @Override
   public double getPrandtlNumber(int phaseNum) {
-    return getBulkSystem().getPhases()[phaseNum].getCp()
-        * getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getViscosity()
-        / getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getConductivity();
+    return getBulkSystem().getPhase(phaseNum).getCp()
+        * getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getViscosity()
+        / getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getConductivity();
   }
 
   /** {@inheritDoc} */
@@ -558,10 +558,10 @@ public void increaseMolarRate(double moles) {}
   /** {@inheritDoc} */
   @Override
   public double calcTotalHeatTransferCoefficient(int phaseNum) {
-    double prandtlNumber = getBulkSystem().getPhases()[phaseNum].getCp()
-        / getBulkSystem().getPhases()[phaseNum].getNumberOfMolesInPhase()
-        * getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getViscosity()
-        / getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getConductivity();
+    double prandtlNumber = getBulkSystem().getPhase(phaseNum).getCp()
+        / getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
+        * getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getViscosity()
+        / getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getConductivity();
     double temp = 1.0 / (1.0
         / interphaseTransportCoefficient.calcWallHeatTransferCoefficient(phaseNum, prandtlNumber,
             this)
diff --git a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/FluidBoundary.java b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/FluidBoundary.java
index 7dbf7ff5b..086d2d3c5 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/FluidBoundary.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/FluidBoundary.java
@@ -242,28 +242,28 @@ public ThermodynamicOperations getBulkSystemOpertions() {
    */
   public void calcFluxTypeCorrectionMatrix(int phaseNum, int k) {
     fluxTypeCorrectionMatrixV[phaseNum] =
-        new Matrix(bulkSystem.getPhases()[phaseNum].getNumberOfComponents(), 1);
+        new Matrix(bulkSystem.getPhase(phaseNum).getNumberOfComponents(), 1);
     fluxTypeCorrectionMatrix[phaseNum] =
-        new Matrix(bulkSystem.getPhases()[phaseNum].getNumberOfComponents() - 1, 1);
+        new Matrix(bulkSystem.getPhase(phaseNum).getNumberOfComponents() - 1, 1);
     double temp = 0;
 
     double sum = 0;
-    for (int i = 0; i < bulkSystem.getPhases()[phaseNum].getNumberOfComponents() - 1; i++) {
+    for (int i = 0; i < bulkSystem.getPhase(phaseNum).getNumberOfComponents() - 1; i++) {
       temp = (i == k) ? 1.0 : 0.0;
       fluxTypeCorrectionMatrixV[phaseNum].set(i, 0, temp);
       sum += fluxTypeCorrectionMatrixV[phaseNum].get(i, 0)
-          * bulkSystem.getPhases()[phaseNum].getComponent(i).getx();
+          * bulkSystem.getPhase(phaseNum).getComponent(i).getx();
     }
 
     sum += fluxTypeCorrectionMatrixV[phaseNum]
-        .get(bulkSystem.getPhases()[phaseNum].getNumberOfComponents() - 1, 0)
-        * bulkSystem.getPhases()[phaseNum]
-            .getComponents()[bulkSystem.getPhases()[phaseNum].getNumberOfComponents() - 1].getx();
+        .get(bulkSystem.getPhase(phaseNum).getNumberOfComponents() - 1, 0)
+        * bulkSystem.getPhase(phaseNum)
+            .getComponents()[bulkSystem.getPhase(phaseNum).getNumberOfComponents() - 1].getx();
 
-    for (int i = 0; i < bulkSystem.getPhases()[phaseNum].getNumberOfComponents() - 1; i++) {
+    for (int i = 0; i < bulkSystem.getPhase(phaseNum).getNumberOfComponents() - 1; i++) {
       fluxTypeCorrectionMatrix[phaseNum].set(i, 0,
           (fluxTypeCorrectionMatrixV[phaseNum].get(i, 0) - fluxTypeCorrectionMatrixV[phaseNum]
-              .get(bulkSystem.getPhases()[phaseNum].getNumberOfComponents() - 1, 0)) / sum);
+              .get(bulkSystem.getPhase(phaseNum).getNumberOfComponents() - 1, 0)) / sum);
     }
   }
 
@@ -276,16 +276,16 @@ public void calcFluxTypeCorrectionMatrix(int phaseNum, int k) {
    */
   public void calcNonIdealCorrections(int phaseNum) {
     nonIdealCorrections[phaseNum] =
-        new Matrix(bulkSystem.getPhases()[phaseNum].getNumberOfComponents() - 1,
-            bulkSystem.getPhases()[phaseNum].getNumberOfComponents() - 1);
+        new Matrix(bulkSystem.getPhase(phaseNum).getNumberOfComponents() - 1,
+            bulkSystem.getPhase(phaseNum).getNumberOfComponents() - 1);
     double temp = 0;
-    for (int i = 0; i < bulkSystem.getPhases()[phaseNum].getNumberOfComponents() - 1; i++) {
-      for (int j = 0; j < bulkSystem.getPhases()[phaseNum].getNumberOfComponents() - 1; j++) {
+    for (int i = 0; i < bulkSystem.getPhase(phaseNum).getNumberOfComponents() - 1; i++) {
+      for (int j = 0; j < bulkSystem.getPhase(phaseNum).getNumberOfComponents() - 1; j++) {
         temp = (i == j) ? 1.0 : 0.0;
         nonIdealCorrections[phaseNum].set(i, j,
-            temp + bulkSystem.getPhases()[phaseNum].getComponent(i).getx()
-                * bulkSystem.getPhases()[phaseNum].getComponent(i).getdfugdn(j)
-                * bulkSystem.getPhases()[phaseNum].getNumberOfMolesInPhase());
+            temp + bulkSystem.getPhase(phaseNum).getComponent(i).getx()
+                * bulkSystem.getPhase(phaseNum).getComponent(i).getdfugdn(j)
+                * bulkSystem.getPhase(phaseNum).getNumberOfMolesInPhase());
         // her må det fylles inn
       }
     }
diff --git a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/NonEquilibriumFluidBoundary.java b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/NonEquilibriumFluidBoundary.java
index 4457f384f..32a161670 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/NonEquilibriumFluidBoundary.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/NonEquilibriumFluidBoundary.java
@@ -441,10 +441,10 @@ public void calcMolFractionDifference() {
    * @param phaseNum a int
    */
   public void calcHeatTransferCoefficients(int phaseNum) {
-    prandtlNumber[phaseNum] = getBulkSystem().getPhases()[phaseNum].getCp()
-        / getBulkSystem().getPhases()[phaseNum].getNumberOfMolesInPhase()
-        * getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getViscosity()
-        / getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getConductivity();
+    prandtlNumber[phaseNum] = getBulkSystem().getPhase(phaseNum).getCp()
+        / getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
+        * getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getViscosity()
+        / getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getConductivity();
     // System.out.println("prandtlNumber " + prandtlNumber[phase] + " interface " +
     // flowNode.getInterphaseTransportCoefficient().calcInterphaseHeatTransferCoefficient(phase,
     // prandtlNumber[phase], flowNode));
@@ -465,8 +465,8 @@ public void calcHeatTransferCorrection(int phaseNum) {
       if (Math.abs(nFlux.get(i, 0)) < 1e-20) {
         temp += 0.0;
       } else {
-        temp += nFlux.get(i, 0) * getBulkSystem().getPhases()[phaseNum].getCp()
-            / getBulkSystem().getPhases()[phaseNum].getNumberOfMolesInPhase()
+        temp += nFlux.get(i, 0) * getBulkSystem().getPhase(phaseNum).getCp()
+            / getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
             / heatTransferCoefficient[phaseNum]; // bulkSystem.getPhases()[0].getComponent(i).getNumberOfMolesInPhase()*
       }
     }
diff --git a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/KrishnaStandartFilmModel.java b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/KrishnaStandartFilmModel.java
index 3f9a00dd5..442794db6 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/KrishnaStandartFilmModel.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/KrishnaStandartFilmModel.java
@@ -91,17 +91,17 @@ public KrishnaStandartFilmModel clone() {
    * @return a double
    */
   public double calcBinarySchmidtNumbers(int phaseNum) {
-    for (int i = 0; i < getBulkSystem().getPhases()[phaseNum].getNumberOfComponents(); i++) {
-      for (int j = 0; j < getBulkSystem().getPhases()[phaseNum].getNumberOfComponents(); j++) {
+    for (int i = 0; i < getBulkSystem().getPhase(phaseNum).getNumberOfComponents(); i++) {
+      for (int j = 0; j < getBulkSystem().getPhase(phaseNum).getNumberOfComponents(); j++) {
         binarySchmidtNumber[phaseNum][i][j] =
-            getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity()
-                / getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
+            getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
+                / getBulkSystem().getPhase(phaseNum).getPhysicalProperties()
                     .getDiffusionCoefficient(i, j);
         // System.out.println("i j " + i +" j " + j);
         // System.out.println("phase " + phase + " diff" +
-        // getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getDiffusionCoefficient(i,j));
+        // getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDiffusionCoefficient(i,j));
         // System.out.println("phase " + phase + " visk" +
-        // getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity());
+        // getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity());
       }
     }
     return 1;
@@ -116,8 +116,8 @@ public double calcBinarySchmidtNumbers(int phaseNum) {
    * @return a double
    */
   public double calcBinaryMassTransferCoefficients(int phaseNum) {
-    for (int i = 0; i < getBulkSystem().getPhases()[phaseNum].getNumberOfComponents(); i++) {
-      for (int j = 0; j < getBulkSystem().getPhases()[phaseNum].getNumberOfComponents(); j++) {
+    for (int i = 0; i < getBulkSystem().getPhase(phaseNum).getNumberOfComponents(); i++) {
+      for (int j = 0; j < getBulkSystem().getPhase(phaseNum).getNumberOfComponents(); j++) {
         binaryMassTransferCoefficient[phaseNum][i][j] =
             flowNode.getInterphaseTransportCoefficient().calcInterphaseMassTransferCoefficient(
                 phaseNum, binarySchmidtNumber[phaseNum][i][j], flowNode);
@@ -135,24 +135,24 @@ public double calcBinaryMassTransferCoefficients(int phaseNum) {
    * @return a double
    */
   public double calcMassTransferCoefficients(int phaseNum) {
-    int n = getBulkSystem().getPhases()[phaseNum].getNumberOfComponents() - 1;
+    int n = getBulkSystem().getPhase(phaseNum).getNumberOfComponents() - 1;
 
-    for (int i = 0; i < getBulkSystem().getPhases()[phaseNum].getNumberOfComponents() - 1; i++) {
+    for (int i = 0; i < getBulkSystem().getPhase(phaseNum).getNumberOfComponents() - 1; i++) {
       double tempVar = 0;
-      for (int j = 0; j < getBulkSystem().getPhases()[phaseNum].getNumberOfComponents(); j++) {
+      for (int j = 0; j < getBulkSystem().getPhase(phaseNum).getNumberOfComponents(); j++) {
         if (i != j) {
-          tempVar += getBulkSystem().getPhases()[phaseNum].getComponent(j).getx()
+          tempVar += getBulkSystem().getPhase(phaseNum).getComponent(j).getx()
               / binaryMassTransferCoefficient[phaseNum][i][j];
         }
         if (j < n) {
           massTransferCoefficientMatrix[phaseNum].set(i, j,
-              -getBulkSystem().getPhases()[phaseNum].getComponent(i).getx()
+              -getBulkSystem().getPhase(phaseNum).getComponent(i).getx()
                   * (1.0 / binaryMassTransferCoefficient[phaseNum][i][j]
                       - 1.0 / binaryMassTransferCoefficient[phaseNum][i][n]));
         }
       }
       massTransferCoefficientMatrix[phaseNum].set(i, i,
-          tempVar + getBulkSystem().getPhases()[phaseNum].getComponent(i).getx()
+          tempVar + getBulkSystem().getPhase(phaseNum).getComponent(i).getx()
               / binaryMassTransferCoefficient[phaseNum][i][n]);
     }
     massTransferCoefficientMatrix[phaseNum] = massTransferCoefficientMatrix[phaseNum].inverse();
@@ -167,28 +167,28 @@ public double calcMassTransferCoefficients(int phaseNum) {
    * @param phaseNum a int
    */
   public void calcPhiMatrix(int phaseNum) {
-    int n = getBulkSystem().getPhases()[phaseNum].getNumberOfComponents() - 1;
+    int n = getBulkSystem().getPhase(phaseNum).getNumberOfComponents() - 1;
 
-    for (int i = 0; i < getBulkSystem().getPhases()[phaseNum].getNumberOfComponents() - 1; i++) {
+    for (int i = 0; i < getBulkSystem().getPhase(phaseNum).getNumberOfComponents() - 1; i++) {
       double tempVar = 0;
-      for (int j = 0; j < getBulkSystem().getPhases()[phaseNum].getNumberOfComponents(); j++) {
+      for (int j = 0; j < getBulkSystem().getPhase(phaseNum).getNumberOfComponents(); j++) {
         if (i != j || i == n) {
-          tempVar += nFlux.get(i, 0)
-              / (1.0 / (getBulkSystem().getPhases()[phaseNum].getMolarVolume() * 1e-5)
+          tempVar +=
+              nFlux.get(i, 0) / (1.0 / (getBulkSystem().getPhase(phaseNum).getMolarVolume() * 1e-5)
                   * binaryMassTransferCoefficient[phaseNum][i][j]);
         }
         if (j < n) {
           phiMatrix.set(i, j,
               -nFlux.get(i, 0) * (1.0
-                  / (1.0 / (getBulkSystem().getPhases()[phaseNum].getMolarVolume() * 1e-5)
+                  / (1.0 / (getBulkSystem().getPhase(phaseNum).getMolarVolume() * 1e-5)
                       * binaryMassTransferCoefficient[phaseNum][i][j])
-                  - 1.0 / (1.0 / (getBulkSystem().getPhases()[phaseNum].getMolarVolume() * 1e-5)
+                  - 1.0 / (1.0 / (getBulkSystem().getPhase(phaseNum).getMolarVolume() * 1e-5)
                       * binaryMassTransferCoefficient[phaseNum][i][n])));
         }
       }
       phiMatrix.set(i, i,
           tempVar + nFlux.get(i, 0)
-              / (1.0 / (getBulkSystem().getPhases()[phaseNum].getMolarVolume() * 1e-5)
+              / (1.0 / (getBulkSystem().getPhase(phaseNum).getMolarVolume() * 1e-5)
                   * binaryMassTransferCoefficient[phaseNum][i][n]));
     }
   }
@@ -212,7 +212,7 @@ public void calcRedPhiMatrix(int phase) {
    * @param phaseNum a int
    */
   public void calcRedCorrectionMatrix(int phaseNum) {
-    for (int i = 0; i < getBulkSystem().getPhases()[phaseNum].getNumberOfComponents() - 1; i++) {
+    for (int i = 0; i < getBulkSystem().getPhase(phaseNum).getNumberOfComponents() - 1; i++) {
       redCorrectionMatrix.set(i, 0, (redPhiMatrix.get(0, i) * Math.exp(redPhiMatrix.get(0, i)))
           / (Math.exp(redPhiMatrix.get(0, i)) - (1.0 - 1e-15)));
     }
@@ -228,9 +228,9 @@ public void calcRedCorrectionMatrix(int phaseNum) {
   public void calcCorrectionMatrix(int phaseNum) {
     Matrix modalPhiMatrix = phiMatrix.eig().getV();
     Matrix diagonalRedCorrectionMatrix =
-        new Matrix(getBulkSystem().getPhases()[phaseNum].getNumberOfComponents() - 1,
-            getBulkSystem().getPhases()[phaseNum].getNumberOfComponents() - 1);
-    for (int i = 0; i < getBulkSystem().getPhases()[phaseNum].getNumberOfComponents() - 1; i++) {
+        new Matrix(getBulkSystem().getPhase(phaseNum).getNumberOfComponents() - 1,
+            getBulkSystem().getPhase(phaseNum).getNumberOfComponents() - 1);
+    for (int i = 0; i < getBulkSystem().getPhase(phaseNum).getNumberOfComponents() - 1; i++) {
       diagonalRedCorrectionMatrix.set(i, i, redCorrectionMatrix.get(i, 0));
     }
     rateCorrectionMatrix[phaseNum] =
diff --git a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/ReactiveFluidBoundary.java b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/ReactiveFluidBoundary.java
index 8586a5d2c..6738f5a36 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/ReactiveFluidBoundary.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/ReactiveFluidBoundary.java
@@ -410,9 +410,9 @@ public void calcMolFractionDifference() {
   /** {@inheritDoc} */
   @Override
   public void calcHeatTransferCoefficients(int phaseNum) {
-    prandtlNumber[phaseNum] = getBulkSystem().getPhases()[phaseNum].getCp()
-        * getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getViscosity()
-        / getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getConductivity();
+    prandtlNumber[phaseNum] = getBulkSystem().getPhase(phaseNum).getCp()
+        * getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getViscosity()
+        / getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getConductivity();
     heatTransferCoefficient[phaseNum] = flowNode.getInterphaseTransportCoefficient()
         .calcInterphaseHeatTransferCoefficient(phaseNum, prandtlNumber[phaseNum], flowNode);
   }
@@ -423,7 +423,7 @@ public void calcHeatTransferCorrection(int phaseNum) {
     double temp = 0;
     for (int i = 0; i < bulkSystem.getPhases()[0].getNumberOfMolecularComponents(); i++) {
       temp += bulkSystem.getPhases()[0].getComponent(i).getNumberOfMolesInPhase()
-          * getBulkSystem().getPhases()[phaseNum].getCp() / heatTransferCoefficient[phaseNum];
+          * getBulkSystem().getPhase(phaseNum).getCp() / heatTransferCoefficient[phaseNum];
     }
     heatTransferCorrection[phaseNum] = temp;
   }
diff --git a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/enhancementfactor/EnhancementFactor.java b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/enhancementfactor/EnhancementFactor.java
index a2f0431e8..29a31c257 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/enhancementfactor/EnhancementFactor.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/enhancementfactor/EnhancementFactor.java
@@ -110,7 +110,7 @@ public void setEnhancementVec(double[] enhancementVec) {
    * @param phaseNum a int
    */
   public void setOnesVec(int phaseNum) {
-    for (int j = 0; j < fluidBoundary.getBulkSystem().getPhases()[phaseNum]
+    for (int j = 0; j < fluidBoundary.getBulkSystem().getPhase(phaseNum)
         .getNumberOfComponents(); j++) {
       enhancementVec[j] = 1.0;
     }
diff --git a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/enhancementfactor/EnhancementFactorAlg.java b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/enhancementfactor/EnhancementFactorAlg.java
index 32e2ca141..c4e537a64 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/enhancementfactor/EnhancementFactorAlg.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/enhancementfactor/EnhancementFactorAlg.java
@@ -29,15 +29,15 @@ public EnhancementFactorAlg(FluidBoundaryInterface fluidBoundary) {
   @Override
   public void calcEnhancementVec(int phaseNum) {
     double hatta = 0.0;
-    for (int j = 0; j < fluidBoundary.getBulkSystem().getPhases()[phaseNum]
+    for (int j = 0; j < fluidBoundary.getBulkSystem().getPhase(phaseNum)
         .getNumberOfComponents(); j++) {
-      if (fluidBoundary.getBulkSystem().getPhases()[phaseNum].getComponent(j).getName()
-          .equals("CO2") && phaseNum == 1) {
+      if (fluidBoundary.getBulkSystem().getPhase(phaseNum).getComponent(j).getName().equals("CO2")
+          && phaseNum == 1) {
         enhancementVec[j] = fluidBoundary.getBulkSystem().getChemicalReactionOperations()
             .solveKinetics(phaseNum, fluidBoundary.getInterphaseSystem().getPhase(phaseNum), j);
         // System.out.println("enh " + enhancementVec[j]);
         hatta = Math
-            .sqrt(enhancementVec[j] * fluidBoundary.getBulkSystem().getPhases()[phaseNum]
+            .sqrt(enhancementVec[j] * fluidBoundary.getBulkSystem().getPhase(phaseNum)
                 .getPhysicalProperties().getEffectiveDiffusionCoefficient(j))
             / fluidBoundary.getEffectiveMassTransferCoefficient(phaseNum, j);
         hattaNumber[j] = hatta;
@@ -65,28 +65,28 @@ public void calcEnhancementVec(int phaseNum) {
   // fluidBoundary.getBulkSystem().getChemicalReactionOperations().getReactionList().getReacMatrix();
   // double[][] stocMatrix =
   // fluidBoundary.getBulkSystem().getChemicalReactionOperations().getReactionList().getStocMatrix();
-  // fluidBoundary.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().calcEffectiveDiffusionCoefficients();
+  // fluidBoundary.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().calcEffectiveDiffusionCoefficients();
 
   // for(int
-  // i=0;i<fluidBoundary.getBulkSystem().getPhases()[phaseNum].getNumberOfComponents();i++){
+  // i=0;i<fluidBoundary.getBulkSystem().getPhase(phaseNum).getNumberOfComponents();i++){
   // for(int
-  // j=0;j<fluidBoundary.getBulkSystem().getPhases()[phaseNum].getNumberOfComponents();j++){
+  // j=0;j<fluidBoundary.getBulkSystem().getPhase(phaseNum).getNumberOfComponents();j++){
   // if(i==j || (stocMatrix[i][i]*stocMatrix[i][j])<=0)
   // enhancementFactor[phase].set(i,j, 1.0);
   // else{
   // double pseudoReacRate =
-  // reacMatrix[i][j]*fluidBoundary.getBulkSystem().getPhases()[phaseNum].getComponent(j).getx()*fluidBoundary.getBulkSystem().getPhases()[phaseNum].getDensity()/fluidBoundary.getBulkSystem().getPhases()[phaseNum].getMolarMass()*Math.abs(stocMatrix[i][i]);
-  // hattaNumber[phase].set(i,j,Math.sqrt(pseudoReacRate*fluidBoundary.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getEffectiveDiffusionCoefficient(i))/fluidBoundary.getBinaryMassTransferCoefficient(phase,i,j));
+  // reacMatrix[i][j]*fluidBoundary.getBulkSystem().getPhase(phaseNum).getComponent(j).getx()*fluidBoundary.getBulkSystem().getPhase(phaseNum).getDensity()/fluidBoundary.getBulkSystem().getPhase(phaseNum).getMolarMass()*Math.abs(stocMatrix[i][i]);
+  // hattaNumber[phase].set(i,j,Math.sqrt(pseudoReacRate*fluidBoundary.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getEffectiveDiffusionCoefficient(i))/fluidBoundary.getBinaryMassTransferCoefficient(phase,i,j));
   // double phiVal = 1.0 +
-  // (fluidBoundary.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getEffectiveDiffusionCoefficient(j)/fluidBoundary.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getEffectiveDiffusionCoefficient(i))*
+  // (fluidBoundary.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getEffectiveDiffusionCoefficient(j)/fluidBoundary.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getEffectiveDiffusionCoefficient(i))*
   // (stocMatrix[i][i]/stocMatrix[i][j]) *
-  // (fluidBoundary.getBulkSystem().getPhases()[phaseNum].getComponent(j).getx()*fluidBoundary.getBulkSystem().getPhases()[phaseNum].getDensity()/fluidBoundary.getBulkSystem().getPhases()[phaseNum].getMolarMass())/
-  // (fluidBoundary.getInterphaseSystem().getPhases()[phaseNum].getComponent(i).getx()*fluidBoundary.getInterphaseSystem().getPhases()[phaseNum].getDensity()/fluidBoundary.getInterphaseSystem().getPhases()[phaseNum].getMolarMass());
+  // (fluidBoundary.getBulkSystem().getPhase(phaseNum).getComponent(j).getx()*fluidBoundary.getBulkSystem().getPhase(phaseNum).getDensity()/fluidBoundary.getBulkSystem().getPhase(phaseNum).getMolarMass())/
+  // (fluidBoundary.getInterphaseSystem().getPhase(phaseNum).getComponent(i).getx()*fluidBoundary.getInterphaseSystem().getPhase(phaseNum).getDensity()/fluidBoundary.getInterphaseSystem().getPhase(phaseNum).getMolarMass());
 
   // System.out.println("components " +
-  // fluidBoundary.getBulkSystem().getPhases()[phaseNum].getComponent(i).getComponentName()
+  // fluidBoundary.getBulkSystem().getPhase(phaseNum).getComponent(i).getComponentName()
   // + " " +
-  // fluidBoundary.getBulkSystem().getPhases()[phaseNum].getComponent(j).getComponentName());
+  // fluidBoundary.getBulkSystem().getPhase(phaseNum).getComponent(j).getComponentName());
   // System.out.println("hatta : " + hattaNumber[phase].get(i,j));
   // System.out.println("phi : " + phiVal);
 
diff --git a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/enhancementfactor/EnhancementFactorNumeric.java b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/enhancementfactor/EnhancementFactorNumeric.java
index e01375eac..5cb09b2b0 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/enhancementfactor/EnhancementFactorNumeric.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/heatmasstransfercalc/nonequilibriumfluidboundary/filmmodelboundary/reactivefilmmodel/enhancementfactor/EnhancementFactorNumeric.java
@@ -45,9 +45,9 @@ public void calcEnhancementMatrix(int phaseNum) {
     nonReactiveInterface.createSystem();
     nonReactiveInterface.solve();
     reactiveInterface.solve();
-    for (int i = 0; i < fluidBoundary.getBulkSystem().getPhases()[phaseNum]
+    for (int i = 0; i < fluidBoundary.getBulkSystem().getPhase(phaseNum)
         .getNumberOfComponents(); i++) {
-      for (int j = 0; j < fluidBoundary.getBulkSystem().getPhases()[phaseNum]
+      for (int j = 0; j < fluidBoundary.getBulkSystem().getPhase(phaseNum)
           .getNumberOfComponents(); j++) {
         // enhancementFactor[1].set(i,j,0);
         // System.out.println("num enhancement " + enhancementFactor[1].get(i,j));
diff --git a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphaseonephase/interphasepipeflow/InterphasePipeFlow.java b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphaseonephase/interphasepipeflow/InterphasePipeFlow.java
index bebacb28a..7196d69a0 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphaseonephase/interphasepipeflow/InterphasePipeFlow.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphaseonephase/interphasepipeflow/InterphasePipeFlow.java
@@ -53,14 +53,14 @@ public double calcWallHeatTransferCoefficient(int phaseNum, double prandtlNumber
       FlowNodeInterface node) {
     if (Math.abs(node.getReynoldsNumber()) < 2000) {
       return 3.66 / node.getGeometry().getDiameter()
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getConductivity();
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getConductivity();
     }
     // if turbulent - use chilton colburn analogy
     else {
-      double temp = node.getBulkSystem().getPhases()[phaseNum].getCp()
-          / node.getBulkSystem().getPhases()[phaseNum].getMolarMass()
-          / node.getBulkSystem().getPhases()[phaseNum].getNumberOfMolesInPhase()
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getDensity()
+      double temp = node.getBulkSystem().getPhase(phaseNum).getCp()
+          / node.getBulkSystem().getPhase(phaseNum).getMolarMass()
+          / node.getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity()
           * node.getVelocity();
       return 0.5 * this.calcWallFrictionFactor(phaseNum, node) * Math.pow(prandtlNumber, -2.0 / 3.0)
           * temp;
@@ -73,8 +73,7 @@ public double calcWallMassTransferCoefficient(int phaseNum, double schmidtNumber
       FlowNodeInterface node) {
     if (Math.abs(node.getReynoldsNumber()) < 2000) {
       return 3.66 / node.getGeometry().getDiameter() / schmidtNumber
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-              .getKinematicViscosity();
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity();
     } else {
       double temp = node.getVelocity();
       return 0.5 * this.calcWallFrictionFactor(phaseNum, node) * Math.pow(schmidtNumber, -2.0 / 3.0)
diff --git a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/interphasepipeflow/InterphaseDropletFlow.java b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/interphasepipeflow/InterphaseDropletFlow.java
index ae46e2d00..43d8fffab 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/interphasepipeflow/InterphaseDropletFlow.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/interphasepipeflow/InterphaseDropletFlow.java
@@ -56,13 +56,13 @@ public double calcWallHeatTransferCoefficient(int phaseNum, double prandtlNumber
       FlowNodeInterface node) {
     if (Math.abs(node.getReynoldsNumber(phaseNum)) < 2000) {
       return 3.66 / node.getHydraulicDiameter(phaseNum)
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getConductivity();
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getConductivity();
     } else {
       // if turbulent - use chilton colburn analogy
-      double temp = node.getBulkSystem().getPhases()[phaseNum].getCp()
-          / node.getBulkSystem().getPhases()[phaseNum].getMolarMass()
-          / node.getBulkSystem().getPhases()[phaseNum].getNumberOfMolesInPhase()
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getDensity()
+      double temp = node.getBulkSystem().getPhase(phaseNum).getCp()
+          / node.getBulkSystem().getPhase(phaseNum).getMolarMass()
+          / node.getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity()
           * node.getVelocity(phaseNum);
       return 0.5 * this.calcWallFrictionFactor(phaseNum, node) * Math.pow(prandtlNumber, -2.0 / 3.0)
           * temp;
@@ -76,13 +76,13 @@ public double calcInterphaseHeatTransferCoefficient(int phaseNum, double prandtl
     // System.out.println("velocity " + node.getVelocity(phase));
     if (Math.abs(node.getReynoldsNumber()) < 2000) {
       return 3.66 / node.getHydraulicDiameter(phaseNum)
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getConductivity();
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getConductivity();
     } else {
       // if turbulent - use chilton colburn analogy
-      double temp = node.getBulkSystem().getPhases()[phaseNum].getCp()
-          / node.getBulkSystem().getPhases()[phaseNum].getMolarMass()
-          / node.getBulkSystem().getPhases()[phaseNum].getNumberOfMolesInPhase()
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getDensity()
+      double temp = node.getBulkSystem().getPhase(phaseNum).getCp()
+          / node.getBulkSystem().getPhase(phaseNum).getMolarMass()
+          / node.getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity()
           * node.getVelocity(phaseNum);
 
       return 0.5 * this.calcWallFrictionFactor(phaseNum, node) * Math.pow(prandtlNumber, -2.0 / 3.0)
@@ -96,8 +96,7 @@ public double calcWallMassTransferCoefficient(int phaseNum, double schmidtNumber
       FlowNodeInterface node) {
     if (Math.abs(node.getReynoldsNumber()) < 2000) {
       return 3.66 / node.getHydraulicDiameter(phaseNum) / schmidtNumber
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-              .getKinematicViscosity();
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity();
     } else {
       double temp = node.getVelocity(phaseNum);
       return 0.5 * this.calcWallFrictionFactor(phaseNum, node) * Math.pow(schmidtNumber, -2.0 / 3.0)
@@ -128,10 +127,9 @@ public double calcInterphaseMassTransferCoefficient(int phaseNum, double schmidt
       // er usikker paa denne korreksjonen med 1e-2 - maa sjekkes opp mot artikkel av
       // Yih og Chen (1982) - satser paa at de ga den med enhet cm/sek
       massTrans = redMassTrans
-          * Math.pow(Math.pow(node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
+          * Math.pow(Math.pow(node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties()
               .getKinematicViscosity(), 2.0) / gravity, -1.0 / 3.0)
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-              .getKinematicViscosity()
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
           / schmidtNumber;
     }
     if (phaseNum == 0) {
@@ -140,15 +138,14 @@ public double calcInterphaseMassTransferCoefficient(int phaseNum, double schmidt
         // System.out.println("hyd diam " + node.getHydraulicDiameter(phaseNum) +" phaseNum "
         // + phaseNum);
         // System.out.println("kin visk " +
-        // node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity()
+        // node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
         // +" phaseNum " + phaseNum);
         // System.out.println("mas " +3.66 / node.getHydraulicDiameter(phaseNum) /
         // schmidtNumber *
-        // node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity()
+        // node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
         // +" phaseNum " + phaseNum);
-        massTrans = 3.66 / node.getHydraulicDiameter(phaseNum) / schmidtNumber
-            * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-                .getKinematicViscosity();
+        massTrans = 3.66 / node.getHydraulicDiameter(phaseNum) / schmidtNumber * node
+            .getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity();
       } else {
         double temp = node.getVelocity(phaseNum);
         // System.out.println("mass " + 1e-5* 0.5 * this.calcWallFrictionFactor(phaseNum,
@@ -172,16 +169,16 @@ public double calcInterphaseMassTransferCoefficient(int phaseNum, double schmidt
   // //System.out.println("red " + redMassTrans/Math.pow(node.getReynoldsNumber(phaseNum),c3));
   // //System.out.println("sc " + schmidtNumber);
   // massTrans =
-  // redMassTrans*node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity()
+  // redMassTrans*node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
   // / schmidtNumber / node.getGeometry().getDiameter();
   // }
   // if(phaseNum==0){
   // //System.out.println("diff " +
-  // node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity() /
+  // node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity() /
   // schmidtNumber);
 
   // //massTrans = 3.66 / node.getHydraulicDiameter(phaseNum) / schmidtNumber *
-  // node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity();
+  // node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity();
   // massTrans=0.010;
   // }
   // //System.out.println("mass trans " +massTrans + " phaseNum " + phaseNum);
diff --git a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/interphasepipeflow/InterphaseStratifiedFlow.java b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/interphasepipeflow/InterphaseStratifiedFlow.java
index 07f19c9f5..5374e811c 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/interphasepipeflow/InterphaseStratifiedFlow.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/interphasepipeflow/InterphaseStratifiedFlow.java
@@ -56,14 +56,14 @@ public double calcWallHeatTransferCoefficient(int phaseNum, double prandtlNumber
       FlowNodeInterface node) {
     if (Math.abs(node.getReynoldsNumber(phaseNum)) < 2000) {
       return 3.66 / node.getHydraulicDiameter(phaseNum)
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getConductivity();
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getConductivity();
     } else {
       // if turbulent - use chilton colburn analogy
 
-      double temp = node.getBulkSystem().getPhases()[phaseNum].getCp()
-          / node.getBulkSystem().getPhases()[phaseNum].getMolarMass()
-          / node.getBulkSystem().getPhases()[phaseNum].getNumberOfMolesInPhase()
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getDensity()
+      double temp = node.getBulkSystem().getPhase(phaseNum).getCp()
+          / node.getBulkSystem().getPhase(phaseNum).getMolarMass()
+          / node.getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity()
           * node.getVelocity(phaseNum);
       return 0.5 * this.calcWallFrictionFactor(phaseNum, node) * Math.pow(prandtlNumber, -2.0 / 3.0)
           * temp;
@@ -77,13 +77,13 @@ public double calcInterphaseHeatTransferCoefficient(int phaseNum, double prandtl
     // System.out.println("velocity " + node.getVelocity(phase));
     if (Math.abs(node.getReynoldsNumber()) < 2000) {
       return 3.66 / node.getHydraulicDiameter(phaseNum)
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getConductivity();
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getConductivity();
     } else {
       // if turbulent - use chilton colburn analogy
-      double temp = node.getBulkSystem().getPhases()[phaseNum].getCp()
-          / node.getBulkSystem().getPhases()[phaseNum].getMolarMass()
-          / node.getBulkSystem().getPhases()[phaseNum].getNumberOfMolesInPhase()
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getDensity()
+      double temp = node.getBulkSystem().getPhase(phaseNum).getCp()
+          / node.getBulkSystem().getPhase(phaseNum).getMolarMass()
+          / node.getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity()
           * node.getVelocity(phaseNum);
 
       return 0.5 * this.calcWallFrictionFactor(phaseNum, node) * Math.pow(prandtlNumber, -2.0 / 3.0)
@@ -97,8 +97,7 @@ public double calcWallMassTransferCoefficient(int phaseNum, double schmidtNumber
       FlowNodeInterface node) {
     if (Math.abs(node.getReynoldsNumber()) < 2000) {
       return 3.66 / node.getHydraulicDiameter(phaseNum) / schmidtNumber
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-              .getKinematicViscosity();
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity();
     } else {
       double temp = node.getVelocity(phaseNum);
       return 0.5 * this.calcWallFrictionFactor(phaseNum, node) * Math.pow(schmidtNumber, -2.0 / 3.0)
@@ -129,10 +128,9 @@ public double calcInterphaseMassTransferCoefficient(int phaseNum, double schmidt
       // er usikker paa denne korreksjonen med 1e-2 - maa sjekkes opp mot artikkel av
       // Yih og Chen (1982) - satser paa at de ga den med enhet cm/sek
       massTrans = redMassTrans
-          * Math.pow(Math.pow(node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
+          * Math.pow(Math.pow(node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties()
               .getKinematicViscosity(), 2.0) / gravity, -1.0 / 3.0)
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-              .getKinematicViscosity()
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
           / schmidtNumber;
     }
     if (phaseNum == 0) {
@@ -141,15 +139,14 @@ public double calcInterphaseMassTransferCoefficient(int phaseNum, double schmidt
         // System.out.println("hyd diam " + node.getHydraulicDiameter(phaseNum) +" phaseNum "
         // + phaseNum);
         // System.out.println("kin visk " +
-        // node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity()
+        // node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
         // +" phaseNum " + phaseNum);
         // System.out.println("mas " +3.66 / node.getHydraulicDiameter(phaseNum) /
         // schmidtNumber *
-        // node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity()
+        // node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
         // +" phaseNum " + phaseNum);
-        massTrans = 3.66 / node.getHydraulicDiameter(phaseNum) / schmidtNumber
-            * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-                .getKinematicViscosity();
+        massTrans = 3.66 / node.getHydraulicDiameter(phaseNum) / schmidtNumber * node
+            .getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity();
       } else {
         double temp = node.getVelocity(phaseNum);
         // System.out.println("mass " + 1e-5* 0.5 * this.calcWallFrictionFactor(phaseNum,
@@ -173,16 +170,16 @@ public double calcInterphaseMassTransferCoefficient(int phaseNum, double schmidt
   // //System.out.println("red " + redMassTrans/Math.pow(node.getReynoldsNumber(phaseNum),c3));
   // //System.out.println("sc " + schmidtNumber);
   // massTrans =
-  // redMassTrans*node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity()
+  // redMassTrans*node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
   // / schmidtNumber / node.getGeometry().getDiameter();
   // }
   // if(phaseNum==0){
   // //System.out.println("diff " +
-  // node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity() /
+  // node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity() /
   // schmidtNumber);
 
   // //massTrans = 3.66 / node.getHydraulicDiameter(phaseNum) / schmidtNumber *
-  // node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity();
+  // node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity();
   // massTrans=0.010;
   // }
   // //System.out.println("mass trans " +massTrans + " phaseNum " + phaseNum);
diff --git a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/interphasereactorflow/InterphasePackedBed.java b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/interphasereactorflow/InterphasePackedBed.java
index c1cf811c2..53cd71b8c 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/interphasereactorflow/InterphasePackedBed.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/interphasereactorflow/InterphasePackedBed.java
@@ -82,7 +82,7 @@ public double calcInterphaseMassTransferCoefficient(int phaseNum, double schmidt
               * Math.pow(node.getGeometry().getPacking().getSurfaceAreaPrVolume()
                   * node.getGeometry().getPacking().getSize(), 0.4);
       massTrans = redMassTrans * Math.pow(
-          node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getKinematicViscosity()
+          node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
               * gravity,
           1.0 / 3.0);
       System.out.println("mas trans liq " + massTrans);
@@ -93,8 +93,7 @@ public double calcInterphaseMassTransferCoefficient(int phaseNum, double schmidt
               * Math.pow(node.getGeometry().getPacking().getSurfaceAreaPrVolume()
                   * node.getGeometry().getPacking().getSize(), -2.0);
       massTrans = redMassTrans * node.getGeometry().getPacking().getSurfaceAreaPrVolume()
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-              .getKinematicViscosity()
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
           / schmidtNumber;
       System.out.println("mas trans gas " + massTrans);
     }
diff --git a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/stirredcell/InterphaseStirredCellFlow.java b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/stirredcell/InterphaseStirredCellFlow.java
index 4e958af34..925e31606 100644
--- a/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/stirredcell/InterphaseStirredCellFlow.java
+++ b/src/main/java/neqsim/fluidmechanics/flownode/fluidboundary/interphasetransportcoefficient/interphasetwophase/stirredcell/InterphaseStirredCellFlow.java
@@ -38,13 +38,13 @@ public double calcWallHeatTransferCoefficient(int phaseNum, double prandtlNumber
       FlowNodeInterface node) {
     if (Math.abs(node.getReynoldsNumber(phaseNum)) < 2000) {
       return 3.66 / node.getHydraulicDiameter(phaseNum)
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getConductivity();
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getConductivity();
     } else {
       // if turbulent - use chilton colburn analogy
-      double temp = node.getBulkSystem().getPhases()[phaseNum].getCp()
-          / node.getBulkSystem().getPhases()[phaseNum].getMolarMass()
-          / node.getBulkSystem().getPhases()[phaseNum].getNumberOfMolesInPhase()
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getDensity()
+      double temp = node.getBulkSystem().getPhase(phaseNum).getCp()
+          / node.getBulkSystem().getPhase(phaseNum).getMolarMass()
+          / node.getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity()
           * node.getVelocity(phaseNum);
       return 0.5 * this.calcWallFrictionFactor(phaseNum, node) * Math.pow(prandtlNumber, -2.0 / 3.0)
           * temp;
@@ -57,14 +57,14 @@ public double calcInterphaseHeatTransferCoefficient(int phaseNum, double prandtl
       FlowNodeInterface node) {
     if (Math.abs(node.getReynoldsNumber()) < 2000) {
       return 3.66 / node.getHydraulicDiameter(phaseNum)
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getConductivity();
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getConductivity();
     }
     // if turbulent - use chilton colburn analogy
     else {
-      double temp = node.getBulkSystem().getPhases()[phaseNum].getCp()
-          / node.getBulkSystem().getPhases()[phaseNum].getMolarMass()
-          / node.getBulkSystem().getPhases()[phaseNum].getNumberOfMolesInPhase()
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getDensity()
+      double temp = node.getBulkSystem().getPhase(phaseNum).getCp()
+          / node.getBulkSystem().getPhase(phaseNum).getMolarMass()
+          / node.getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity()
           * node.getVelocity(phaseNum);
       return 0.5 * this.calcWallFrictionFactor(phaseNum, node) * Math.pow(prandtlNumber, -2.0 / 3.0)
           * temp;
@@ -77,8 +77,7 @@ public double calcWallMassTransferCoefficient(int phaseNum, double schmidtNumber
       FlowNodeInterface node) {
     if (Math.abs(node.getReynoldsNumber()) < 2000) {
       return 3.66 / node.getHydraulicDiameter(phaseNum) / schmidtNumber
-          * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-              .getKinematicViscosity();
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity();
     } else {
       double temp = node.getVelocity(phaseNum);
       return 0.5 * this.calcWallFrictionFactor(phaseNum, node) * Math.pow(schmidtNumber, -2.0 / 3.0)
@@ -98,8 +97,9 @@ public double calcInterphaseMassTransferCoefficient(int phaseNum, double schmidt
       // System.out.println("red gas " +
       // redMassTrans/Math.pow(node.getReynoldsNumber(phase),c3));
       // System.out.println("sc gas " + schmidtNumber);
-      massTrans = redMassTrans * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-          .getKinematicViscosity() / schmidtNumber / node.getGeometry().getDiameter();
+      massTrans = redMassTrans
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
+          / schmidtNumber / node.getGeometry().getDiameter();
     }
     if (phaseNum == 1) {
       double c2 = 0.181, c3 = 0.72, c4 = 0.33;
@@ -108,8 +108,9 @@ public double calcInterphaseMassTransferCoefficient(int phaseNum, double schmidt
       // System.out.println("red liq" +
       // redMassTrans/Math.pow(node.getReynoldsNumber(phase),c3));
       // System.out.println("sc liq " + schmidtNumber);
-      massTrans = redMassTrans * node.getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-          .getKinematicViscosity() / schmidtNumber / node.getGeometry().getDiameter();
+      massTrans = redMassTrans
+          * node.getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getKinematicViscosity()
+          / schmidtNumber / node.getGeometry().getDiameter();
     }
     return massTrans;
   }
diff --git a/src/main/java/neqsim/fluidmechanics/flowsolver/twophaseflowsolver/stirredcellsolver/StirredCellSolver.java b/src/main/java/neqsim/fluidmechanics/flowsolver/twophaseflowsolver/stirredcellsolver/StirredCellSolver.java
index 6bc26a568..95c81f082 100644
--- a/src/main/java/neqsim/fluidmechanics/flowsolver/twophaseflowsolver/stirredcellsolver/StirredCellSolver.java
+++ b/src/main/java/neqsim/fluidmechanics/flowsolver/twophaseflowsolver/stirredcellsolver/StirredCellSolver.java
@@ -214,7 +214,7 @@ public void initPressure(int phaseNum) {
     for (int i = 0; i < numberOfNodes; i++) {
       pipe.getNode(i).init();
       pipe.getNode(i).getBulkSystem()
-          .setPressure(0.8 * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getdPdrho()
+          .setPressure(0.8 * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getdPdrho()
               * diffMatrix.get(i, 0) * 1e-5 + pipe.getNode(i).getBulkSystem().getPressure());
       pipe.getNode(i).init();
     }
@@ -253,10 +253,9 @@ public void initTemperature(int phaseNum) {
       pipe.getNode(i).getBulkSystem()
           .setTemperature(
               pipe.getNode(i).getBulkSystem().getTemperature(phaseNum) + 0.8 * diffMatrix.get(i, 0)
-                  / (pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getCp()
-                      / pipe.getNode(i).getBulkSystem().getPhases()[phaseNum]
-                          .getNumberOfMolesInPhase()
-                      / pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getMolarMass()),
+                  / (pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getCp()
+                      / pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
+                      / pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getMolarMass()),
               phaseNum);
       pipe.getNode(i).init();
     }
@@ -288,39 +287,35 @@ public void initPhaseFraction(int phase) {
    */
   public void initComposition(int phaseNum, int comp) {
     for (int j = 0; j < numberOfNodes; j++) {
-      if ((pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp].getx()
-          + diffMatrix.get(j, 0)
-              * pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getMolarMass()
-              / pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp]
+      if ((pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp].getx()
+          + diffMatrix.get(j, 0) * pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getMolarMass()
+              / pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp]
                   .getMolarMass()) > 1.0) {
-        pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp]
-            .setx(1.0 - 1e-30);
-      } else if (pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp].getx()
-          + diffMatrix.get(j, 0)
-              * pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getMolarMass()
-              / pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp]
+        pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp].setx(1.0 - 1e-30);
+      } else if (pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp].getx()
+          + diffMatrix.get(j, 0) * pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getMolarMass()
+              / pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp]
                   .getMolarMass() < 0.0) {
-        pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp].setx(1e-30);
+        pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp].setx(1e-30);
       } else {
-        pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp]
-            .setx(pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp].getx()
+        pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp]
+            .setx(pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp].getx()
                 + diffMatrix.get(j, 0)
-                    * pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getMolarMass()
-                    / pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp]
+                    * pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getMolarMass()
+                    / pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp]
                         .getMolarMass());
         // pipe.getNode(j).getBulkSystem().getPhases()[0].getComponent(p).getx() +
         // 0.5*diff4Matrix[p].get(j,0));
       }
 
       double xSum = 0.0;
-      for (int i = 0; i < pipe.getNode(j).getBulkSystem().getPhases()[phaseNum]
-          .getNumberOfComponents() - 1; i++) {
-        xSum += pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponent(i).getx();
+      for (int i = 0; i < pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getNumberOfComponents()
+          - 1; i++) {
+        xSum += pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponent(i).getx();
       }
 
-      pipe.getNode(j).getBulkSystem().getPhases()[phaseNum]
-          .getComponents()[pipe.getNode(j).getBulkSystem().getPhases()[phaseNum]
-              .getNumberOfComponents() - 1].setx(1.0 - xSum);
+      pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[pipe.getNode(j)
+          .getBulkSystem().getPhase(phaseNum).getNumberOfComponents() - 1].setx(1.0 - xSum);
       pipe.getNode(j).init();
     }
   }
diff --git a/src/main/java/neqsim/fluidmechanics/flowsolver/twophaseflowsolver/twophasepipeflowsolver/TwoPhaseFixedStaggeredGridSolver.java b/src/main/java/neqsim/fluidmechanics/flowsolver/twophaseflowsolver/twophasepipeflowsolver/TwoPhaseFixedStaggeredGridSolver.java
index d5e1c3101..47cd5d93f 100644
--- a/src/main/java/neqsim/fluidmechanics/flowsolver/twophaseflowsolver/twophasepipeflowsolver/TwoPhaseFixedStaggeredGridSolver.java
+++ b/src/main/java/neqsim/fluidmechanics/flowsolver/twophaseflowsolver/twophasepipeflowsolver/TwoPhaseFixedStaggeredGridSolver.java
@@ -235,7 +235,7 @@ public void initPressure(int phaseNum) {
     for (int i = 0; i < numberOfNodes; i++) {
       pipe.getNode(i).init();
       pipe.getNode(i).getBulkSystem()
-          .setPressure(0.8 * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getdPdrho()
+          .setPressure(0.8 * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getdPdrho()
               * diffMatrix.get(i, 0) * 1e-5 + pipe.getNode(i).getBulkSystem().getPressure());
       pipe.getNode(i).init();
     }
@@ -274,10 +274,9 @@ public void initTemperature(int phaseNum) {
       pipe.getNode(i).getBulkSystem()
           .setTemperature(
               pipe.getNode(i).getBulkSystem().getTemperature(phaseNum) + 0.8 * diffMatrix.get(i, 0)
-                  / (pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getCp()
-                      / pipe.getNode(i).getBulkSystem().getPhases()[phaseNum]
-                          .getNumberOfMolesInPhase()
-                      / pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getMolarMass()),
+                  / (pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getCp()
+                      / pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
+                      / pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getMolarMass()),
               phaseNum);
       pipe.getNode(i).init();
     }
@@ -309,39 +308,35 @@ public void initPhaseFraction(int phase) {
    */
   public void initComposition(int phaseNum, int comp) {
     for (int j = 0; j < numberOfNodes; j++) {
-      if ((pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp].getx()
-          + diffMatrix.get(j, 0)
-              * pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getMolarMass()
-              / pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp]
+      if ((pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp].getx()
+          + diffMatrix.get(j, 0) * pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getMolarMass()
+              / pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp]
                   .getMolarMass()) > 1.0) {
-        pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp]
-            .setx(1.0 - 1e-30);
-      } else if (pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp].getx()
-          + diffMatrix.get(j, 0)
-              * pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getMolarMass()
-              / pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp]
+        pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp].setx(1.0 - 1e-30);
+      } else if (pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp].getx()
+          + diffMatrix.get(j, 0) * pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getMolarMass()
+              / pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp]
                   .getMolarMass() < 0.0) {
-        pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp].setx(1e-30);
+        pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp].setx(1e-30);
       } else {
-        pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp]
-            .setx(pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp].getx()
+        pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp]
+            .setx(pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp].getx()
                 + diffMatrix.get(j, 0)
-                    * pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getMolarMass()
-                    / pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponents()[comp]
+                    * pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getMolarMass()
+                    / pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[comp]
                         .getMolarMass());
         // pipe.getNode(j).getBulkSystem().getPhases()[0].getComponent(p).getx()
         // + 0.5*diff4Matrix[p].get(j,0));
       }
 
       double xSum = 0.0;
-      for (int i = 0; i < pipe.getNode(j).getBulkSystem().getPhases()[phaseNum]
-          .getNumberOfComponents() - 1; i++) {
-        xSum += pipe.getNode(j).getBulkSystem().getPhases()[phaseNum].getComponent(i).getx();
+      for (int i = 0; i < pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getNumberOfComponents()
+          - 1; i++) {
+        xSum += pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponent(i).getx();
       }
 
-      pipe.getNode(j).getBulkSystem().getPhases()[phaseNum]
-          .getComponents()[pipe.getNode(j).getBulkSystem().getPhases()[phaseNum]
-              .getNumberOfComponents() - 1].setx(1.0 - xSum);
+      pipe.getNode(j).getBulkSystem().getPhase(phaseNum).getComponents()[pipe.getNode(j)
+          .getBulkSystem().getPhase(phaseNum).getNumberOfComponents() - 1].setx(1.0 - xSum);
       pipe.getNode(j).init();
     }
   }
@@ -359,8 +354,7 @@ public void setMassConservationMatrix(int phaseNum) {
       a[0] = 0;
       b[0] = 1.0;
       c[0] = 0;
-      SU = pipe.getNode(0).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-          .getDensity();
+      SU = pipe.getNode(0).getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity();
       r[0] = SU;
     } else {
       // double Ae = pipe.getNode(0).getArea(phase);
@@ -467,11 +461,9 @@ public void setPhaseFractionMatrix(int phaseNum) {
       double Ae = pipe.getNode(i).getGeometry().getArea();
       double Aw = pipe.getNode(i - 1).getGeometry().getArea();
       double Fe = pipe.getNode(i).getVelocityOut(phaseNum).doubleValue() * Ae
-          * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-              .getDensity();
-      double Fw = pipe.getNode(i).getVelocityIn(phaseNum).doubleValue() * Aw
-          * pipe.getNode(i - 1).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-              .getDensity();
+          * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity();
+      double Fw = pipe.getNode(i).getVelocityIn(phaseNum).doubleValue() * Aw * pipe.getNode(i - 1)
+          .getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity();
 
       if (dynamic) {
         oldMass[phaseNum][i] = 1.0 / timeStep * pipe.getNode(i).getGeometry().getArea()
@@ -495,11 +487,9 @@ public void setPhaseFractionMatrix(int phaseNum) {
     double Aw = pipe.getNode(i - 1).getGeometry().getArea();
 
     double Fe = pipe.getNode(i).getVelocity(phaseNum) * Ae
-        * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-            .getDensity();
-    double Fw = pipe.getNode(i).getVelocityIn(phaseNum).doubleValue() * Aw
-        * pipe.getNode(i - 1).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-            .getDensity();
+        * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity();
+    double Fw = pipe.getNode(i).getVelocityIn(phaseNum).doubleValue() * Aw * pipe.getNode(i - 1)
+        .getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity();
 
     if (dynamic) {
       oldMass[phaseNum][i] = 1.0 / timeStep * pipe.getNode(i).getGeometry().getArea()
@@ -552,7 +542,7 @@ public void setImpulsMatrixTDMA(int phaseNum) {
       double Aw = pipe.getNode(i - 1).getArea(phaseNum);
       double Amean = pipe.getNode(i - 1).getArea(phaseNum);
       double meanFrik = pipe.getNode(i - 1).getWallFrictionFactor(phaseNum);
-      double meanDensity = pipe.getNode(i - 1).getBulkSystem().getPhases()[phaseNum]
+      double meanDensity = pipe.getNode(i - 1).getBulkSystem().getPhase(phaseNum)
           .getPhysicalProperties().getDensity();
       double oldMeanDensity = oldDensity[phaseNum][i];
       double meanVelocity = pipe.getNode(i - 1).getVelocity(phaseNum);
@@ -583,10 +573,11 @@ public void setImpulsMatrixTDMA(int phaseNum) {
       // System.out.println("su " + SU);
       SP = -pipe.getNode(i - 1).getWallContactLength(phaseNum) * nodeLength * meanDensity * meanFrik
           * meanVelocity / 4.0;
-      Fw = Aw * pipe.getNode(i - 1).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
+      Fw = Aw * pipe.getNode(i - 1).getBulkSystem().getPhase(phaseNum).getPhysicalProperties()
           .getDensity() * pipe.getNode(i - 1).getVelocity(phaseNum);
-      Fe = Ae * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-          .getDensity() * pipe.getNode(i).getVelocity(phaseNum);
+      Fe = Ae
+          * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity()
+          * pipe.getNode(i).getVelocity(phaseNum);
 
       if (dynamic) {
         oldImpuls[phaseNum][i] = 1.0 / timeStep * oldMeanDensity * nodeLength * Amean;
@@ -609,8 +600,8 @@ public void setImpulsMatrixTDMA(int phaseNum) {
     double Amean = pipe.getNode(i - 1).getArea(phaseNum);
     // double meanDiameter = pipe.getNode(i - 1).getGeometry().getDiameter();
     double meanFrik = pipe.getNode(i - 1).getWallFrictionFactor(phaseNum);
-    double meanDensity = pipe.getNode(i - 1).getBulkSystem().getPhases()[phaseNum]
-        .getPhysicalProperties().getDensity();
+    double meanDensity =
+        pipe.getNode(i - 1).getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity();
     double oldMeanDensity = oldDensity[phaseNum][i];
     double meanVelocity = pipe.getNode(i - 1).getVelocity(phaseNum);
     double vertposchange = pipe.getNode(i).getVerticalPositionOfNode()
@@ -631,10 +622,10 @@ public void setImpulsMatrixTDMA(int phaseNum) {
             * (pipe.getNode(i).getVelocity(0) - pipe.getNode(i).getVelocity(1)) / 8.0 * sign;
     SP = -pipe.getNode(i - 1).getWallContactLength(phaseNum) * nodeLength * meanDensity * meanFrik
         * meanVelocity / 4.0;
-    Fw = Aw * pipe.getNode(i - 1).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
+    Fw = Aw * pipe.getNode(i - 1).getBulkSystem().getPhase(phaseNum).getPhysicalProperties()
         .getDensity() * pipe.getNode(i).getVelocityIn(phaseNum).doubleValue();
     Fe = Ae
-        * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getDensity()
+        * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity()
         * pipe.getNode(i).getVelocity(phaseNum);
 
     if (dynamic) {
@@ -667,9 +658,9 @@ public void setEnergyMatrixTDMA(int phaseNum) {
     a[0] = 0;
     b[0] = 1.0;
     c[0] = 0;
-    double SU = pipe.getNode(0).getBulkSystem().getPhases()[phaseNum].getEnthalpy()
-        / pipe.getNode(0).getBulkSystem().getPhases()[phaseNum].getNumberOfMolesInPhase()
-        / pipe.getNode(0).getBulkSystem().getPhases()[phaseNum].getMolarMass();
+    double SU = pipe.getNode(0).getBulkSystem().getPhase(phaseNum).getEnthalpy()
+        / pipe.getNode(0).getBulkSystem().getPhase(phaseNum).getNumberOfMolesInPhase()
+        / pipe.getNode(0).getBulkSystem().getPhase(phaseNum).getMolarMass();
     r[0] = SU;
 
     for (int i = 1; i < numberOfNodes - 1; i++) {
@@ -688,10 +679,10 @@ public void setEnergyMatrixTDMA(int phaseNum) {
               - pipe.getNode(i - 1).getVerticalPositionOfNode());
 
       SU = -pipe.getNode(i).getArea(phaseNum) * gravity
-          * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getDensity()
+          * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getDensity()
           * pipe.getNode(i).getVelocity(phaseNum) * vertposchange
           + pipe.getNode(i).getArea(phaseNum) * 4.0 * 0.02
-              * (278.0 - pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getTemperature())
+              * (278.0 - pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getTemperature())
               / (pipe.getNode(i).getGeometry().getDiameter())
               * pipe.getNode(i).getGeometry().getNodeLength()
           + sign * pipe.getNode(i).getFluidBoundary().getInterphaseHeatFlux(phaseNum)
@@ -702,10 +693,11 @@ public void setEnergyMatrixTDMA(int phaseNum) {
       double SP = 0;
 
       double Fw =
-          Aw * pipe.getNode(i - 1).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
+          Aw * pipe.getNode(i - 1).getBulkSystem().getPhase(phaseNum).getPhysicalProperties()
               .getDensity() * pipe.getNode(i).getVelocityIn(phaseNum).doubleValue();
-      double Fe = Ae * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-          .getDensity() * pipe.getNode(i).getVelocityOut(phaseNum).doubleValue();
+      double Fe = Ae
+          * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity()
+          * pipe.getNode(i).getVelocityOut(phaseNum).doubleValue();
 
       if (dynamic) {
         oldEnergy[phaseNum][i] = 1.0 / timeStep * oldDensity[phaseNum][i]
@@ -739,12 +731,12 @@ public void setEnergyMatrixTDMA(int phaseNum) {
         - pipe.getNode(i - 1).getVerticalPositionOfNode());
 
     SU = -pipe.getNode(i).getArea(phaseNum) * gravity
-        * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getDensity()
+        * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getDensity()
         * pipe.getNode(i).getVelocity(phaseNum) * vertposchange
         + pipe.getNode(i).getWallContactLength(phaseNum)
             / pipe.getNode(i).getGeometry().getCircumference() * pipe.getNode(i).getArea(phaseNum)
             * 4.0 * 0.02
-            * (278.0 - pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getTemperature())
+            * (278.0 - pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getTemperature())
             / (pipe.getNode(i).getGeometry().getDiameter())
             * pipe.getNode(i).getGeometry().getNodeLength()
         + sign * pipe.getNode(i).getFluidBoundary().getInterphaseHeatFlux(phaseNum)
@@ -756,12 +748,11 @@ public void setEnergyMatrixTDMA(int phaseNum) {
     // -pipe.getNode(i).getGeometry().getArea() * 4.0 * 12.0 /
     // (pipe.getNode(i).getGeometry().getDiameter())*pipe.getNode(i).getGeometry().getNodeLength();
 
-    double Fw =
-        Aw * pipe.getNode(i - 1).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-            .getDensity() * pipe.getNode(i).getVelocityIn(phaseNum).doubleValue();
-    double Fe = Ae
-        * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties().getDensity()
-        * pipe.getNode(i).getVelocity(phaseNum);
+    double Fw = Aw * pipe.getNode(i - 1).getBulkSystem().getPhase(phaseNum).getPhysicalProperties()
+        .getDensity() * pipe.getNode(i).getVelocityIn(phaseNum).doubleValue();
+    double Fe =
+        Ae * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity()
+            * pipe.getNode(i).getVelocity(phaseNum);
 
     if (dynamic) {
       oldEnergy[phaseNum][i] = 1.0 / timeStep * oldDensity[phaseNum][i]
@@ -792,8 +783,7 @@ public void setComponentConservationMatrix2(int phaseNum, int componentNumber) {
     a[0] = 0;
     b[0] = 1.0;
     c[0] = 0;
-    SU = pipe.getNode(0).getBulkSystem().getPhases()[phaseNum].getComponents()[componentNumber]
-        .getx();
+    SU = pipe.getNode(0).getBulkSystem().getPhase(phaseNum).getComponents()[componentNumber].getx();
     // System.out.println("phase x0: "
     // +pipe.getNode(0).getBulkSystem().getPhases()[0].getComponents()[componentNumber].getx());
     // System.out.println("phase x1: "
@@ -813,17 +803,16 @@ public void setComponentConservationMatrix2(int phaseNum, int componentNumber) {
           + fw / pipe.getNode(i - 1).getArea(phaseNum));
 
       double Fe = pipe.getNode(i).getVelocityOut(phaseNum).doubleValue()
-          * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
-              .getDensity()
-          * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getComponents()[componentNumber]
+          * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getPhysicalProperties().getDensity()
+          * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getComponents()[componentNumber]
               .getMolarMass()
-          / pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getMolarMass() * Ae;
+          / pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getMolarMass() * Ae;
       double Fw = pipe.getNode(i).getVelocityIn(phaseNum).doubleValue()
-          * pipe.getNode(i - 1).getBulkSystem().getPhases()[phaseNum].getPhysicalProperties()
+          * pipe.getNode(i - 1).getBulkSystem().getPhase(phaseNum).getPhysicalProperties()
               .getDensity()
-          * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getComponents()[componentNumber]
+          * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getComponents()[componentNumber]
               .getMolarMass()
-          / pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getMolarMass() * Aw;
+          / pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getMolarMass() * Aw;
 
       // System.out.println("vel: " +
       // pipe.getNode(i).getVelocityOut(phase).doubleValue() + " fe " + Fe);
@@ -842,12 +831,12 @@ public void setComponentConservationMatrix2(int phaseNum, int componentNumber) {
     a[numberOfNodes - 1] = -1.0; // -1.0;
     b[numberOfNodes - 1] = 1.0;
     c[numberOfNodes - 1] = 0;
-    SU = pipe.getNode(numberOfNodes - 2).getBulkSystem().getPhases()[phaseNum]
-        .getPhysicalProperties().getDensity()
+    SU = pipe.getNode(numberOfNodes - 2).getBulkSystem().getPhase(phaseNum).getPhysicalProperties()
+        .getDensity()
         * pipe.getNode(numberOfNodes - 2).getVelocityIn(phaseNum).doubleValue()
-        * pipe.getNode(numberOfNodes - 2).getBulkSystem().getPhases()[phaseNum]
+        * pipe.getNode(numberOfNodes - 2).getBulkSystem().getPhase(phaseNum)
             .getComponents()[componentNumber].getx()
-        / (pipe.getNode(numberOfNodes - 1).getBulkSystem().getPhases()[phaseNum]
+        / (pipe.getNode(numberOfNodes - 1).getBulkSystem().getPhase(phaseNum)
             .getPhysicalProperties().getDensity()
             * pipe.getNode(numberOfNodes - 1).getVelocityIn(phaseNum).doubleValue());
     r[numberOfNodes - 1] = 0; // SU;
@@ -867,11 +856,10 @@ public void setComponentConservationMatrix(int phaseNum, int componentNumber) {
     a[0] = 0;
     b[0] = 1.0;
     c[0] = 0;
-    SU = pipe.getNode(0).getBulkSystem().getPhases()[phaseNum].getComponents()[componentNumber]
-        .getx()
-        * pipe.getNode(0).getBulkSystem().getPhases()[phaseNum].getComponents()[componentNumber]
+    SU = pipe.getNode(0).getBulkSystem().getPhase(phaseNum).getComponents()[componentNumber].getx()
+        * pipe.getNode(0).getBulkSystem().getPhase(phaseNum).getComponents()[componentNumber]
             .getMolarMass()
-        / pipe.getNode(0).getBulkSystem().getPhases()[phaseNum].getMolarMass();
+        / pipe.getNode(0).getBulkSystem().getPhase(phaseNum).getMolarMass();
     r[0] = SU;
 
     for (int i = 1; i < numberOfNodes - 1; i++) {
@@ -879,21 +867,21 @@ public void setComponentConservationMatrix(int phaseNum, int componentNumber) {
       double Aw = pipe.getNode(i - 1).getArea(phaseNum);
 
       double Fe = pipe.getNode(i).getVelocityOut(phaseNum).doubleValue()
-          * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getDensity() * Ae;
+          * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getDensity() * Ae;
       double Fw = pipe.getNode(i).getVelocityIn(phaseNum).doubleValue()
-          * pipe.getNode(i - 1).getBulkSystem().getPhases()[phaseNum].getDensity() * Aw;
+          * pipe.getNode(i - 1).getBulkSystem().getPhase(phaseNum).getDensity() * Aw;
       // System.out.println("vel: " +
       // pipe.getNode(i).getVelocityOut(phase).doubleValue() + " fe " + Fe);
       if (dynamic) {
         oldComp[phaseNum][i] = 1.0 / timeStep * pipe.getNode(i).getArea(phaseNum)
             * pipe.getNode(i).getGeometry().getNodeLength()
-            * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getDensity();
+            * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getDensity();
       } else {
         oldComp[phaseNum][i] = 0.0;
       }
 
       SU = +sign * pipe.getNode(i).getFluidBoundary().getInterphaseMolarFlux(componentNumber)
-          * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getComponents()[componentNumber]
+          * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getComponents()[componentNumber]
               .getMolarMass()
           * pipe.getNode(i).getGeometry().getNodeLength()
           * pipe.getNode(i).getInterphaseContactLength(phaseNum)
@@ -922,12 +910,12 @@ public void setComponentConservationMatrix(int phaseNum, int componentNumber) {
     // double Aw = pipe.getNode(i - 1).getArea(phaseNum);
 
     // double Fe = pipe.getNode(i).getVelocity(phaseNum) *
-    // pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getDensity() * Ae;
+    // pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getDensity() * Ae;
     // double Fw = pipe.getNode(i).getVelocityIn(phaseNum).doubleValue() * pipe.getNode(i -
-    // 1).getBulkSystem().getPhases()[phaseNum].getDensity() * Aw;
+    // 1).getBulkSystem().getPhase(phaseNum).getDensity() * Aw;
 
     SU = +sign * pipe.getNode(i).getFluidBoundary().getInterphaseMolarFlux(componentNumber)
-        * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getComponents()[componentNumber]
+        * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getComponents()[componentNumber]
             .getMolarMass()
         * pipe.getNode(i).getGeometry().getNodeLength()
         * pipe.getNode(i).getInterphaseContactLength(phaseNum)
@@ -939,7 +927,7 @@ public void setComponentConservationMatrix(int phaseNum, int componentNumber) {
     if (dynamic) {
       oldComp[phaseNum][i] = 1.0 / timeStep * pipe.getNode(i).getArea(phaseNum)
           * pipe.getNode(i).getGeometry().getNodeLength()
-          * pipe.getNode(i).getBulkSystem().getPhases()[phaseNum].getDensity();
+          * pipe.getNode(i).getBulkSystem().getPhase(phaseNum).getDensity();
     } else {
       oldComp[phaseNum][i] = 0.0;
     }
diff --git a/src/main/java/neqsim/thermo/ThermodynamicModelTest.java b/src/main/java/neqsim/thermo/ThermodynamicModelTest.java
index 70931b97d..728f74912 100644
--- a/src/main/java/neqsim/thermo/ThermodynamicModelTest.java
+++ b/src/main/java/neqsim/thermo/ThermodynamicModelTest.java
@@ -362,7 +362,7 @@ public boolean checkNumerically() {
             .info("component name " + system.getPhase(phaseNum).getComponent(k).getComponentName());
         logger.info("dn " + dn);
         // logger.info("component name " +
-        // system.getPhases()[phaseNum].getComponent(k).getComponentName());
+        // system.getPhase(phaseNum).getComponent(k).getComponentName());
         // logger.info("dn " + dn);
         if (dn < 1e-12) {
           dn = 1e-12;
diff --git a/src/main/java/neqsim/thermodynamicoperations/flashops/TPflash.java b/src/main/java/neqsim/thermodynamicoperations/flashops/TPflash.java
index 0b65d6ce3..eeb6aca42 100644
--- a/src/main/java/neqsim/thermodynamicoperations/flashops/TPflash.java
+++ b/src/main/java/neqsim/thermodynamicoperations/flashops/TPflash.java
@@ -328,10 +328,11 @@ public void run() {
           }
         } else if (tpdy < 0) {
           for (i = 0; i < system.getPhase(0).getNumberOfComponents(); i++) {
-            system.getPhase(0).getComponent(i).setK(Math
-                .exp(minGibsLogFugCoef[i]
-                    - Math.log(system.getPhase(0).getComponent(i).getFugacityCoefficient()))
-                * presdiff);
+            system.getPhase(0).getComponent(i)
+                .setK(Math
+                    .exp(minGibsLogFugCoef[i]
+                        - Math.log(system.getPhase(0).getComponent(i).getFugacityCoefficient()))
+                    * presdiff);
             system.getPhase(1).getComponent(i).setK(system.getPhase(0).getComponent(i).getK());
           }
         } else {