MAST
piston_theory_boundary_condition.cpp
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2019 Manav Bhatia
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 
21 // MAST includes
23 
24 
25 
27 PistonTheoryPressure(unsigned int order,
30  const MAST::FieldFunction<Real>& rho,
31  const MAST::FieldFunction<Real>& gamma,
32  const MAST::FieldFunction<Real>& dwdx,
33  const MAST::FieldFunction<Real>& dwdt):
34 MAST::FieldFunction<Real>("pressure"),
35 _order(order),
36 _V_inf(V),
37 _M_inf(M),
38 _rho_inf(rho),
39 _gamma(gamma),
40 _dwdx(dwdx),
41 _dwdt(dwdt) {
42 
43  _functions.insert(&_V_inf);
44  _functions.insert(&_M_inf);
45  _functions.insert(&_rho_inf);
46  _functions.insert(&_gamma);
47  _functions.insert(&_dwdx);
48  _functions.insert(&_dwdt);
49 }
50 
51 
52 
53 
55 
56 
57 
58 void
60  const Real t,
61  Real& m) const {
62 
63  Real V, M, rho, gamma, dwdx, dwdt;
64  _V_inf (p, t, V);
65  _M_inf (p, t, M);
66  _rho_inf (p, t, rho);
67  _gamma (p, t, gamma);
68  _dwdx (p, t, dwdx);
69  _dwdt (p, t, dwdt);
70 
71 
72  Real
73  v0 = dwdx + (M*M-2)/(M*M-1)*dwdt/V,
74  v1 = V*V*dwdx + (M*M-2)/(M*M-1)*V*dwdt,
75  v2 = pow(V*dwdx + (M*M-2)/(M*M-1)*dwdt, 2),
76  v3 = v0*v2;
77 
78  m = v1; // linear
79  switch (_order) {
80 
81  case 3:
82  m += (gamma+1)/12.*M*M*v3;
83  case 2:
84  m += (gamma+1)/4.*M*v2;
85  }
86 
87  m *= rho/sqrt(M*M-1.);
88 }
89 
90 
91 
92 
93 void
95  const libMesh::Point& p,
96  const Real t,
97  Real& m) const {
98 
99  Real V, M, rho, gamma, dwdx, dwdt, dV, dM, drho, dgamma, ddwdx, ddwdt;
100  _V_inf (p, t, V); _V_inf.derivative (f, p, t, dV);
101  _M_inf (p, t, M); _M_inf.derivative (f, p, t, dM);
102  _rho_inf (p, t, rho); _rho_inf.derivative (f, p, t, drho);
103  _gamma (p, t, gamma); _gamma.derivative (f, p, t, dgamma);
104  _dwdx (p, t, dwdx); _dwdx.derivative (f, p, t, ddwdx);
105  _dwdt (p, t, dwdt); _dwdt.derivative (f, p, t, ddwdt);
106 
107  Real
108  Mf = (M*M-2)/(M*M-1),
109  dMf = 2*M*dM/(M*M-1) - (M*M-2)/pow(M*M-1,2)*2*M*dM,
110 
111  v0 = dwdx + Mf*dwdt/V,
112  v1 = V*V*dwdx + Mf*V*dwdt,
113  v2 = pow(V*dwdx + Mf*dwdt, 2),
114  v3 = v0*v2,
115 
116  v0dw = ddwdx,
117  v1dw = V*V*ddwdx,
118  v2dw = 2.*(V*dwdx + Mf*dwdt)*V*ddwdx,
119  v3dw = v0*v2dw+v0dw*v2,
120 
121  v0dwt = Mf/V*ddwdt,
122  v1dwt = Mf*V*ddwdt,
123  v2dwt = 2.*(V*dwdx + Mf*dwdt)*Mf*ddwdt,
124  v3dwt = v0dwt*v2+v0*v2dwt,
125 
126  dv0dV = -Mf*dwdt/(V*V)*dV,
127  dv1dV = (2.*V*dwdx + Mf*dwdt)*dV,
128  dv2dV = 2.*(V*dwdx + Mf*dwdt)*dwdx*dV,
129  dv3dV = dv0dV*v2+v0*dv2dV,
130 
131  dv0dM = dMf*dwdt/V,
132  dv1dM = dMf*V*dwdt,
133  dv2dM = 2.*(V*dwdx + Mf*dwdt)*dMf*dwdt,
134  dv3dM = dv0dM*v2+v0*dv2dM,
135 
136  f1 = rho/sqrt(M*M-1.),
137  df1drho = drho/sqrt(M*M-1.),
138  df1dM = -.5*rho/pow(M*M-1., 1.5)*2*M*dM;
139 
140 
141  //
142  // linear part: f1 * v1,
143  // where f1 = rho/sqrt(M*M-1.) and
144  // v1 = V*V*dwdx + Mf*V*dwdt,
145  //
146  m = f1*(dv1dV + dv1dM + v1dw + v1dwt) + (df1drho + df1dM) * v1; // linear
147 
148  switch (_order) {
149 
150  case 3:
151  //
152  // cubic part: f1 * (gamma+1)/12.*M*M*v3
153  //
154  m +=
155  (df1drho + df1dM) *M*M*v3 * (gamma+1)/12.+
156  f1*M*(2*dM*v3 + M*(dv3dV + dv3dM + v3dw + v3dwt)) * (gamma+1)/12. +
157  f1 * M*M*v3/12.*dgamma;
158  case 2:
159  //
160  // quadratic part: f1 * (gamma+1)/4.*M*v2
161  //
162  m +=
163  (df1drho + df1dM) * M*v2 * (gamma+1)/4.+
164  f1*(dM*v2+M*(dv2dV + dv2dM + v2dw + v2dwt)) * (gamma+1)/4. +
165  f1*M*v2/4.*dgamma;
166  }
167 }
168 
169 
170 
171 
172 
175  const MAST::FieldFunction<Real>& V,
176  const MAST::FieldFunction<Real>& M,
177  const MAST::FieldFunction<Real>& rho,
178  const MAST::FieldFunction<Real>& gamma,
179  const MAST::FieldFunction<Real>& dwdx,
180  const MAST::FieldFunction<Real>& dwdt):
181 MAST::FieldFunction<Real>("pressure"),
182 _order(order),
183 _V_inf(V),
184 _M_inf(M),
185 _rho_inf(rho),
186 _gamma(gamma),
187 _dwdx(dwdx),
188 _dwdt(dwdt) {
189 
190  _functions.insert(&_V_inf);
191  _functions.insert(&_M_inf);
192  _functions.insert(&_rho_inf);
193  _functions.insert(&_gamma);
194  _functions.insert(&_dwdx);
195  _functions.insert(&_dwdt);
196 }
197 
198 
199 
200 
201 
203 
204 
205 
206 void
208  const Real t,
209  Real& m) const {
210 
211  Real V, M, rho, gamma, dwdx, dwdt;
212  _V_inf (p, t, V);
213  _M_inf (p, t, M);
214  _rho_inf (p, t, rho);
215  _gamma (p, t, gamma);
216  _dwdx (p, t, dwdx);
217  _dwdt (p, t, dwdt);
218 
219 
220  Real
221  Mf = (M*M-2)/(M*M-1),
222 
223  v0 = dwdx + (M*M-2)/(M*M-1)*dwdt/V,
224  v2 = pow(V*dwdx + (M*M-2)/(M*M-1)*dwdt, 2),
225 
226  v0dw = 1.,
227  v1dw = V*V,
228  v2dw = 2.*(V*dwdx + Mf*dwdt)*V,
229  v3dw = v0*v2dw+v0dw*v2,
230 
231  f1 = rho/sqrt(M*M-1.);
232 
233 
234  //
235  // linear part: f1 * v1,
236  // where f1 = rho/sqrt(M*M-1.) and
237  // v1 = V*V*dwdx + Mf*V*dwdt,
238  //
239  m = f1*v1dw; // linear
240 
241  switch (_order) {
242 
243  case 3:
244  //
245  // cubic part: f1 * (gamma+1)/12.*M*M*v3
246  //
247  m += f1*M*M*v3dw * (gamma+1)/12.;
248  case 2:
249  //
250  // quadratic part: f1 * (gamma+1)/4.*M*v2
251  //
252  m += f1*M*v2dw * (gamma+1)/4.;
253  }
254 }
255 
256 
257 
258 
259 void
261  const libMesh::Point& p,
262  const Real t,
263  Real& m) const {
264 
265  Real V, M, rho, gamma, dwdx, dwdt, dV, dM, drho, dgamma;
266  _V_inf (p, t, V); _V_inf.derivative (f, p, t, dV);
267  _M_inf (p, t, M); _M_inf.derivative (f, p, t, dM);
268  _rho_inf (p, t, rho); _rho_inf.derivative (f, p, t, drho);
269  _gamma (p, t, gamma); _gamma.derivative (f, p, t, dgamma);
270  _dwdx (p, t, dwdx);
271  _dwdt (p, t, dwdt);
272 
273  Real
274  Mf = (M*M-2)/(M*M-1),
275  dMfdM = 2*M*dM/(M*M-1) - (M*M-2)/pow(M*M-1,2)*2*M*dM,
276 
277  v0 = dwdx + Mf*dwdt/V,
278  v2 = pow(V*dwdx + Mf*dwdt, 2),
279 
280  v0dw = 1.,
281  v1dw = V*V,
282  v2dw = 2.*(V*dwdx + Mf*dwdt)*V,
283  v3dw = v0*v2dw+v0dw*v2,
284 
285  dv0dV = -Mf*dwdt/(V*V)*dV,
286  dv2dV = 2.*(V*dwdx + Mf*dwdt)*dwdx*dV,
287 
288  dv0dM = dMfdM*dwdt/V,
289  dv2dM = 2.*(V*dwdx + Mf*dwdt)*dMfdM*dwdt,
290 
291  v0dwdV = 0.,
292  v1dwdV = 2.*V*dV,
293  v2dwdV = 2.*(2.*V*dV*dwdx + Mf*dV*dwdt),
294  v3dwdV = dv0dV*v2dw+v0*v2dwdV + v0dwdV*v2+v0dw*dv2dV,
295 
296  v0dwdM = 0.,
297  v1dwdM = 0.,
298  v2dwdM = 2.*dMfdM*dwdt*V,
299  v3dwdM = dv0dM*v2dw+v0*v2dwdM + v0dwdM*v2+v0dw*dv2dM,
300 
301  f1 = rho/sqrt(M*M-1.),
302  df1drho = drho/sqrt(M*M-1.),
303  df1dM = -.5*rho/pow(M*M-1., 1.5)*2*M*dM;
304 
305 
306  //
307  // linear part: d/dp (f1 * dv1/dw),
308  // where f1 = rho/sqrt(M*M-1.) and
309  // dv1dw = V*V*dwdx,
310  //
311  m = f1*(v1dwdV + v1dwdM) + (df1drho + df1dM) * v1dw; // linear
312 
313  switch (_order) {
314 
315  case 3:
316  //
317  // cubic part: d/dp (f1 * (gamma+1)/12.*M*M* dv3/dw)
318  //
319  m +=
320  (df1drho + df1dM) *M*M*v3dw * (gamma+1)/12.+
321  f1*M* (2*dM*v3dw + M*(v3dwdV + v3dwdM)) * (gamma+1)/12. +
322  f1 * M*M*v3dw/12.*dgamma;
323  case 2:
324  //
325  // quadratic part: d/dp (f1 * (gamma+1)/4.*M* dv2/dw)
326  //
327  m +=
328  (df1drho + df1dM) * M*v2dw * (gamma+1)/4.+
329  f1*(dM*v2dw + M*(v2dwdV + v2dwdM)) * (gamma+1)/4. +
330  f1*M*v2dw/4.*dgamma;
331  }
332 }
333 
334 
335 
336 
339  const MAST::FieldFunction<Real>& V,
340  const MAST::FieldFunction<Real>& M,
341  const MAST::FieldFunction<Real>& rho,
342  const MAST::FieldFunction<Real>& gamma,
343  const MAST::FieldFunction<Real>& dwdx,
344  const MAST::FieldFunction<Real>& dwdt):
345 MAST::FieldFunction<Real>("pressure"),
346 _order(order),
347 _V_inf(V),
348 _M_inf(M),
349 _rho_inf(rho),
350 _gamma(gamma),
351 _dwdx(dwdx),
352 _dwdt(dwdt) {
353 
354  _functions.insert(&_V_inf);
355  _functions.insert(&_M_inf);
356  _functions.insert(&_rho_inf);
357  _functions.insert(&_gamma);
358  _functions.insert(&_dwdx);
359  _functions.insert(&_dwdt);
360 }
361 
362 
363 
364 
366 
367 
368 
369 void
371  const Real t,
372  Real& m) const {
373 
374  Real V, M, rho, gamma, dwdx, dwdt;
375  _V_inf (p, t, V);
376  _M_inf (p, t, M);
377  _rho_inf (p, t, rho);
378  _gamma (p, t, gamma);
379  _dwdx (p, t, dwdx);
380  _dwdt (p, t, dwdt);
381 
382 
383  Real
384  Mf = (M*M-2)/(M*M-1),
385  f1 = rho/sqrt(M*M-1.),
386 
387  v0 = dwdx + (M*M-2)/(M*M-1)*dwdt/V,
388  v2 = pow(V*dwdx + (M*M-2)/(M*M-1)*dwdt, 2),
389 
390  v0dwt = Mf/V,
391  v1dwt = Mf*V,
392  v2dwt = 2.*(V*dwdx + Mf*dwdt)*Mf,
393  v3dwt = v0dwt*v2+v0*v2dwt;
394 
395  //
396  // linear part: f1 * dv1/dwt,
397  // where f1 = rho/sqrt(M*M-1.) and
398  // v1 = V*V*dwdx + Mf*V*dwdt,
399  //
400  m = f1*v1dwt; // linear
401 
402  switch (_order) {
403 
404  case 3:
405  //
406  // cubic part: f1 * (gamma+1)/12.*M*M*dv3/dwt
407  //
408  m += f1*M*M*v3dwt * (gamma+1)/12.;
409  case 2:
410  //
411  // quadratic part: f1 * (gamma+1)/4.*M*dv2/dwt
412  //
413  m += f1*M*v2dwt * (gamma+1)/4.;
414  }
415 }
416 
417 
418 
419 
420 void
422  const libMesh::Point& p,
423  const Real t,
424  Real& m) const {
425 
426  Real V, M, rho, gamma, dwdx, dwdt, dV, dM, drho, dgamma;
427  _V_inf (p, t, V); _V_inf.derivative (f, p, t, dV);
428  _M_inf (p, t, M); _M_inf.derivative (f, p, t, dM);
429  _rho_inf (p, t, rho); _rho_inf.derivative (f, p, t, drho);
430  _gamma (p, t, gamma); _gamma.derivative (f, p, t, dgamma);
431  _dwdx (p, t, dwdx);
432  _dwdt (p, t, dwdt);
433 
434  Real
435  Mf = (M*M-2)/(M*M-1),
436  dMfdM = 2*M*dM/(M*M-1) - (M*M-2)/pow(M*M-1,2)*2*M*dM,
437 
438  v0 = dwdx + Mf*dwdt/V,
439  v2 = pow(V*dwdx + Mf*dwdt, 2),
440 
441  dv0dV = -Mf*dwdt/(V*V)*dV,
442  dv2dV = 2.*(V*dwdx + Mf*dwdt)*dwdx*dV,
443 
444  dv0dM = dMfdM*dwdt/V,
445  dv2dM = 2.*(V*dwdx + Mf*dwdt)*dMfdM*dwdt,
446 
447  v0dwt = Mf/V,
448  v1dwt = Mf*V,
449  v2dwt = 2.*(V*dwdx + Mf*dwdt)*Mf,
450  v3dwt = v0dwt*v2+v0*v2dwt,
451 
452  v0dwtdV = -Mf/(V*V)*dV,
453  v1dwtdV = Mf*dV,
454  v2dwtdV = 2.*dV*dwdx *Mf,
455  v3dwtdV = v0dwtdV*v2+v0dwt*dv2dV + dv0dV*v2dwt+v0*v2dwtdV,
456 
457  v0dwtdM = dMfdM/V,
458  v1dwtdM = dMfdM*V,
459  v2dwtdM = 2.*(V*dMfdM*dwdx + 2*Mf*dMfdM*dwdt),
460  v3dwtdM = v0dwtdM*v2+v0dwt*dv2dM + dv0dM*v2dwt+v0*v2dwtdM,
461 
462  f1 = rho/sqrt(M*M-1.),
463  df1drho = drho/sqrt(M*M-1.),
464  df1dM = -.5*rho/pow(M*M-1., 1.5)*2*M*dM;
465 
466 
467  //
468  // linear part: d/dp (f1 * dv1/dwt),
469  // where f1 = rho/sqrt(M*M-1.) and
470  // v1 = V*V*dwdx + Mf*V*dwdt,
471  //
472  m = f1*(v1dwtdV + v1dwtdM) + (df1drho + df1dM) * v1dwt; // linear
473 
474  switch (_order) {
475 
476  case 3:
477  //
478  // cubic part: d/dp (f1 * (gamma+1)/12.*M*M* dv3/dwt)
479  //
480  m +=
481  (df1drho + df1dM) *M*M*v3dwt * (gamma+1)/12.+
482  f1*M*(2*dM*v3dwt + M*(v3dwtdV + v3dwtdM)) * (gamma+1)/12. +
483  f1 * M*M*v3dwt/12.*dgamma;
484  case 2:
485  //
486  // quadratic part: d/dp (f1 * (gamma+1)/4.*M* dv2/dwt)
487  //
488  m +=
489  (df1drho + df1dM) * M*v2dwt * (gamma+1)/4. +
490  f1*(dM*v2dwt + M*(v2dwtdV + v2dwtdM)) * (gamma+1)/4. +
491  f1*M*v2dwt/4.*dgamma;
492  }
493 }
494 
495 
496 
497 
499 PistonTheoryBoundaryCondition(unsigned int order,
500  const RealVectorX& vel_vec):
502 _order(order),
503 _vel_vec(vel_vec) {
504 
505  // add the parameters
506 
507 
508  // scale the velocity vector to unit length
509  _vel_vec.normalize();
510 }
511 
512 
514 
515 }
516 
517 
518 
519 unsigned int
521  return _order;
522 }
523 
524 
525 
526 const RealVectorX&
528  return _vel_vec;
529 }
530 
531 
532 
533 
534 
535 std::unique_ptr<MAST::FieldFunction<Real> >
538  const MAST::FieldFunction<Real>& dwdt) const {
539 
541  *rval = new MAST::PistonTheoryPressure
542  (_order,
543  this->get<MAST::FieldFunction<Real> >("V"),
544  this->get<MAST::FieldFunction<Real> >("mach"),
545  this->get<MAST::FieldFunction<Real> >("rho"),
546  this->get<MAST::FieldFunction<Real> >("gamma"),
547  dwdx,
548  dwdt);
549 
550  return std::unique_ptr<MAST::FieldFunction<Real> >(rval);
551 }
552 
553 
554 
555 
556 std::unique_ptr<MAST::FieldFunction<Real> >
559  const MAST::FieldFunction<Real>& dwdt) const {
560 
563  (_order,
564  this->get<MAST::FieldFunction<Real> >("V"),
565  this->get<MAST::FieldFunction<Real> >("mach"),
566  this->get<MAST::FieldFunction<Real> >("rho"),
567  this->get<MAST::FieldFunction<Real> >("gamma"),
568  dwdx,
569  dwdt);
570 
571  return std::unique_ptr<MAST::FieldFunction<Real> >(rval);
572 }
573 
574 
575 
576 
577 std::unique_ptr<MAST::FieldFunction<Real> >
580  const MAST::FieldFunction<Real>& dwdt) const {
581 
584  (_order,
585  this->get<MAST::FieldFunction<Real> >("V"),
586  this->get<MAST::FieldFunction<Real> >("mach"),
587  this->get<MAST::FieldFunction<Real> >("rho"),
588  this->get<MAST::FieldFunction<Real> >("gamma"),
589  dwdx,
590  dwdt);
591 
592  return std::unique_ptr<MAST::FieldFunction<Real> >(rval);
593 }
594 
595 
virtual void derivative(const MAST::FunctionBase &f, const libMesh::Point &p, const Real t, Real &m) const
calculates the value of the derivative of function with respect to the function f at the specified po...
const MAST::FieldFunction< Real > & _gamma
PistonTheoryPressureXdotDerivative(unsigned int order, const MAST::FieldFunction< Real > &V, const MAST::FieldFunction< Real > &M, const MAST::FieldFunction< Real > &rho, const MAST::FieldFunction< Real > &gamma, const MAST::FieldFunction< Real > &dwdx, const MAST::FieldFunction< Real > &dwdt)
virtual void operator()(const libMesh::Point &p, const Real t, Real &m) const
calculates the value of the function at the specified point, p, and time, t, and returns it in v...
virtual void operator()(const libMesh::Point &p, const Real t, Real &m) const
calculates the value of the function at the specified point, p, and time, t, and returns it in v...
const MAST::FieldFunction< Real > & _dwdt
virtual void derivative(const MAST::FunctionBase &f, const libMesh::Point &p, const Real t, Real &m) const
calculates the value of the derivative of function with respect to the function f at the specified po...
const MAST::FieldFunction< Real > & _M_inf
virtual void derivative(const MAST::FunctionBase &f, const libMesh::Point &p, const Real t, Real &m) const
calculates the value of the derivative of function with respect to the function f at the specified po...
std::set< const MAST::FunctionBase * > _functions
set of functions that this function depends on
PistonTheoryBoundaryCondition(unsigned int order, const RealVectorX &vel_vec)
Constructor for the Piston Theory boundary condition object.
libMesh::Real Real
virtual void derivative(const MAST::FunctionBase &f, ValType &v) const
calculates the value of the function derivative and returns it in v.
const MAST::FieldFunction< Real > & _V_inf
const MAST::FieldFunction< Real > & _dwdx
std::unique_ptr< MAST::FieldFunction< Real > > get_dpdx_function(const MAST::FieldFunction< Real > &dwdx, const MAST::FieldFunction< Real > &dwdt) const
Matrix< Real, Dynamic, 1 > RealVectorX
RealVectorX _vel_vec
Ambient flow velocity vector.
PistonTheoryPressureXDerivative(unsigned int order, const MAST::FieldFunction< Real > &V, const MAST::FieldFunction< Real > &M, const MAST::FieldFunction< Real > &rho, const MAST::FieldFunction< Real > &gamma, const MAST::FieldFunction< Real > &dwdx, const MAST::FieldFunction< Real > &dwdt)
std::unique_ptr< MAST::FieldFunction< Real > > get_dpdxdot_function(const MAST::FieldFunction< Real > &dwdx, const MAST::FieldFunction< Real > &dwdt) const
std::unique_ptr< MAST::FieldFunction< Real > > get_pressure_function(const MAST::FieldFunction< Real > &dwdx, const MAST::FieldFunction< Real > &dwdt) const
const MAST::FieldFunction< Real > & _rho_inf
PistonTheoryPressure(unsigned int order, const MAST::FieldFunction< Real > &V, const MAST::FieldFunction< Real > &M, const MAST::FieldFunction< Real > &rho, const MAST::FieldFunction< Real > &gamma, const MAST::FieldFunction< Real > &dwdx, const MAST::FieldFunction< Real > &dwdt)
unsigned int _order
Order of the boundary condition.
virtual void operator()(const libMesh::Point &p, const Real t, Real &m) const
calculates the value of the function at the specified point, p, and time, t, and returns it in v...