Nektar++
NavierStokesCFEAxisym.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: NavierStokesCFEAxisym.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: Navier Stokes equations in conservative variables
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
37namespace Nektar
38{
39
42 "NavierStokesCFEAxisym", NavierStokesCFEAxisym::create,
43 "Axisymmetric NavierStokes equations in conservative variables.");
44
48 : UnsteadySystem(pSession, pGraph),
49 CompressibleFlowSystem(pSession, pGraph),
50 NavierStokesCFE(pSession, pGraph)
51{
52}
53
55{
56 NavierStokesCFE::v_InitObject(DeclareFields);
57
58 int nVariables = m_fields.size();
59 int npoints = GetNpoints();
61 for (int i = 0; i < nVariables; ++i)
62 {
64 }
65}
66
68 const Array<OneD, Array<OneD, NekDouble>> &inarray,
70 const Array<OneD, Array<OneD, NekDouble>> &pFwd,
71 const Array<OneD, Array<OneD, NekDouble>> &pBwd)
72{
73 int npoints = GetNpoints();
74 int nvariables = inarray.size();
75
76 NavierStokesCFE::v_DoDiffusion(inarray, outarray, pFwd, pBwd);
77
78 for (int i = 0; i < nvariables; ++i)
79 {
80 Vmath::Vadd(npoints, m_viscousForcing[i], 1, outarray[i], 1,
81 outarray[i], 1);
82 }
83}
84
85/**
86 * @brief Return the flux vector for the LDG diffusion problem.
87 * \todo Complete the viscous flux vector
88 */
90 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
91 TensorOfArray3D<NekDouble> &derivativesO1,
92 TensorOfArray3D<NekDouble> &viscousTensor)
93{
94 int i, j;
95 int nVariables = m_fields.size();
96 int nPts = physfield[0].size();
97
98 // 1/r
100 Array<OneD, NekDouble> invR(nPts, 0.0);
101 for (int i = 0; i < 3; i++)
102 {
103 coords[i] = Array<OneD, NekDouble>(nPts);
104 }
105 m_fields[0]->GetCoords(coords[0], coords[1], coords[2]);
106 for (int i = 0; i < nPts; ++i)
107 {
108 if (coords[0][i] < NekConstants::kNekZeroTol)
109 {
110 invR[i] = 0;
111 }
112 else
113 {
114 invR[i] = 1.0 / coords[0][i];
115 }
116 }
117
118 // Stokes hypothesis
119 const NekDouble lambda = -2.0 / 3.0;
120
121 // Auxiliary variables
122 Array<OneD, NekDouble> divVel(nPts, 0.0);
123 Array<OneD, NekDouble> tmp(nPts, 0.0);
124 Array<OneD, NekDouble> mu(nPts, 0.0);
125 Array<OneD, NekDouble> thermalConductivity(nPts, 0.0);
126
127 // Update viscosity and thermal conductivity
128 GetViscosityAndThermalCondFromTemp(physfield[nVariables - 2], mu,
129 thermalConductivity);
130
131 // Velocity divergence = d(u_r)/dr + d(u_z)/dz + u_r/r
132 Vmath::Vadd(nPts, derivativesO1[0][0], 1, derivativesO1[1][1], 1, divVel,
133 1);
134 Vmath::Vvtvp(nPts, physfield[0], 1, invR, 1, divVel, 1, divVel, 1);
135
136 // Velocity divergence scaled by lambda * mu
137 Vmath::Smul(nPts, lambda, divVel, 1, divVel, 1);
138 Vmath::Vmul(nPts, mu, 1, divVel, 1, divVel, 1);
139
140 // Viscous flux vector for the rho equation = 0
141 for (i = 0; i < m_spacedim; ++i)
142 {
143 Vmath::Zero(nPts, viscousTensor[i][0], 1);
144 }
145
146 // Viscous stress tensor (for the momentum equations)
147
148 for (i = 0; i < 2; ++i)
149 {
150 for (j = i; j < 2; ++j)
151 {
152 Vmath::Vadd(nPts, derivativesO1[i][j], 1, derivativesO1[j][i], 1,
153 viscousTensor[i][j + 1], 1);
154
155 Vmath::Vmul(nPts, mu, 1, viscousTensor[i][j + 1], 1,
156 viscousTensor[i][j + 1], 1);
157
158 if (i == j)
159 {
160 // Add divergence term to diagonal
161 Vmath::Vadd(nPts, viscousTensor[i][j + 1], 1, divVel, 1,
162 viscousTensor[i][j + 1], 1);
163 }
164 else
165 {
166 // Copy to make symmetric
167 Vmath::Vcopy(nPts, viscousTensor[i][j + 1], 1,
168 viscousTensor[j][i + 1], 1);
169 }
170 }
171 }
172 // Swirl case
173 if (m_spacedim == 3)
174 {
175 // Tau_theta_theta = mu ( 2*u_r/r - 2/3*div(u))
176 Vmath::Vmul(nPts, physfield[0], 1, invR, 1, viscousTensor[2][3], 1);
177 Vmath::Smul(nPts, 2.0, viscousTensor[2][3], 1, viscousTensor[2][3], 1);
178 Vmath::Vmul(nPts, mu, 1, viscousTensor[2][3], 1, viscousTensor[2][3],
179 1);
180 Vmath::Vadd(nPts, viscousTensor[2][3], 1, divVel, 1,
181 viscousTensor[2][3], 1);
182
183 // Tau_r_theta = mu (-u_theta/r + d(u_theta)/dr )
184 Vmath::Vmul(nPts, physfield[2], 1, invR, 1, viscousTensor[2][1], 1);
185 Vmath::Smul(nPts, -1.0, viscousTensor[2][1], 1, viscousTensor[2][1], 1);
186 Vmath::Vadd(nPts, derivativesO1[0][2], 1, viscousTensor[2][1], 1,
187 viscousTensor[2][1], 1);
188 Vmath::Vmul(nPts, mu, 1, viscousTensor[2][1], 1, viscousTensor[2][1],
189 1);
190 Vmath::Vcopy(nPts, viscousTensor[2][1], 1, viscousTensor[0][3], 1);
191
192 // Tau_z_theta = mu (d(u_theta)/dz )
193 Vmath::Vmul(nPts, mu, 1, derivativesO1[1][2], 1, viscousTensor[2][2],
194 1);
195 Vmath::Vcopy(nPts, viscousTensor[2][2], 1, viscousTensor[1][3], 1);
196 }
197
198 // Terms for the energy equation
199 for (i = 0; i < m_spacedim; ++i)
200 {
201 Vmath::Zero(nPts, viscousTensor[i][m_spacedim + 1], 1);
202 // u_j * tau_ij
203 for (j = 0; j < m_spacedim; ++j)
204 {
205 Vmath::Vvtvp(nPts, physfield[j], 1, viscousTensor[i][j + 1], 1,
206 viscousTensor[i][m_spacedim + 1], 1,
207 viscousTensor[i][m_spacedim + 1], 1);
208 }
209 // Add k*T_i
210 if (i != 2)
211 {
212 Vmath::Vvtvp(nPts, thermalConductivity, 1,
213 derivativesO1[i][m_spacedim], 1,
214 viscousTensor[i][m_spacedim + 1], 1,
215 viscousTensor[i][m_spacedim + 1], 1);
216 }
217 else
218 {
219 Vmath::Vmul(nPts, derivativesO1[i][m_spacedim], 1, invR, 1, tmp, 1);
220 Vmath::Vvtvp(nPts, thermalConductivity, 1, tmp, 1,
221 viscousTensor[i][m_spacedim + 1], 1,
222 viscousTensor[i][m_spacedim + 1], 1);
223 }
224 }
225
226 // Update viscous forcing
227 // r-momentum: F = 1/r * (tau_rr - tau_theta_theta)
228 if (m_spacedim == 3)
229 {
230 Vmath::Vsub(nPts, viscousTensor[0][1], 1, viscousTensor[2][3], 1,
231 m_viscousForcing[1], 1);
232 Vmath::Vmul(nPts, m_viscousForcing[1], 1, invR, 1, m_viscousForcing[1],
233 1);
234 }
235 else
236 {
237 Vmath::Vmul(nPts, viscousTensor[0][1], 1, invR, 1, m_viscousForcing[1],
238 1);
239 }
240
241 // z-momentum: F = 1/r * tau_r_z
242 Vmath::Vmul(nPts, viscousTensor[0][2], 1, invR, 1, m_viscousForcing[2], 1);
243
244 // Theta_momentum: F = 2* tau_r_theta
245 if (m_spacedim == 3)
246 {
247 Vmath::Vmul(nPts, viscousTensor[0][3], 1, invR, 1, m_viscousForcing[3],
248 1);
249 Vmath::Smul(nPts, 2.0, m_viscousForcing[3], 1, m_viscousForcing[3], 1);
250 }
251
252 // Energy: F = 1/r* viscousTensor_T_r
253 Vmath::Vmul(nPts, viscousTensor[0][m_spacedim + 1], 1, invR, 1,
255}
256
257} // namespace Nektar
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void v_InitObject(bool DeclareFields=true) override
Initialization object for CompressibleFlowSystem class.
void v_DoDiffusion(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd) override
Array< OneD, Array< OneD, NekDouble > > m_viscousForcing
void v_GetViscousFluxVector(const Array< OneD, const Array< OneD, NekDouble > > &physfield, TensorOfArray3D< NekDouble > &derivatives, TensorOfArray3D< NekDouble > &viscousTensor) override
Return the flux vector for the LDG diffusion problem.
NavierStokesCFEAxisym(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
void v_InitObject(bool DeclareField=true) override
Initialization object for CompressibleFlowSystem class.
void v_DoDiffusion(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const Array< OneD, Array< OneD, NekDouble > > &pFwd, const Array< OneD, Array< OneD, NekDouble > > &pBwd) override
void GetViscosityAndThermalCondFromTemp(const Array< OneD, NekDouble > &temperature, Array< OneD, NekDouble > &mu, Array< OneD, NekDouble > &thermalCond)
Update viscosity todo: add artificial viscosity here.
int m_spacedim
Spatial dimension (>= expansion dim).
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNpoints()
Base class for unsteady solvers.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
static const NekDouble kNekZeroTol
EquationSystemFactory & GetEquationSystemFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
double NekDouble
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.hpp:72
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.hpp:366
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.hpp:180
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.hpp:100
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
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.hpp:220