bubble-dynamics

commit ad5f515d16442179b65474df548d82d0d594ba84
parent 58431e995aa4ea27b52bf1a0258d11f077310ba7
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date:   Fri,  9 Jun 2017 17:09:24 -0700

tidy-up KM based cases with correct bubble pressures and dpext issue

corrected coefficients in all cases for the forcing pressure derivative
term.  Tidy-up individual terms

Diffstat:
M
src/kernels/KMClusterPositions_D.h
 | 
21
+++++++++++----------
M
src/kernels/KMCluster_FC.h
 | 
18
+++++++++---------
M
src/kernels/KMCluster_TY.h
 | 
16
++++++++--------
M
src/kernels/KellerMiksis.h
 | 
22
++++++++++++----------
M
src/kernels/RayleighPlesset.h
 | 
7
+++++--
M
src/kernels/RayleighPlesset_HBGL.h
 | 
9
++++++---

6 files changed, 51 insertions(+), 42 deletions(-)
diff --git a/src/kernels/KMClusterPositions_D.h b/src/kernels/KMClusterPositions_D.h
@@ -58,7 +58,7 @@ public:
                 Bubble& b = data.coords[i];
                 infile >> b.pos[0] >> b.pos[1] >> b.pos[2];
                 infile >> data.R0[i];
-                infile >> data.Rdot0[i]; 
+                infile >> data.Rdot0[i];
 
                 State8& IC = U[i];
                 IC[0] = data.R0[i];
@@ -126,7 +126,7 @@ public:
                       A(ip,_N+3*i+(rr+2)%3) = static_cast<Real>(0.0);
                       A(ip,j) = static_cast<Real>(0.0);
                     }
-                }  
+                }
                 else
                 {
                     // compute distance between bubbles
@@ -233,17 +233,18 @@ public:
             }
 
             // compute bi
-            const Real pstat = bd.p0 - bd.pv;
             const Real Rdc = Ui[1]*clInv;
-            const Real pbs = (pstat + static_cast<Real>(2.0)*bd.sigma/bd.R0[i])*std::pow(bd.R0[i]/Ui[0], 3.0*bd.gamma);
+            const Real pB = (bd.p0 - bd.pv + static_cast<Real>(2.0)*bd.sigma/bd.R0[i])*std::pow(bd.R0[i]/Ui[0], 3.0*bd.gamma); // bubble pressure
 
             b(i)  = (static_cast<Real>(0.5)*Rdc - static_cast<Real>(1.5))*Ui[1]*Ui[1];
-            b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t + Ui[0]*clInv));
-            // b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t));
-            b(i) += bd.rInv*(pbs - static_cast<Real>(2.0)*bd.sigma/Ui[0]) - static_cast<Real>(4.0)*bd.nuL*Ui[1]/Ui[0];
-            b(i) += bd.rInv*static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc*pbs;
-            b(i) -= bd.rInv*clInv*Ui[0]*bd.dpext(t + Ui[0]*clInv);
-            // b(i) -= bd.rInv*Rdc*bd.dpext(t);
+            b(i) += bd.rInv*(static_cast<Real>(1.0) + static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc)*pB;
+            b(i) -= bd.rInv*static_cast<Real>(2.0)*bd.sigma/Ui[0];
+            b(i) -= static_cast<Real>(4.0)*bd.nuL*Ui[1]/Ui[0];
+            b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t + Ui[0]*clInv));
+            // b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t));
+            b(i) -= bd.rInv*Ui[0]*clInv*bd.dpext(t + Ui[0]*clInv);
+            // b(i) -= bd.rInv*Ui[0]*clInv*bd.dpext(t);
+
             // coupling terms
             b(i) -= static_cast<Real>(2.0)*bnbr;
             // translation term
diff --git a/src/kernels/KMCluster_FC.h b/src/kernels/KMCluster_FC.h
@@ -110,7 +110,7 @@ public:
                 {
                     const Real UDij = Uj[0]/bd.Dij[j + i*_N];
                     A(i,j) = (static_cast<Real>(1) + Ui[1]*clInv)*Uj[0]*UDij + static_cast<Real>(1.0)*Uj[1]*UDij*Ui[0]*clInv;
-                    // extension of Matrix: keeping all terms would yield a factor of 6 for the second term, 
+                    // extension of Matrix: keeping all terms would yield a factor of 6 for the second term,
                     // following simplifications of Fuster and Colonius results in a factor of 1
                     bnbr += Uj[1]*Uj[1]*UDij;
                     bnbr2 += Uj[1]*Uj[1]*Uj[1]/bd.Dij[j + i*_N];
@@ -118,17 +118,17 @@ public:
             }
 
             // compute bi
-            const Real pstat = bd.p0 - bd.pv;
             const Real Rdc = Ui[1]*clInv;
-            const Real pbs = (pstat + static_cast<Real>(2.0)*bd.sigma/bd.R0[i])*std::pow(bd.R0[i]/Ui[0], 3.0*bd.gamma);
+            const Real pB = (bd.p0 - bd.pv + static_cast<Real>(2.0)*bd.sigma/bd.R0[i])*std::pow(bd.R0[i]/Ui[0], 3.0*bd.gamma); // bubble pressure
 
             b(i)  = (static_cast<Real>(0.5)*Rdc - static_cast<Real>(1.5))*Ui[1]*Ui[1];
-            b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t + Ui[0]*clInv));
-            // b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t));
-            b(i) += bd.rInv*(pbs - static_cast<Real>(2.0)*bd.sigma/Ui[0]) - static_cast<Real>(4.0)*bd.nuL*Ui[1]/Ui[0];
-            b(i) += bd.rInv*static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc*pbs;
-            b(i) -= bd.rInv*Rdc*bd.dpext(t + Ui[0]*clInv);
-            // b(i) -= bd.rInv*Rdc*bd.dpext(t);
+            b(i) += bd.rInv*(static_cast<Real>(1.0) + static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc)*pB;
+            b(i) -= bd.rInv*static_cast<Real>(2.0)*bd.sigma/Ui[0];
+            b(i) -= static_cast<Real>(4.0)*bd.nuL*Ui[1]/Ui[0];
+            b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t + Ui[0]*clInv));
+            // b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t));
+            b(i) -= bd.rInv*Ui[0]*clInv*bd.dpext(t + Ui[0]*clInv);
+            // b(i) -= bd.rInv*Ui[0]*clInv*bd.dpext(t);
             b(i) -= static_cast<Real>(2.0)*(static_cast<Real>(1.0) + Rdc)*bnbr;
             b(i) -= static_cast<Real>(2.0)*Ui[0]*clInv*bnbr2;
 
diff --git a/src/kernels/KMCluster_TY.h b/src/kernels/KMCluster_TY.h
@@ -114,17 +114,17 @@ public:
             }
 
             // compute bi
-            const Real pstat = bd.p0 - bd.pv;
             const Real Rdc = Ui[1]*clInv;
-            const Real pbs = (pstat + static_cast<Real>(2.0)*bd.sigma/bd.R0[i])*std::pow(bd.R0[i]/Ui[0], 3.0*bd.gamma);
+            const Real pB = (bd.p0 - bd.pv + static_cast<Real>(2.0)*bd.sigma/bd.R0[i])*std::pow(bd.R0[i]/Ui[0], 3.0*bd.gamma); // bubble pressure
 
             b(i)  = (static_cast<Real>(0.5)*Rdc - static_cast<Real>(1.5))*Ui[1]*Ui[1];
-            b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t + Ui[0]*clInv));
-            // b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t));
-            b(i) += bd.rInv*(pbs - static_cast<Real>(2.0)*bd.sigma/Ui[0]) - static_cast<Real>(4.0)*bd.nuL*Ui[1]/Ui[0];
-            b(i) += bd.rInv*static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc*pbs;
-            b(i) -= bd.rInv*clInv*Ui[0]*bd.dpext(t + Ui[0]*clInv);
-            // b(i) -= bd.rInv*Rdc*bd.dpext(t);
+            b(i) += bd.rInv*(static_cast<Real>(1.0) + static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc)*pB;
+            b(i) -= bd.rInv*static_cast<Real>(2.0)*bd.sigma/Ui[0];
+            b(i) -= static_cast<Real>(4.0)*bd.nuL*Ui[1]/Ui[0];
+            b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t + Ui[0]*clInv));
+            // b(i) -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t));
+            b(i) -= bd.rInv*Ui[0]*clInv*bd.dpext(t + Ui[0]*clInv);
+            // b(i) -= bd.rInv*Ui[0]*clInv*bd.dpext(t);
             b(i) -= static_cast<Real>(2.0)*bnbr;
 
             assert(!isnan(b(i)));
diff --git a/src/kernels/KellerMiksis.h b/src/kernels/KellerMiksis.h
@@ -46,17 +46,19 @@ public:
 
         r[0] = u[1];
 
-        const Real pstat = bd.p0 - bd.pv;
-        const Real Rdc = u[1]/bd.cl;
-        const Real pbs = (pstat + static_cast<Real>(2.0)*bd.sigma/bd.R0[0])*std::pow(bd.R0[0]/u[0], 3.0*bd.gamma);
+        const Real clInv = 1.0/bd.cl;
+        const Real Rdc = u[1]*clInv;
+        const Real pB = (bd.p0 - bd.pv + static_cast<Real>(2.0)*bd.sigma/bd.R0[0])*std::pow(bd.R0[0]/u[0], 3.0*bd.gamma); // bubble pressure
+
         r[1]  = (static_cast<Real>(0.5)*Rdc - static_cast<Real>(1.5))*u[1]*u[1];
-        r[1] -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t + u[0]/bd.cl));
-        // r[1] -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(pstat + bd.pext(t));
-        r[1] += bd.rInv*(pbs - static_cast<Real>(2.0)*bd.sigma/u[0]) - static_cast<Real>(4.0)*bd.nuL*u[1]/u[0];
-        r[1] += bd.rInv*static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc*pbs;
-        r[1] -= bd.rInv*Rdc*bd.dpext(t + u[0]/bd.cl);
-        // r[1] -= bd.rInv*Rdc*bd.dpext(t);
-        r[1] /= (static_cast<Real>(1.0 - Rdc)*u[0] + static_cast<Real>(4.0)*bd.nuL/bd.cl);
+        r[1] += bd.rInv*(static_cast<Real>(1.0) + static_cast<Real>(1.0 - 3.0*bd.gamma)*Rdc)*pB;
+        r[1] -= bd.rInv*static_cast<Real>(2.0)*bd.sigma/u[0];
+        r[1] -= static_cast<Real>(4.0)*bd.nuL*u[1]/u[0];
+        r[1] -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t + u[0]*clInv));
+        // r[1] -= bd.rInv*(static_cast<Real>(1.0) + Rdc)*(bd.p0 - bd.pv + bd.pext(t));
+        r[1] -= bd.rInv*u[0]*clInv*bd.dpext(t + u[0]*clInv);
+        // r[1] -= bd.rInv*u[0]*clInv*bd.dpext(t);
+        r[1] /= (static_cast<Real>(1.0 - Rdc)*u[0] + static_cast<Real>(4.0)*bd.nuL*clInv);
 
         assert(!isnan(r[0]));
         assert(!isnan(r[1]));
diff --git a/src/kernels/RayleighPlesset.h b/src/kernels/RayleighPlesset.h
@@ -45,7 +45,7 @@ public:
 
         r[0] = u[1];
 
-        r[1]  = bd.rInv*(bd.pv - bd.pext(t));
+        r[1]  = bd.rInv*(bd.pv - bd.p0 - bd.pext(t));
         r[1] += bd.rInv*(bd.p0 - bd.pv + static_cast<Real>(2.0)*bd.sigma/bd.R0[0])*std::pow(bd.R0[0]/u[0], 3.0*bd.gamma);
         r[1] -= static_cast<Real>(1.5)*u[1]*u[1];
         r[1] -= static_cast<Real>(4.0)*bd.nuL*u[1]/u[0];
@@ -58,7 +58,10 @@ public:
 
     virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
     {
-        m_dumper->write(step, t, dt, U.data(), Tinput::DataType::SIZE*U.size());
+        const BubbleData& bd = *(BubbleData const* const)data;
+        const Real odat[3] = {U[0][0], U[0][1], bd.pext(t)};
+        m_dumper->write(step, t, dt, &odat[0], 3);
+        // m_dumper->write(step, t, dt, U.data(), Tinput::DataType::SIZE*U.size());
     }
 };
 
diff --git a/src/kernels/RayleighPlesset_HBGL.h b/src/kernels/RayleighPlesset_HBGL.h
@@ -48,8 +48,8 @@ public:
 
         const Real dR03 = (bd.R0[0]*bd.R0[0]*bd.R0[0] - bd.h*bd.h*bd.h);
         const Real dR3  = (u[0]*u[0]*u[0] - bd.h*bd.h*bd.h);
-        r[1]  = bd.rInv*(bd.p0 + static_cast<Real>(2.0)*bd.sigma/bd.R0[0])*std::pow(dR03/dR3,bd.gamma)*(static_cast<Real>(1.0) - static_cast<Real>(3.0)*bd.gamma*u[0]*u[0]*u[0]*u[1]/(bd.cl*dR3));
-        r[1] -= bd.rInv*bd.pext(t);
+        r[1]  = bd.rInv*(bd.pv - bd.p0 - bd.pext(t));
+        r[1] += bd.rInv*(bd.p0 - bd.pv + static_cast<Real>(2.0)*bd.sigma/bd.R0[0])*std::pow(dR03/dR3,bd.gamma)*(static_cast<Real>(1.0) - static_cast<Real>(3.0)*bd.gamma*u[0]*u[0]*u[0]*u[1]/(bd.cl*dR3));
         r[1] -= static_cast<Real>(1.5)*u[1]*u[1];
         r[1] -= static_cast<Real>(4.0)*bd.nuL*u[1]/u[0];
         r[1] -= static_cast<Real>(2.0)*bd.sigma*bd.rInv/u[0];
@@ -61,7 +61,10 @@ public:
 
     virtual void write(const size_t step, const Real t, const Real dt, const Tinput& U, const void * const data=nullptr)
     {
-        m_dumper->write(step, t, dt, U.data(), Tinput::DataType::SIZE*U.size());
+        const BubbleData& bd = *(BubbleData const* const)data;
+        const Real odat[3] = {U[0][0], U[0][1], bd.pext(t)};
+        m_dumper->write(step, t, dt, &odat[0], 3);
+        // m_dumper->write(step, t, dt, U.data(), Tinput::DataType::SIZE*U.size());
     }
 };