Line data Source code
1 : #include "crpropa/advectionField/TimeDependentAdvectionField.h"
2 : #include "crpropa/Common.h"
3 :
4 :
5 : namespace crpropa {
6 :
7 1 : OneDimensionalTimeDependentShock::OneDimensionalTimeDependentShock(double v_sh, double v1, double v0, double l_sh){
8 :
9 1 : setShockSpeed(v_sh);
10 1 : setSpeeds(v1, v0);
11 1 : setShockWidth(l_sh);
12 1 : setShockPosition(0.);
13 1 : setShockTime(0.);
14 1 : }
15 :
16 3 : Vector3d OneDimensionalTimeDependentShock::getField(const Vector3d &position, const double &time) const {
17 :
18 3 : double A = 0.5 * ( v1 + v0 );
19 3 : double B = 0.5 * ( v1 - v0 );
20 :
21 3 : double x = position.x;
22 3 : double xsh = v_sh * (time-t_sh0) + x_sh0; // shock position
23 :
24 : Vector3d v(0.);
25 3 : v.x = A - B * tanh( (x - xsh) / l_sh );
26 :
27 3 : return v;
28 : }
29 :
30 1 : double OneDimensionalTimeDependentShock::getDivergence(const Vector3d &position, const double &time) const {
31 :
32 : double dvdx = 0.;
33 1 : double x = position.x;
34 :
35 1 : if (time < t_sh0){
36 : // no shock
37 : dvdx = 0;
38 : }
39 : else{
40 :
41 1 : double B = 0.5 * ( v1 - v0 );
42 1 : double xsh = v_sh * (time-t_sh0) + x_sh0; // shock position
43 1 : dvdx = - B / l_sh * ( 1 - tanh( (x - xsh) / l_sh ) * tanh( (x - xsh) / l_sh ) );
44 : }
45 :
46 1 : return dvdx;
47 : }
48 :
49 :
50 1 : void OneDimensionalTimeDependentShock::setShockSpeed(double v) {
51 1 : v_sh = v;
52 1 : }
53 :
54 1 : void OneDimensionalTimeDependentShock::setSpeeds(double u1, double u0) {
55 1 : v1 = u1;
56 1 : v0 = u0;
57 1 : }
58 :
59 1 : void OneDimensionalTimeDependentShock::setShockWidth(double w) {
60 1 : l_sh = w;
61 1 : }
62 :
63 1 : void OneDimensionalTimeDependentShock::setShockPosition(double x) {
64 1 : x_sh0 = x;
65 1 : }
66 :
67 1 : void OneDimensionalTimeDependentShock::setShockTime(double t) {
68 1 : t_sh0 = t;
69 1 : }
70 :
71 1 : double OneDimensionalTimeDependentShock::getVshock() const {
72 1 : return v_sh;
73 : }
74 :
75 1 : double OneDimensionalTimeDependentShock::getV1() const {
76 1 : return v1;
77 : }
78 :
79 1 : double OneDimensionalTimeDependentShock::getV0() const {
80 1 : return v0;
81 : }
82 :
83 1 : double OneDimensionalTimeDependentShock::getShockWidth() const {
84 1 : return l_sh;
85 : }
86 :
87 2 : double OneDimensionalTimeDependentShock::getShockPosition(double time) const {
88 2 : return x_sh0 + v_sh * (time - t_sh0);
89 : }
90 :
91 1 : double OneDimensionalTimeDependentShock::getShockTime() const {
92 1 : return t_sh0;
93 : }
94 :
95 : //----------------------------------------------------------------
96 :
97 1 : SedovTaylorBlastWave::SedovTaylorBlastWave(double E0, double rho0, double l_sh){
98 1 : setEnergy(E0);
99 1 : setDensity(rho0);
100 1 : setShockWidth(l_sh);
101 1 : }
102 :
103 3 : Vector3d SedovTaylorBlastWave::getField(const Vector3d &position, const double &time) const {
104 : double r = position.getR();
105 3 : Vector3d e_r = position.getUnitVector();
106 :
107 3 : if (time == 0)
108 : return 0. * e_r ;
109 0 : double A = rho0;
110 :
111 0 : double R = pow(E0 / A, 1./5.) * pow(time, 2. / 5. ); // position of shock front
112 0 : double vs = 2. / 5. * pow(E0 / A, 1. / 5.) * pow(time, -3. / 5.); // shock speed
113 :
114 0 : double xi = r / R; // dimensionless radius
115 0 : double V = 3 * (pow_integer<8>(xi) + 1) / (3 * pow_integer<8>(xi) + 5);
116 0 : double u = 0.5 * xi * vs * V * ( 1 - tanh( (xi - 1) * R / l_sh) );
117 :
118 : return u * e_r;
119 : }
120 :
121 2 : double SedovTaylorBlastWave::getDivergence(const Vector3d &position, const double &time) const {
122 2 : if (time == 0)
123 : return 0;
124 : double r = position.getR();
125 0 : double A = rho0;
126 :
127 0 : double R = pow(E0 / A, 1./5.) * pow(time, 2. / 5.);
128 0 : double vs = 2. / 5. * pow(E0 / A, 1. / 5.) * pow(time, -3. / 5.);
129 :
130 0 : double a = (-3 * pow_integer<3>(r) * (1 + pow_integer<8>(r/R)) * 1./pow_integer<2>(cosh((r - R)/l_sh)))
131 0 : / (l_sh * (5 + (3 * pow_integer<8>(r))/pow_integer<8>(R)))
132 0 : + (9 * pow_integer<2>(r) * (1 + pow_integer<8>(r/R)) * (1 - tanh((r - R)/l_sh)))
133 0 : / (5 + (3 * pow_integer<8>(r))/pow_integer<8>(R))
134 0 : - (72 *pow_integer<10>(r) * (1 + pow_integer<8>(r/R)) * (1 - tanh((r - R)/l_sh)))
135 0 : / (pow_integer<2>(5 + (3 * pow_integer<8>(r))/pow_integer<8>(R)) * pow_integer<8>(R))
136 0 : + (24 * pow_integer<10>(r) * (1 - tanh((r - R)/l_sh)))
137 0 : / ((5 + (3 * pow_integer<8>(r)) / pow_integer<8>(R)) * pow(R,8));
138 :
139 0 : double dudr = 0.5 * vs / pow_integer<2>(r) * 1./ R * a;
140 :
141 0 : return dudr;
142 : }
143 :
144 :
145 1 : void SedovTaylorBlastWave::setShockWidth(double w) {
146 1 : l_sh = w;
147 1 : }
148 :
149 1 : void SedovTaylorBlastWave::setDensity(double rho) {
150 1 : rho0 = rho;
151 1 : }
152 :
153 1 : void SedovTaylorBlastWave::setEnergy(double e) {
154 1 : E0 = e;
155 1 : }
156 :
157 1 : double SedovTaylorBlastWave::getShockRadius(double time) const{
158 1 : double A = rho0;
159 1 : return pow(E0 / A, 1./5.) * pow(time, 2./5.);
160 : }
161 :
162 1 : double SedovTaylorBlastWave::getShockSpeed(double time) const{
163 1 : double A = rho0;
164 1 : return 2. / 5. * pow(E0 / A, 1. / 5.) * pow(time, -3. / 5.);
165 : }
166 :
167 1 : double SedovTaylorBlastWave::getShockWidth() const {
168 1 : return l_sh;
169 : }
170 :
171 1 : double SedovTaylorBlastWave::getDensity() const {
172 1 : return rho0;
173 : }
174 :
175 1 : double SedovTaylorBlastWave::getEnergy() const {
176 1 : return E0;
177 : }
178 :
179 : //----------------------------------------------------------------
180 :
181 : } // namespace crpropa
|
