atmBoundaryLayer.C 7.52 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
OpenFOAM bot's avatar
OpenFOAM bot committed
6
7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2014-2016 OpenFOAM Foundation
9
    Copyright (C) 2018-2020 OpenCFD Ltd.
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
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM 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 3 of the License, or
    (at your option) any later version.

    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "atmBoundaryLayer.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

38
atmBoundaryLayer::atmBoundaryLayer(const Time& time, const polyPatch& pp)
39
:
40
41
42
43
44
45
    initABL_(false),
    kappa_(0.41),
    Cmu_(0.09),
    C1_(0.0),
    C2_(1.0),
    ppMin_((boundBox(pp.points())).min()),
46
47
48
49
50
51
    time_(time),
    patch_(pp),
    flowDir_(time, "flowDir"),
    zDir_(time, "zDir"),
    Uref_(time, "Uref"),
    Zref_(time, "Zref"),
52
53
    z0_(nullptr),
    d_(nullptr)
54
55
56
{}


57
58
59
60
61
62
atmBoundaryLayer::atmBoundaryLayer
(
    const Time& time,
    const polyPatch& pp,
    const dictionary& dict
)
63
:
64
    initABL_(dict.getOrDefault<bool>("initABL", true)),
65
66
67
68
69
70
71
72
    kappa_
    (
        dict.getCheckOrDefault<scalar>("kappa", 0.41, scalarMinMax::ge(SMALL))
    ),
    Cmu_(dict.getCheckOrDefault<scalar>("Cmu", 0.09, scalarMinMax::ge(SMALL))),
    C1_(dict.getOrDefault("C1", 0.0)),
    C2_(dict.getOrDefault("C2", 1.0)),
    ppMin_((boundBox(pp.points())).min()),
73
74
75
76
77
78
79
    time_(time),
    patch_(pp),
    flowDir_(TimeFunction1<vector>(time, "flowDir", dict)),
    zDir_(TimeFunction1<vector>(time, "zDir", dict)),
    Uref_(TimeFunction1<scalar>(time, "Uref", dict)),
    Zref_(TimeFunction1<scalar>(time, "Zref", dict)),
    z0_(PatchFunction1<scalar>::New(pp, "z0", dict)),
80
    d_(PatchFunction1<scalar>::New(pp, "d", dict))
81
{}
82
83
84
85


atmBoundaryLayer::atmBoundaryLayer
(
86
87
    const atmBoundaryLayer& abl,
    const fvPatch& patch,
88
89
90
    const fvPatchFieldMapper& mapper
)
:
91
92
93
94
95
96
    initABL_(abl.initABL_),
    kappa_(abl.kappa_),
    Cmu_(abl.Cmu_),
    C1_(abl.C1_),
    C2_(abl.C2_),
    ppMin_(abl.ppMin_),
97
98
99
100
101
102
103
    time_(abl.time_),
    patch_(patch.patch()),
    flowDir_(abl.flowDir_),
    zDir_(abl.zDir_),
    Uref_(abl.Uref_),
    Zref_(abl.Zref_),
    z0_(abl.z0_.clone(patch_)),
104
    d_(abl.d_.clone(patch_))
105
106
107
{}


108
atmBoundaryLayer::atmBoundaryLayer(const atmBoundaryLayer& abl)
109
:
110
111
112
113
114
115
    initABL_(abl.initABL_),
    kappa_(abl.kappa_),
    Cmu_(abl.Cmu_),
    C1_(abl.C1_),
    C2_(abl.C2_),
    ppMin_(abl.ppMin_),
116
117
118
119
120
121
122
    time_(abl.time_),
    patch_(abl.patch_),
    flowDir_(abl.flowDir_),
    zDir_(abl.zDir_),
    Uref_(abl.Uref_),
    Zref_(abl.Zref_),
    z0_(abl.z0_.clone(patch_)),
123
    d_(abl.d_.clone(patch_))
124
125
126
127
128
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

129
130
131
vector atmBoundaryLayer::flowDir() const
{
    const scalar t = time_.timeOutputValue();
132
    const vector dir(flowDir_.value(t));
133
134
135
136
137
    const scalar magDir = mag(dir);

    if (magDir < SMALL)
    {
        FatalErrorInFunction
138
            << "magnitude of " << flowDir_.name() << " = " << magDir
139
140
141
142
143
144
145
146
147
148
149
            << " vector must be greater than zero"
            << abort(FatalError);
    }

    return dir/magDir;
}


vector atmBoundaryLayer::zDir() const
{
    const scalar t = time_.timeOutputValue();
150
    const vector dir(zDir_.value(t));
151
152
153
154
155
    const scalar magDir = mag(dir);

    if (magDir < SMALL)
    {
        FatalErrorInFunction
156
            << "magnitude of " << zDir_.name() << " = " << magDir
157
158
159
160
161
162
163
164
165
166
167
168
169
170
            << " vector must be greater than zero"
            << abort(FatalError);
    }

    return dir/magDir;
}


tmp<scalarField> atmBoundaryLayer::Ustar(const scalarField& z0) const
{
    const scalar t = time_.timeOutputValue();
    const scalar Uref = Uref_.value(t);
    const scalar Zref = Zref_.value(t);

171
172
173
174
175
176
177
178
179
    if (Zref < 0)
    {
        FatalErrorInFunction
            << "Negative entry in " << Zref_.name() << " = " << Zref
            << abort(FatalError);
    }

    // (derived from RH:Eq. 6, HW:Eq. 5)
    return kappa_*Uref/log((Zref + z0)/z0);
180
181
182
183
}


void atmBoundaryLayer::autoMap(const fvPatchFieldMapper& mapper)
184
{
185
186
187
188
189
190
191
192
    if (z0_)
    {
        z0_->autoMap(mapper);
    }
    if (d_)
    {
        d_->autoMap(mapper);
    }
193
194
195
196
197
}


void atmBoundaryLayer::rmap
(
198
    const atmBoundaryLayer& abl,
199
200
201
    const labelList& addr
)
{
202
203
204
205
206
207
208
209
    if (z0_)
    {
        z0_->rmap(abl.z0_(), addr);
    }
    if (d_)
    {
        d_->rmap(abl.d_(), addr);
    }
210
211
212
}


213
tmp<vectorField> atmBoundaryLayer::U(const vectorField& pCf) const
214
{
215
    const scalar t = time_.timeOutputValue();
216
    const scalarField d(d_->value(t));
217
    const scalarField z0(max(z0_->value(t), ROOTVSMALL));
218
    const scalar groundMin = zDir() & ppMin_;
219

220
221
222
223
224
    // (YGCJ:Table 1, RH:Eq. 6, HW:Eq. 5)
    scalarField Un
    (
        (Ustar(z0)/kappa_)*log(((zDir() & pCf) - groundMin - d + z0)/z0)
    );
225

226
    return flowDir()*Un;
227
228
229
}


230
tmp<scalarField> atmBoundaryLayer::k(const vectorField& pCf) const
231
{
232
    const scalar t = time_.timeOutputValue();
233
    const scalarField d(d_->value(t));
234
    const scalarField z0(max(z0_->value(t), ROOTVSMALL));
235
    const scalar groundMin = zDir() & ppMin_;
236

237
238
239
240
    // (YGCJ:Eq. 21; RH:Eq. 7, HW:Eq. 6 when C1=0 and C2=1)
    return
        sqr(Ustar(z0))/sqrt(Cmu_)
       *sqrt(C1_*log(((zDir() & pCf) - groundMin - d + z0)/z0) + C2_);
241
242
243
}


244
tmp<scalarField> atmBoundaryLayer::epsilon(const vectorField& pCf) const
245
{
246
    const scalar t = time_.timeOutputValue();
247
    const scalarField d(d_->value(t));
248
    const scalarField z0(max(z0_->value(t), ROOTVSMALL));
249
    const scalar groundMin = zDir() & ppMin_;
250

251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
    // (YGCJ:Eq. 22; RH:Eq. 8, HW:Eq. 7 when C1=0 and C2=1)
    return
        pow3(Ustar(z0))/(kappa_*((zDir() & pCf) - groundMin - d + z0))
       *sqrt(C1_*log(((zDir() & pCf) - groundMin - d + z0)/z0) + C2_);
}


tmp<scalarField> atmBoundaryLayer::omega(const vectorField& pCf) const
{
    const scalar t = time_.timeOutputValue();
    const scalarField d(d_->value(t));
    const scalarField z0(max(z0_->value(t), ROOTVSMALL));
    const scalar groundMin = zDir() & ppMin_;

    // (YGJ:Eq. 13)
    return Ustar(z0)/(kappa_*sqrt(Cmu_)*((zDir() & pCf) - groundMin - d + z0));
267
268
269
270
271
}


void atmBoundaryLayer::write(Ostream& os) const
{
272
    os.writeEntry("initABL", initABL_);
273
274
    os.writeEntry("kappa", kappa_);
    os.writeEntry("Cmu", Cmu_);
275
276
277
278
    os.writeEntry("C1", C1_);
    os.writeEntry("C2", C2_);
    flowDir_.writeData(os);
    zDir_.writeData(os);
279
280
    Uref_.writeData(os);
    Zref_.writeData(os);
281
282
283
284
285
286
287
288
    if (z0_)
    {
        z0_->writeData(os) ;
    }
    if (d_)
    {
        d_->writeData(os);
    }
289
290
291
292
293
294
295
296
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //