BubbleBase.h (5605B)
1 /* File: BubbleBase.h */ 2 /* Date: Fri Feb 12 12:39:08 2016 */ 3 /* Author: Fabian Wermelinger */ 4 /* Tag: Bubble base class */ 5 /* Copyright 2016 ETH Zurich. All Rights Reserved. */ 6 #ifndef BUBBLEBASE_H_ZVSZD3E8 7 #define BUBBLEBASE_H_ZVSZD3E8 8 9 #include <ODETB/ArgumentParser.h> 10 #include <ODETB/common.h> 11 12 #include <cmath> 13 #include <cstdio> 14 #include <functional> 15 #include <string> 16 #include <vector> 17 18 struct Bubble 19 { 20 double pos[3]; 21 double distance(const Bubble& b) const 22 { 23 const double dx = b.pos[0] - pos[0]; 24 const double dy = b.pos[1] - pos[1]; 25 const double dz = b.pos[2] - pos[2]; 26 return std::sqrt(dx*dx + dy*dy + dz*dz); 27 } 28 }; 29 30 struct BubbleData 31 { 32 BubbleData(ArgumentParser& p) 33 { 34 Nbubbles = p("-N").asInt(1); 35 R0.reserve(Nbubbles); 36 Rdot0.reserve(Nbubbles); 37 coords.reserve(Nbubbles); 38 Dij.reserve(Nbubbles*Nbubbles); 39 40 R0[0] = p("-R0").asDouble(4.0e-6); // m 41 Rdot0[0] = p("-Rdot0").asDouble(0.0); // m/s 42 coords[0].pos[0] = coords[0].pos[1] = coords[0].pos[2] = 0.0; 43 44 // gas 45 gamma = p("-gamma").asDouble(1.4); 46 47 // liquid 48 rhoL = p("-rho").asDouble(1000.0); // kg/m^3 49 nuL = p("-nu").asDouble(1.0e-6); // m^2/s 50 pv = p("-pv").asDouble(2400.0); // Pa 51 p0 = p("-p0").asDouble(1.0e5); // Pa 52 sigma = p("-sigma").asDouble(7.3e-2); // kg/s^2 53 h = p("-h").asDouble(4.514673e-7); 54 cl = p("-c").asDouble(1481.0); // m/s 55 56 // external forcing 57 pamp = p("-pamp").asDouble(1.4); // external pressure amplification 58 frequency = p("-frequency").asDouble(26500.0); // s^-1 59 omega = 2.0*M_PI*frequency; 60 support = p("-support").asDouble(1.0/omega); 61 p_pre = p("-p_pre").asDouble(1.0e5); 62 p_post = p("-p_post").asDouble(1.0e5); 63 t_0 = p("-t_0").asDouble(5.0/omega); 64 t_12 = p("-t_12").asDouble(1.0e-6); 65 dt_smooth = p("-dt_smooth").asDouble(1.0e-9); 66 67 // drivers 68 if (p("-pext").asString("const") == "const") 69 // implemented as (p0 + pext) therefore the -1. 70 pext = [this](const Real) { return ((this->pamp - 1) * this->p0); }; 71 else if (p("-pext").asString("const") == "sin") 72 pext = [this] (const Real t) { return (this->p0*this->pamp*std::sin(this->omega * t)); }; 73 else if (p("-pext").asString("const") == "coshbgl") 74 pext = [this] (const Real t) { return (-this->p0*this->pamp*std::sin(this->omega * t)); }; 75 else if (p("-pext").asString("const") == "cosramp") 76 pext = [this] (const Real t) 77 { 78 const double phi = std::min(1.0,std::max(0.0,(t-this->t_12)/this->dt_smooth+0.5)); 79 return this->p_post + 0.5*(this->p_pre-this->p_post)*(1.0+std::cos(M_PI*phi)); 80 }; 81 else if (p("-pext").asString("const") == "sinpulse") 82 pext = [this] (const Real t) 83 { 84 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; 85 }; 86 else if (p("-pext").asString("const") == "gaussian") 87 pext = [this] (const Real t) 88 { 89 return this->pamp*this->p0*std::exp(-0.5*(t-this->t_0)*(t-this->t_0)/(this->support*this->support)); 90 }; 91 92 // pressure derivative 93 dpext = [this](const Real) { return 0.0; }; 94 95 // other 96 filename = p("-bubbleFile").asString("none"); 97 bVTK = p("-vtk").asBool(false); 98 99 // helpers 100 assert(rhoL > 0); 101 rInv = 1.0/rhoL; 102 } 103 104 // Number of bubbles and their initial state 105 size_t Nbubbles; 106 std::vector<Real> R0; 107 std::vector<Real> Rdot0; 108 std::vector<Bubble> coords; 109 110 // Bubble properties 111 Real gamma; // ratio of specific heats for gas 112 113 // Liquid properties 114 Real rhoL; // density 115 Real nuL; // viscosity 116 Real pv; // vaporization pressure 117 Real p0; // far field pressure 118 Real sigma; // surface tension 119 Real h; // van der Waals hard-core radius 120 Real cl; // speed of sound 121 122 // external forcing 123 Real pamp; // forcing pressure amplitude 124 Real frequency; 125 Real omega; 126 Real support; 127 Real p_pre, p_post; 128 Real t_0, t_12, dt_smooth; 129 130 // driver functions 131 std::function<Real(const Real)> pext; 132 std::function<Real(const Real)> dpext; 133 134 // helpers 135 std::string filename; 136 std::vector<Real> Dij; 137 Real rInv; 138 bool bVTK; 139 140 template <typename T> 141 inline Real pBubble(const T& U, const size_t i=0) const 142 { 143 // return (p0 - pv + 2.0*sigma/R0[i])*std::pow(R0[i]/U[i][0], 3.0*gamma) - 2.0*sigma/U[i][0] - 4.0*nuL*rhoL/U[i][0]*U[i][1]; 144 return (p0 - pv + 2.0*sigma/R0[i])*std::pow(R0[i]/U[i][0], 3.0*gamma); 145 } 146 147 void print() const 148 { 149 printf("Bubble Data:\n"); 150 printf("\tNbubbles = %zu\n", Nbubbles); 151 printf("\tBubbles:\n"); 152 for (size_t i = 0; i < Nbubbles; ++i) 153 printf("\t\tBubble %zu: R0 = %e, Rdot0 = %e\n", i, R0[i], Rdot0[i]); 154 printf("\tgamma = %e\n", gamma); 155 printf("\trho = %e\n", rhoL); 156 printf("\tnu = %e\n", nuL); 157 printf("\tpv = %e\n", pv); 158 printf("\tp0 = %e\n", p0); 159 printf("\tsigma = %e\n", sigma); 160 printf("\th = %e\n", h); 161 printf("\tcl = %e\n", cl); 162 printf("\tpamp = %e\n", pamp); 163 printf("\tomega = %e\n", omega); 164 } 165 }; 166 167 #endif /* BUBBLEBASE_H_ZVSZD3E8 */