diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
index 99f14113f3..885b1f080e 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
@@ -11,8 +11,9 @@
     PtrList<radiation::radiationModel> radiation(fluidRegions.size());
 
     List<scalar> initialMassFluid(fluidRegions.size());
-    List<label> pRefCellFluid(fluidRegions.size(),0);
-    List<scalar> pRefValueFluid(fluidRegions.size(),0.0);
+    List<label> pRefCellFluid(fluidRegions.size(), 0);
+    List<scalar> pRefValueFluid(fluidRegions.size(), 0.0);
+    List<bool> frozenFlowFluid(fluidRegions.size(), false);
 
     PtrList<dimensionedScalar> rhoMax(fluidRegions.size());
     PtrList<dimensionedScalar> rhoMin(fluidRegions.size());
@@ -160,38 +161,22 @@
 
         initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
 
+        const dictionary& simpleDict =
+            fluidRegions[i].solutionDict().subDict("SIMPLE");
+
         setRefCell
         (
             thermoFluid[i].p(),
             p_rghFluid[i],
-            fluidRegions[i].solutionDict().subDict("SIMPLE"),
+            simpleDict,
             pRefCellFluid[i],
             pRefValueFluid[i]
         );
 
-        rhoMax.set
-        (
-            i,
-            new dimensionedScalar
-            (
-                fluidRegions[i].solutionDict().subDict("SIMPLE").lookup
-                (
-                    "rhoMax"
-                )
-            )
-        );
+        simpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
 
-        rhoMin.set
-        (
-            i,
-            new dimensionedScalar
-            (
-                fluidRegions[i].solutionDict().subDict("SIMPLE").lookup
-                (
-                    "rhoMin"
-                )
-            )
-        );
+        rhoMax.set(i, new dimensionedScalar(simpleDict.lookup("rhoMax")));
+        rhoMin.set(i, new dimensionedScalar(simpleDict.lookup("rhoMin")));
 
         Info<< "    Adding fvOptions\n" << endl;
         fluidFvOptions.set
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
index 34aea7643e..fe5ce1e823 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
@@ -25,6 +25,7 @@
 
     const label pRefCell = pRefCellFluid[i];
     const scalar pRefValue = pRefValueFluid[i];
+    const bool frozenFlow = frozenFlowFluid[i];
 
     volScalarField& p_rgh = p_rghFluid[i];
     const volScalarField& gh = ghFluid[i];
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H
index a17893fb83..b43a949043 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H
@@ -1,11 +1,19 @@
 //  Pressure-velocity SIMPLE corrector
 
-    p_rgh.storePrevIter();
-    rho.storePrevIter();
     {
-        #include "UEqn.H"
-        #include "EEqn.H"
-        #include "pEqn.H"
-    }
+        if (frozenFlow)
+        {
+            #include "EEqn.H"
+        }
+        else
+        {
+            p_rgh.storePrevIter();
+            rho.storePrevIter();
 
-    turb.correct();
+            #include "UEqn.H"
+            #include "EEqn.H"
+            #include "pEqn.H"
+
+            turb.correct();
+        }
+    }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index 9afa705c11..98ef8d7f61 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -13,6 +13,7 @@
     PtrList<volScalarField> dpdtFluid(fluidRegions.size());
 
     List<scalar> initialMassFluid(fluidRegions.size());
+    List<bool> frozenFlowFluid(fluidRegions.size(), false);
 
     PtrList<fv::IOoptionList> fluidFvOptions(fluidRegions.size());
 
@@ -102,7 +103,7 @@
             )
         );
 
-        Info<< "    Adding turbulence \n" << endl;
+        Info<< "    Adding to turbulence\n" << endl;
         turbulence.set
         (
             i,
@@ -190,6 +191,10 @@
             )
         );
 
+        const dictionary& pimpleDict =
+            fluidRegions[i].solutionDict().subDict("PIMPLE");
+        pimpleDict.readIfPresent("frozenFlow", frozenFlowFluid[i]);
+
         Info<< "    Adding fvOptions\n" << endl;
         fluidFvOptions.set
         (
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 867585f43a..ef3a3df78f 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -28,3 +28,5 @@
         dimMass,
         initialMassFluid[i]
     );
+
+    const bool frozenFlow = frozenFlowFluid[i];
\ No newline at end of file
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
index 354c4a137c..6f511e0292 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
@@ -3,24 +3,31 @@ if (finalIter)
     mesh.data::add("finalIteration", true);
 }
 
-if (oCorr == 0)
+if (frozenFlow)
 {
-    #include "rhoEqn.H"
+    #include "EEqn.H"
 }
-
-#include "UEqn.H"
-#include "EEqn.H"
-
-// --- PISO loop
-for (int corr=0; corr<nCorr; corr++)
+else
 {
-    #include "pEqn.H"
+    if (oCorr == 0)
+    {
+        #include "rhoEqn.H"
+    }
+
+    #include "UEqn.H"
+    #include "EEqn.H"
+
+    // --- PISO loop
+    for (int corr=0; corr<nCorr; corr++)
+    {
+        #include "pEqn.H"
+    }
+
+    turb.correct();
+
+    rho = thermo.rho();
 }
 
-turb.correct();
-
-rho = thermo.rho();
-
 if (finalIter)
 {
     mesh.data::remove("finalIteration");