1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
|
/* This file is part of the OpenLB library
*
* Copyright (C) 2016-2017 Davide Dapelo, Mathias J. Krause
* OpenLB e-mail contact: info@openlb.net
* The most recent release of OpenLB can be downloaded at
* <http://www.openlb.net/>
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public
* License along with this program; if not, write to the Free
* Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301, USA.
*/
/** \file
* Helper functions for the implementation of the
* Guo-ZhapLB dynamics.
*/
#ifndef LB_GUOZHAO_HELPERS_H
#define LB_GUOZHAO_HELPERS_H
#include "dynamics/guoZhaoLatticeDescriptors.h"
#include "core/cell.h"
#include "core/util.h"
namespace olb {
// Forward declarations
template<typename T, typename DESCRIPTORBASE> struct GuoZhaoLbDynamicsHelpers;
template<typename T, typename DESCRIPTOR> struct GuoZhaoLbExternalHelpers;
/// This structure forwards the calls to the appropriate Guo Zhao helper class
template<typename T, typename DESCRIPTOR>
struct GuoZhaoLbHelpers {
static T equilibrium(int iPop, T epsilon, T rho, const T u[DESCRIPTOR::d], const T uSqr) {
return GuoZhaoLbDynamicsHelpers<T,typename DESCRIPTOR::BaseDescriptor>
::equilibrium(iPop, epsilon, rho, u, uSqr);
}
static T bgkCollision(Cell<T,DESCRIPTOR>& cell, T const& epsilon, T const& rho, const T u[DESCRIPTOR::d], T const& omega) {
return GuoZhaoLbDynamicsHelpers<T,typename DESCRIPTOR::BaseDescriptor>
::bgkCollision(cell, epsilon, rho, u, omega);
}
static void updateGuoZhaoForce(Cell<T,DESCRIPTOR>& cell, const T u[DESCRIPTOR::d]) {
GuoZhaoLbExternalHelpers<T,DESCRIPTOR>::updateGuoZhaoForce(cell, u);
}
static void addExternalForce(Cell<T,DESCRIPTOR>& cell, const T u[DESCRIPTOR::d], T omega, T rho, T epsilon=(T)1)
{
GuoZhaoLbExternalHelpers<T,DESCRIPTOR>::addExternalForce(cell, u, omega, rho, epsilon);
}
}; // struct GuoZhaoLbHelpers
/// All Guo Zhao helper functions are inside this structure
template<typename T, typename DESCRIPTORBASE>
struct GuoZhaoLbDynamicsHelpers {
/// Computation of Guo Zhao equilibrium distribution - original (compressible) formulation following Guo and Zhao (2002).
static T forceEquilibrium(int iPop, T epsilon, T rho, const T u[DESCRIPTORBASE::d], const T force[DESCRIPTORBASE::d], T nu) {
}
/// Computation of Guo Zhao equilibrium distribution - original (compressible) formulation following Guo and Zhao (2002).
static T equilibrium(int iPop, T epsilon, T rho, const T u[DESCRIPTORBASE::d], const T uSqr) {
T c_u = T();
for (int iD=0; iD < DESCRIPTORBASE::d; ++iD) {
c_u += descriptors::c<DESCRIPTORBASE>(iPop,iD)*u[iD];
}
return rho * descriptors::t<T,DESCRIPTORBASE>(iPop) * (
(T)1 + descriptors::invCs2<T,DESCRIPTORBASE>() * c_u +
descriptors::invCs2<T,DESCRIPTORBASE>() * descriptors::invCs2<T,DESCRIPTORBASE>()/((T)2*epsilon) * c_u*c_u -
descriptors::invCs2<T,DESCRIPTORBASE>()/((T)2*epsilon) * uSqr
) - descriptors::t<T,DESCRIPTORBASE>(iPop);
}
/// Guo Zhao BGK collision step
static T bgkCollision(CellBase<T,DESCRIPTORBASE>& cell, T const& epsilon, T const& rho, const T u[DESCRIPTORBASE::d], T const& omega) {
const T uSqr = util::normSqr<T,DESCRIPTORBASE::d>(u);
for (int iPop=0; iPop < DESCRIPTORBASE::q; ++iPop) {
cell[iPop] *= (T)1-omega;
cell[iPop] += omega * GuoZhaoLbDynamicsHelpers<T,DESCRIPTORBASE>::equilibrium (
iPop, epsilon, rho, u, uSqr );
}
return uSqr;
}
}; // struct GuoZhaoLbDynamicsHelpers
/// Helper functions for dynamics that access external field
template<typename T, typename DESCRIPTOR>
/// Updates Guo Zhao porous force
struct GuoZhaoLbExternalHelpers {
static void updateGuoZhaoForce(Cell<T,DESCRIPTOR>& cell, const T u[DESCRIPTOR::d]) {
T epsilon = *cell.template getFieldPointer<descriptors::EPSILON>();
T k = *cell.template getFieldPointer<descriptors::K>();
T nu = *cell.template getFieldPointer<descriptors::NU>();
T* force[DESCRIPTOR::d];
for (int iDim=0; iDim <DESCRIPTOR::d; iDim++) {
force[iDim] = cell.template getFieldPointer<descriptors::FORCE>() + iDim;
}
T bodyF[DESCRIPTOR::d];
for (int iDim=0; iDim <DESCRIPTOR::d; iDim++) {
bodyF[iDim] = cell.template getFieldPointer<descriptors::BODY_FORCE>()[iDim];
}
T uMag = 0.;
for (int iDim=0; iDim <DESCRIPTOR::d; iDim++) {
uMag += u[iDim]*u[iDim];
}
uMag = sqrt(uMag);
T Fe = 0.; //1.75/sqrt(150.*pow(epsilon,3));
// Linear Darcy term, nonlinear Forchheimer term and body force
for (int iDim=0; iDim <DESCRIPTOR::d; iDim++) {
*force[iDim] = -u[iDim]*epsilon*nu/k - epsilon*Fe/sqrt(k)*uMag*u[iDim] + bodyF[iDim]*epsilon;
}
}
/// Add a force term scaled by physical porosity epsilon after BGK collision
static void addExternalForce(Cell<T,DESCRIPTOR>& cell, const T u[DESCRIPTOR::d], T omega, T rho, T epsilon)
{
T* force = cell.template getFieldPointer<descriptors::FORCE>();
for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
T c_u = T();
for (int iD=0; iD < DESCRIPTOR::d; ++iD) {
c_u += descriptors::c<DESCRIPTOR>(iPop,iD)*u[iD];
}
c_u *= descriptors::invCs2<T,DESCRIPTOR>()*descriptors::invCs2<T,DESCRIPTOR>()/epsilon;
T forceTerm = T();
for (int iD=0; iD < DESCRIPTOR::d; ++iD) {
forceTerm +=
( (epsilon*(T)descriptors::c<DESCRIPTOR>(iPop,iD)-u[iD]) * descriptors::invCs2<T,DESCRIPTOR>()/epsilon
+ c_u * descriptors::c<DESCRIPTOR>(iPop,iD)
)
* force[iD];
}
forceTerm *= descriptors::t<T,DESCRIPTOR>(iPop);
forceTerm *= T(1) - omega/T(2);
forceTerm *= rho;
cell[iPop] += forceTerm;
}
}
}; // struct GuoZhaoLbExternalHelpers
} // namespace olb
#endif
|
