Nektar++
ProcessWSS.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessWSS.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Computes wss field.
32 //
33 ////////////////////////////////////////////////////////////////////////////////
34 
35 #include <iostream>
36 #include <string>
37 
38 #include "ProcessWSS.h"
39 
41 #include <MultiRegions/ExpList.h>
42 
43 using namespace std;
44 
45 namespace Nektar
46 {
47 namespace FieldUtils
48 {
49 
50 ModuleKey ProcessWSS::className = GetModuleFactory().RegisterCreatorFunction(
51  ModuleKey(eProcessModule, "wss"), ProcessWSS::create,
52  "Computes wall shear stress field.");
53 
54 ProcessWSS::ProcessWSS(FieldSharedPtr f) : ProcessBoundaryExtract(f)
55 {
56 }
57 
59 {
60 }
61 
62 void ProcessWSS::Process(po::variables_map &vm)
63 {
65 
66  int i, j;
67  int nfields = m_f->m_variables.size();
68  int expdim = m_f->m_graph->GetSpaceDimension();
69  m_spacedim = expdim + m_f->m_numHomogeneousDir;
70 
71  if (m_f->m_exp[0]->GetNumElmts() == 0)
72  {
73  return;
74  }
75 
76  if (m_spacedim == 1)
77  {
78  ASSERTL0(false, "Error: wss for a 1D problem cannot "
79  "be computed");
80  }
81 
82  // Declare arrays
83  int nshear = m_spacedim + 1;
84  int nstress = m_spacedim * m_spacedim;
85  int ngrad = m_spacedim * m_spacedim;
86 
87  Array<OneD, Array<OneD, NekDouble>> velocity(nfields);
89  Array<OneD, Array<OneD, NekDouble>> stress(nstress), fstress(nstress);
90  Array<OneD, Array<OneD, NekDouble>> fshear(nshear);
91 
94 
95  // will resuse nfields expansions to write shear components.
96  if (nshear > nfields)
97  {
98  m_f->m_exp.resize(nshear);
99  for (i = nfields; i < nshear; ++i)
100  {
101  m_f->m_exp[nfields + i] =
102  m_f->AppendExpList(m_f->m_numHomogeneousDir);
103  }
104  }
105 
106  // Create map of boundary ids for partitioned domains
108  m_f->m_exp[0]->GetGraph());
110  bcs.GetBoundaryRegions();
111  map<int, int> BndRegionMap;
112  int cnt = 0;
113  for (auto &breg_it : bregions)
114  {
115  BndRegionMap[breg_it.first] = cnt++;
116  }
117 
118  // Loop over boundaries to Write
119  for (int b = 0; b < m_f->m_bndRegionsToWrite.size(); ++b)
120  {
121  if (BndRegionMap.count(m_f->m_bndRegionsToWrite[b]) == 1)
122  {
123  int bnd = BndRegionMap[m_f->m_bndRegionsToWrite[b]];
124  // Get expansion list for boundary and for elements containing this
125  // bnd
126  for (i = 0; i < nshear; i++)
127  {
128  BndExp[i] = m_f->m_exp[i]->UpdateBndCondExpansion(bnd);
129  }
130  for (i = 0; i < nfields; i++)
131  {
132  m_f->m_exp[i]->GetBndElmtExpansion(bnd, BndElmtExp[i]);
133  }
134 
135  // Get number of points in expansions
136  int nqb = BndExp[0]->GetTotPoints();
137  int nqe = BndElmtExp[0]->GetTotPoints();
138 
139  // Initialise local arrays for the velocity gradients, and
140  // stress components size of total number of quadrature
141  // points for elements in this bnd
142  for (i = 0; i < ngrad; ++i)
143  {
144  grad[i] = Array<OneD, NekDouble>(nqe);
145  }
146 
147  for (i = 0; i < nstress; ++i)
148  {
149  stress[i] = Array<OneD, NekDouble>(nqe);
150  }
151 
152  Array<OneD, NekDouble> div(nqe, 0.0);
153 
154  // initialise arrays in the boundary
155  for (i = 0; i < nstress; ++i)
156  {
157  fstress[i] = Array<OneD, NekDouble>(nqb);
158  }
159 
160  for (i = 0; i < nshear; ++i)
161  {
162  fshear[i] = Array<OneD, NekDouble>(nqb, 0.0);
163  }
164 
165  // Extract Velocities
166  GetVelocity(BndElmtExp, velocity);
167 
168  // Extract viscosity coefficients
169  NekDouble lambda;
170  Array<OneD, NekDouble> mu(nqe, 0.0);
171  GetViscosity(BndElmtExp, mu, lambda);
172 
173  // Compute gradients
174  for (i = 0; i < m_spacedim; ++i)
175  {
176  if (m_spacedim == 2)
177  {
178  BndElmtExp[i]->PhysDeriv(velocity[i],
179  grad[i * m_spacedim + 0],
180  grad[i * m_spacedim + 1]);
181  }
182  else
183  {
184  BndElmtExp[i]->PhysDeriv(
185  velocity[i], grad[i * m_spacedim + 0],
186  grad[i * m_spacedim + 1], grad[i * m_spacedim + 2]);
187  }
188  // Add contribution to div(u)
189  Vmath::Vadd(nqe, grad[i * m_spacedim + i], 1, div, 1, div, 1);
190  }
191 
192  // Velocity divergence scaled by lambda * mu
193  Vmath::Smul(nqe, lambda, div, 1, div, 1);
194  Vmath::Vmul(nqe, mu, 1, div, 1, div, 1);
195 
196  // Compute stress component terms
197  // tau_ij = mu*(u_i,j + u_j,i) + mu*lambda*delta_ij*div(u)
198  for (i = 0; i < m_spacedim; ++i)
199  {
200  for (j = i; j < m_spacedim; ++j)
201  {
202  Vmath::Vadd(nqe, grad[i * m_spacedim + j], 1,
203  grad[j * m_spacedim + i], 1,
204  stress[i * m_spacedim + j], 1);
205 
206  Vmath::Vmul(nqe, mu, 1, stress[i * m_spacedim + j], 1,
207  stress[i * m_spacedim + j], 1);
208 
209  if (i == j)
210  {
211  // Add divergence term to diagonal
212  Vmath::Vadd(nqe, stress[i * m_spacedim + j], 1, div, 1,
213  stress[i * m_spacedim + j], 1);
214  }
215  else
216  {
217  // Copy to make symmetric
218  Vmath::Vcopy(nqe, stress[i * m_spacedim + j], 1,
219  stress[j * m_spacedim + i], 1);
220  }
221  }
222  }
223 
224  // Get boundary stress values.
225  for (j = 0; j < nstress; ++j)
226  {
227  m_f->m_exp[0]->ExtractElmtToBndPhys(bnd, stress[j], fstress[j]);
228  }
229 
230  // Get normals
232  m_f->m_exp[0]->GetBoundaryNormals(bnd, normals);
233  // Reverse normals, to get correct orientation for the body
234  for (i = 0; i < m_spacedim; ++i)
235  {
236  Vmath::Neg(nqb, normals[i], 1);
237  }
238 
239  // calculate wss, and update coeffs in the boundary expansion
240  // S = tau_ij * n_j
241  for (i = 0; i < m_spacedim; ++i)
242  {
243  for (j = 0; j < m_spacedim; ++j)
244  {
245  Vmath::Vvtvp(nqb, normals[j], 1,
246  fstress[i * m_spacedim + j], 1, fshear[i], 1,
247  fshear[i], 1);
248  }
249  }
250 
251  // T = S - (S.n)n
252  for (i = 0; i < m_spacedim; ++i)
253  {
254  Vmath::Vvtvp(nqb, normals[i], 1, fshear[i], 1,
255  fshear[nshear - 1], 1, fshear[nshear - 1], 1);
256  }
257  Vmath::Smul(nqb, -1.0, fshear[nshear - 1], 1, fshear[nshear - 1],
258  1);
259 
260  for (i = 0; i < m_spacedim; i++)
261  {
262  Vmath::Vvtvp(nqb, normals[i], 1, fshear[nshear - 1], 1,
263  fshear[i], 1, fshear[i], 1);
264  BndExp[i]->FwdTransLocalElmt(fshear[i],
265  BndExp[i]->UpdateCoeffs());
266  }
267 
268  // Tw
269  Vmath::Zero(nqb, fshear[nshear - 1], 1);
270  for (i = 0; i < m_spacedim; ++i)
271  {
272  Vmath::Vvtvp(nqb, fshear[i], 1, fshear[i], 1,
273  fshear[nshear - 1], 1, fshear[nshear - 1], 1);
274  }
275  Vmath::Vsqrt(nqb, fshear[nshear - 1], 1, fshear[nshear - 1], 1);
276  BndExp[nshear - 1]->FwdTransLocalElmt(
277  fshear[nshear - 1], BndExp[nshear - 1]->UpdateCoeffs());
278  }
279  }
280 
281  if (m_spacedim == 2)
282  {
283  m_f->m_variables[0] = "Shear_s";
284  m_f->m_variables[1] = "Shear_n";
285  m_f->m_variables[2] = "Shear_mag";
286  }
287  else
288  {
289  m_f->m_variables[0] = "Shear_s";
290  m_f->m_variables[1] = "Shear_n";
291  m_f->m_variables[2] = "Shear_z";
292  m_f->m_variables[3] = "Shear_mag";
293  }
294 }
295 
298  Array<OneD, NekDouble> &mu, NekDouble &lambda)
299 {
300  NekDouble m_mu;
301  int npoints = exp[0]->GetNpoints();
302 
303  if (boost::iequals(m_f->m_variables[0], "u"))
304  {
305  // IncNavierStokesSolver
306  m_mu = m_f->m_session->GetParameter("Kinvis");
307  Vmath::Fill(npoints, m_mu, mu, 1);
308  lambda = 0;
309  }
310  else if (boost::iequals(m_f->m_variables[0], "rho") &&
311  boost::iequals(m_f->m_variables[1], "rhou"))
312  {
313  // CompressibleFlowSolver
314  std::string m_ViscosityType;
315  m_f->m_session->LoadParameter("mu", m_mu, 1.78e-05);
316  m_f->m_session->LoadParameter("lambda", lambda, -2.0 / 3.0);
317  m_f->m_session->LoadSolverInfo("ViscosityType", m_ViscosityType,
318  "Constant");
319 
320  if (m_ViscosityType == "Variable")
321  {
322  // Check equation of state
323  std::string eosType;
324  bool m_idealGas;
325  m_f->m_session->LoadSolverInfo("EquationOfState", eosType,
326  "IdealGas");
327  m_idealGas = boost::iequals(eosType, "IdealGas");
328  ASSERTL0(
329  m_idealGas,
330  "Only IdealGas EOS implemented for Variable ViscosityType");
331 
332  // Get relevant parameters
333  NekDouble m_gamma;
334  NekDouble m_pInf;
336  NekDouble m_gasConstant;
337  NekDouble cv_inv;
338  NekDouble m_Tref;
339  NekDouble m_TRatioSutherland;
340  m_f->m_session->LoadParameter("Gamma", m_gamma, 1.4);
341  m_f->m_session->LoadParameter("pInf", m_pInf, 101325);
342  m_f->m_session->LoadParameter("rhoInf", m_rhoInf, 1.225);
343  m_f->m_session->LoadParameter("GasConstant", m_gasConstant,
344  287.058);
345  m_f->m_session->LoadParameter("Tref", m_Tref, 288.15);
346  m_TRatioSutherland = 110.0 / m_Tref;
347 
348  // Get temperature from flowfield
349  cv_inv = (m_gamma - 1.0) / m_gasConstant;
350  // e = 1/rho ( E - 1/2 ( rhou^2/rho + ... ) )
351  Array<OneD, NekDouble> tmp(npoints, 0.0);
352  Array<OneD, NekDouble> energy(npoints, 0.0);
353  Array<OneD, NekDouble> temperature(npoints, 0.0);
354  Vmath::Vcopy(npoints, exp[m_spacedim + 1]->GetPhys(), 1, energy, 1);
355  for (int i = 0; i < m_spacedim; i++)
356  {
357  // rhou^2
358  Vmath::Vmul(npoints, exp[i + 1]->GetPhys(), 1,
359  exp[i + 1]->GetPhys(), 1, tmp, 1);
360  // rhou^2/rho
361  Vmath::Vdiv(npoints, tmp, 1, exp[0]->GetPhys(), 1, tmp, 1);
362  // 0.5 rhou^2/rho
363  Vmath::Smul(npoints, 0.5, tmp, 1, tmp, 1);
364  // E - 0.5 rhou^2/rho - ...
365  Vmath::Vsub(npoints, energy, 1, tmp, 1, energy, 1);
366  }
367  // rhoe/rho
368  Vmath::Vdiv(npoints, energy, 1, exp[0]->GetPhys(), 1, energy, 1);
369  // T = e/Cv
370  Vmath::Smul(npoints, cv_inv, energy, 1, temperature, 1);
371 
372  // Variable viscosity through the Sutherland's law
373  //
374  // WARNING, if this routine is modified the same must be done in the
375  // CompressibleFlowSolver function in VariableConverter.cpp
376  // (this class should be restructured).
377 
378  const NekDouble C = m_TRatioSutherland;
379  NekDouble mu_star = m_mu;
380  NekDouble T_star = m_pInf / (m_rhoInf * m_gasConstant);
381  NekDouble ratio;
382  for (int i = 0; i < npoints; ++i)
383  {
384  ratio = temperature[i] / T_star;
385  mu[i] = mu_star * ratio * sqrt(ratio) * (1 + C) / (ratio + C);
386  }
387  }
388  else
389  {
390  Vmath::Fill(npoints, m_mu, mu, 1);
391  }
392  }
393  else
394  {
395  // Unknown
396  ASSERTL0(false, "Invalid variables for WSS");
397  }
398 }
399 
403 {
404  int npoints = exp[0]->GetNpoints();
405  if (boost::iequals(m_f->m_variables[0], "u"))
406  {
407  // IncNavierStokesSolver
408  for (int i = 0; i < m_spacedim; ++i)
409  {
410  vel[i] = Array<OneD, NekDouble>(npoints);
411  Vmath::Vcopy(npoints, exp[i]->GetPhys(), 1, vel[i], 1);
412  }
413  }
414  else if (boost::iequals(m_f->m_variables[0], "rho") &&
415  boost::iequals(m_f->m_variables[1], "rhou"))
416  {
417  // CompressibleFlowSolver
418  for (int i = 0; i < m_spacedim; ++i)
419  {
420  vel[i] = Array<OneD, NekDouble>(npoints);
421  Vmath::Vdiv(npoints, exp[i + 1]->GetPhys(), 1, exp[0]->GetPhys(), 1,
422  vel[i], 1);
423  }
424  }
425  else
426  {
427  // Unknown
428  ASSERTL0(false, "Could not identify velocity for WSS");
429  }
430 }
431 
432 } // namespace FieldUtils
433 } // namespace Nektar
NekDouble m_mu
NekDouble m_rhoInf
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
FieldSharedPtr m_f
Field object.
Definition: Module.h:225
This processing module sets up for the boundary field to be extracted.
virtual void Process(po::variables_map &vm)
void GetViscosity(const Array< OneD, MultiRegions::ExpListSharedPtr > exp, Array< OneD, NekDouble > &mu, NekDouble &lambda)
Definition: ProcessWSS.cpp:296
virtual void Process(po::variables_map &vm)
Write mesh to output file.
Definition: ProcessWSS.cpp:62
void GetVelocity(const Array< OneD, MultiRegions::ExpListSharedPtr > exp, Array< OneD, Array< OneD, NekDouble >> &vel)
Definition: ProcessWSS.cpp:400
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:234
std::shared_ptr< Field > FieldSharedPtr
Definition: Field.hpp:989
std::pair< ModuleType, std::string > ModuleKey
Definition: Module.h:285
ModuleFactory & GetModuleFactory()
Definition: Module.cpp:49
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:210
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
double NekDouble
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:534
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:209
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:574
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:359
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:284
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:419
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291