diff --git a/src/regionFaModels/KirchhoffShell/KirchhoffShell.C b/src/regionFaModels/KirchhoffShell/KirchhoffShell.C
index 0b8cfd0a98dd65afb27e7bca54ddf222433516b6..d19871c2fe9d6f414d8e7bce4e1a2c85e0f5396a 100644
--- a/src/regionFaModels/KirchhoffShell/KirchhoffShell.C
+++ b/src/regionFaModels/KirchhoffShell/KirchhoffShell.C
@@ -49,7 +49,6 @@ bool KirchhoffShell::init(const dictionary& dict)
     return true;
 }
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 void KirchhoffShell::solveDisplacement()
 {
@@ -57,20 +56,27 @@ void KirchhoffShell::solveDisplacement()
 
     const Time& time = primaryMesh().time();
 
-    areaScalarField solidMass(rho()*h_);
-    areaScalarField solidD(D()/solidMass);
+    // Create operand fields for solid physics
+    const areaScalarField solidMass(rho()*h_);
+    const areaScalarField solidD(D()/solidMass);
 
-    // Save old times
-    areaScalarField w0(w_.oldTime());
-    areaScalarField w00(w_.oldTime().oldTime());
+    // Deep copy old times of shell displacement field
+    auto tw0 = tmp<areaScalarField>::New(w_.oldTime());
+    auto tw00 = tmp<areaScalarField>::New(w_.oldTime().oldTime());
 
+    // Flag terms to avoid redundant time-derivative calculations
+    const bool f0_enabled = (f0_.value() != scalar(0));
+    const bool f1_enabled = (f1_.value() != scalar(0));
+    const bool f2_enabled = (f2_.value() != scalar(0));
+
+    // Restore various old fields in sub-cycling, if need be
     if (nSubCycles_ > 1)
     {
-        // Restore the oldTime in sub-cycling
         w_.oldTime() = w0_;
         w_.oldTime().oldTime() = w00_;
-        laplaceW_.oldTime() = laplaceW0_;
-        laplace2W_.oldTime() = laplace2W0_;
+
+        if (f0_enabled) laplaceW_.oldTime() = laplaceW0_;
+        if (f2_enabled) laplace2W_.oldTime() = laplace2W0_;
     }
 
     for
@@ -89,25 +95,29 @@ void KirchhoffShell::solveDisplacement()
         faScalarMatrix wEqn
         (
             fam::d2dt2(w_)
-         +  f1_*fam::ddt(w_)
-         -  f0_*sqrt(solidD)*fac::ddt(laplaceW_)
-         +  solidD*(laplace2W_ + f2_*fac::ddt(laplace2W_))
+          + solidD*laplace2W_
         ==
             ps_/solidMass
           + faOptions()(solidMass, w_, dimLength/sqr(dimTime))
         );
 
+        // Avoid time-derivative calculations for f terms, if possible
+        if (f0_enabled) wEqn -= f0_*sqrt(solidD)*fac::ddt(laplaceW_);
+        if (f1_enabled) wEqn += f1_*fam::ddt(w_);
+        if (f2_enabled) wEqn += f2_*solidD*fac::ddt(laplace2W_);
+
         faOptions().constrain(wEqn);
 
         wEqn.solve();
 
+        // Cache various old fields inside the sub-cycling
         if (wSubCycle.index() >= wSubCycle.nSubCycles())
         {
-            // Cache oldTimes inside the sub-cycling
             w0_ = w_.oldTime();
             w00_ = w_.oldTime().oldTime();
-            laplaceW0_ = laplaceW_.oldTime();
-            laplace2W0_ = laplace2W_.oldTime();
+
+            if (f0_enabled) laplaceW0_ = laplaceW_.oldTime();
+            if (f2_enabled) laplace2W0_ = laplace2W_.oldTime();
 
             // Update shell acceleration
             a_ = fac::d2dt2(w_);
@@ -116,9 +126,9 @@ void KirchhoffShell::solveDisplacement()
 
     Info<< w_.name() << " min/max   = " << gMinMax(w_) << endl;
 
-    // Restore old time in main time
-    w_.oldTime() = w0;
-    w_.oldTime().oldTime() = w00;
+    // Steal the deep-copy of old times to restore the shell displacement
+    w_.oldTime() = tw0;
+    w_.oldTime().oldTime() = tw00;
 
     faOptions().correct(w_);
 }
@@ -134,11 +144,6 @@ KirchhoffShell::KirchhoffShell
 )
 :
     vibrationShellModel(modelType, mesh, dict),
-    f0_("f0", dimless, dict),
-    f1_("f1", inv(dimTime), dict),
-    f2_("f2", dimTime, dict),
-    nNonOrthCorr_(1),
-    nSubCycles_(1),
     ps_
     (
         IOobject
@@ -241,7 +246,12 @@ KirchhoffShell::KirchhoffShell
         ),
         regionMesh(),
         dimensionedScalar(inv(pow3(dimLength)), Zero)
-    )
+    ),
+    f0_("f0", dimless, dict),
+    f1_("f1", inv(dimTime), dict),
+    f2_("f2", dimTime, dict),
+    nNonOrthCorr_(1),
+    nSubCycles_(1)
 {
     init(dict);
 }
@@ -254,10 +264,10 @@ void KirchhoffShell::preEvolveRegion()
 
 void KirchhoffShell::evolveRegion()
 {
-    nNonOrthCorr_ = solution().get<label>("nNonOrthCorr");
-    nSubCycles_ = solution().get<label>("nSubCycles");
+    nNonOrthCorr_ = solution().getLabel("nNonOrthCorr");
+    nSubCycles_ = solution().getLabel("nSubCycles");
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr_; ++nonOrth)
+    for (int nonOrth=0; nonOrth <= nNonOrthCorr_; ++nonOrth)
     {
         solveDisplacement();
     }
diff --git a/src/regionFaModels/KirchhoffShell/KirchhoffShell.H b/src/regionFaModels/KirchhoffShell/KirchhoffShell.H
index 38460fb83ea74cd57c28c867eaf9eece58b36a7c..cf09f8dd670ceb59f646681d9af01934c7007bee 100644
--- a/src/regionFaModels/KirchhoffShell/KirchhoffShell.H
+++ b/src/regionFaModels/KirchhoffShell/KirchhoffShell.H
@@ -88,64 +88,54 @@ class KirchhoffShell
 {
     // Private Data
 
-        //- Damping coefficients [1/s]
-        dimensionedScalar f0_;
-        dimensionedScalar f1_;
-        dimensionedScalar f2_;
-
-
-    // Private Member Functions
-
-        //- Initialise KirchhoffShell
-        bool init(const dictionary& dict);
-
-
-protected:
-
-    // Protected Data
+        // Source term fields
 
-        // Solution parameters
+        //- External surface source [Pa]
+        const areaScalarField ps_;
 
-            //- Number of non orthogonal correctors
-            label nNonOrthCorr_;
+        //- Thickness [m]
+        areaScalarField h_;
 
-            //- Sub cycles
-            label nSubCycles_;
+        //- Laplace of the displacement
+        areaScalarField laplaceW_;
 
+        //- Laplace of the Laplace for the displacement
+        areaScalarField laplace2W_;
 
-        // Source term fields
+        //- Cache w.oldTime() in sub-cycling
+        areaScalarField w0_;
 
-            //- External surface source [Pa]
-            areaScalarField ps_;
+        //- Cache w.oldTime.oldTime() in sub-cycling
+        areaScalarField w00_;
 
-            //- Thickness [m]
-            areaScalarField h_;
+        //- Cache laplaceW.oldTime() in sub-cycling
+        areaScalarField laplaceW0_;
 
-            //- Laplace of the displacement
-            areaScalarField laplaceW_;
+        //- Cache laplace2.oldTime() in sub-cycling
+        areaScalarField laplace2W0_;
 
-            //- Laplace of the Laplace for the displacement
-            areaScalarField laplace2W_;
 
-            //- Cache w.oldTime() in sub-cycling
-            areaScalarField w0_;
+        // Solution parameters
 
-            //- Cache w.oldTime.oldTime() in sub-cycling
-            areaScalarField w00_;
+        //- Damping coefficients [1/s]
+        const dimensionedScalar f0_;
+        const dimensionedScalar f1_;
+        const dimensionedScalar f2_;
 
-            //- Cache laplaceW.oldTime() in sub-cycling
-            areaScalarField laplaceW0_;
+        //- Number of non orthogonal correctors
+        label nNonOrthCorr_;
 
-            //- Cache laplace2.oldTime() in sub-cycling
-            areaScalarField laplace2W0_;
+        //- Sub cycles
+        label nSubCycles_;
 
 
-    // Protected Member Functions
+    // Private Member Functions
 
-        // Equations
+        //- Initialise Kirchhoff shell model
+        bool init(const dictionary& dict);
 
-            //- Solve energy equation
-            void solveDisplacement();
+        //- Solve energy equation
+        void solveDisplacement();
 
 
 public:
diff --git a/src/regionFaModels/vibrationShellModel/vibrationShellModel.C b/src/regionFaModels/vibrationShellModel/vibrationShellModel.C
index bc952ad3a77ada5c7d555c2c925299c2b3ab8986..7ce1aaf7eecc5ac7a6bb98a6c1af4f1e705ab9f6 100644
--- a/src/regionFaModels/vibrationShellModel/vibrationShellModel.C
+++ b/src/regionFaModels/vibrationShellModel/vibrationShellModel.C
@@ -49,8 +49,6 @@ vibrationShellModel::vibrationShellModel
 )
 :
     regionFaModel(mesh, "vibratingShell", modelType, dict, true),
-    pName_(dict.get<word>("p")),
-    pa_(mesh.lookupObject<volScalarField>(pName_)),
     w_
     (
         IOobject
@@ -76,8 +74,10 @@ vibrationShellModel::vibrationShellModel
         regionMesh(),
         dimensionedScalar(dimAcceleration, Zero)
     ),
-    faOptions_(Foam::fa::options::New(mesh)),
-    solid_(dict.subDict("solid"))
+    solid_(dict.subDict("solid")),
+    pName_(dict.get<word>("p")),
+    pa_(mesh.lookupObject<volScalarField>(pName_)),
+    faOptions_(Foam::fa::options::New(mesh))
 {
     if (faOptions_.optionList::empty())
     {
diff --git a/src/regionFaModels/vibrationShellModel/vibrationShellModel.H b/src/regionFaModels/vibrationShellModel/vibrationShellModel.H
index c21c928d5f2f746113f89c92dcd27cb1d1b0e786..e8053ec526dd210ac1cd5674d9f8312397586476 100644
--- a/src/regionFaModels/vibrationShellModel/vibrationShellModel.H
+++ b/src/regionFaModels/vibrationShellModel/vibrationShellModel.H
@@ -93,24 +93,24 @@ protected:
 
     // Protected Data
 
-        //- Name of the coupled field in the primary region
-        word pName_;
-
-        //- Primary region acoustic pressure
-        const volScalarField& pa_;
-
         //- Shell displacement
         areaScalarField w_;
 
         //- Shell acceleration
         areaScalarField a_;
 
-        //- Pointer to faOptions
-        Foam::fa::options& faOptions_;
-
         //- Solid properties
         solidProperties solid_;
 
+        //- Name of the coupled field in the primary region
+        word pName_;
+
+        //- Primary region acoustic pressure
+        const volScalarField& pa_;
+
+        //- Reference to faOptions
+        Foam::fa::options& faOptions_;
+
 
 public: