bubble-dynamics

commit 58431e995aa4ea27b52bf1a0258d11f077310ba7
parent 411bf7e617fbb5c7ccf46eb363362996db576f86
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date:   Fri,  9 Jun 2017 17:07:44 -0700

added Gaussian and sinusoidal pulse forcing

Diffstat:
M
src/BubbleBase.h
 | 
39
++++++++++++++++++++++++++-------------

1 file changed, 26 insertions(+), 13 deletions(-)
diff --git a/src/BubbleBase.h b/src/BubbleBase.h
@@ -52,26 +52,39 @@ struct BubbleData
         cl    = p("-c").asDouble(1481.0); // m/s
 
         // external forcing
-        pamp  = p("-pamp").asDouble(1.4); // external pressure amplification
-        omega = p("-omega").asDouble(2.0*M_PI*26500.0); // s^-1
-        p_pre = p("-p_pre").asDouble(1.0e5);
-        p_post= p("-p_post").asDouble(1.0e5);
-        t_12  = p("-t_12").asDouble(1.0e-6);
+        pamp      = p("-pamp").asDouble(1.4); // external pressure amplification
+        frequency = p("-frequency").asDouble(26500.0); // s^-1
+        omega     = 2.0*M_PI*frequency;
+        support   = p("-support").asDouble(1.0/omega);
+        p_pre     = p("-p_pre").asDouble(1.0e5);
+        p_post    = p("-p_post").asDouble(1.0e5);
+        t_0       = p("-t_0").asDouble(5.0/omega);
+        t_12      = p("-t_12").asDouble(1.0e-6);
         dt_smooth = p("-dt_smooth").asDouble(1.0e-9);
 
         // drivers
         if (p("-pext").asString("const") == "const")
             pext = [this] (const Real t) { return (this->pamp * this->p0); };
-        if (p("-pext").asString("const") == "sin")
+        else if (p("-pext").asString("const") == "sin")
+            pext = [this] (const Real t) { return (this->p0*this->pamp*std::sin(this->omega * t)); };
+        else if (p("-pext").asString("const") == "coshbgl")
             pext = [this] (const Real t) { return (-this->p0*this->pamp*std::sin(this->omega * t)); };
-        if (p("-pext").asString("const") == "coshbgl")
-            pext = [this] (const Real t) { return (this->p0 * (1.0 - this->pamp*std::cos(this->omega * (t - t_star)))); };
-        if (p("-pext").asString("const") == "cosramp")
+        else if (p("-pext").asString("const") == "cosramp")
             pext = [this] (const Real t)
             {
             const double phi = std::min(1.0,std::max(0.0,(t-this->t_12)/this->dt_smooth+0.5));
             return this->p_post + 0.5*(this->p_pre-this->p_post)*(1.0+std::cos(M_PI*phi));
             };
+        else if (p("-pext").asString("const") == "sinpulse")
+            pext = [this] (const Real t)
+            {
+            return (t>this->t_0 && t<=(this->t_0+1.0/this->frequency)) ? this->pamp*this->p0*std::sin(this->omega*(t-this->t_0)) : 0.0;
+            };
+        else if (p("-pext").asString("const") == "gaussian")
+            pext = [this] (const Real t)
+            {
+            return this->pamp*this->p0*std::exp(-0.5*(t-this->t_0)*(t-this->t_0)/(this->support*this->support));
+            };
 
         // pressure derivative
         dpext = [this] (const Real t) { return 0.0; };
@@ -83,7 +96,6 @@ struct BubbleData
         // helpers
         assert(rhoL > 0);
         rInv = 1.0/rhoL;
-        t_star = 1.0/omega * std::acos(1.0/pamp);
     }
 
     // Number of bubbles and their initial state
@@ -106,9 +118,11 @@ struct BubbleData
 
     // external forcing
     Real pamp;  // forcing pressure amplitude
-    Real omega; // external frequency
+    Real frequency;
+    Real omega;
+    Real support;
     Real p_pre, p_post;
-    Real t_12, dt_smooth;
+    Real t_0, t_12, dt_smooth;
 
     // driver functions
     std::function<Real(const Real)> pext;
@@ -118,7 +132,6 @@ struct BubbleData
     std::string filename;
     std::vector<Real> Dij;
     Real rInv;
-    Real t_star;
     bool bVTK;
 
     template <typename T>