Skip to content

Commit 41be0c8

Browse files
author
Bart Doekemeijer
authored
Merge pull request #9 from NREL/master
Merge changes from master branch
2 parents 220a6f1 + 4d290da commit 41be0c8

File tree

12 files changed

+343
-26
lines changed

12 files changed

+343
-26
lines changed

‎applications/solvers/incompressible/windEnergy/ABLSolver/UEqn.H

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
// Solve the momentum equation
22

33
#include "computeCoriolisForce.H"
4+
#include "computeSpongeForce.H"
45

56
#include "computeBuoyancyTerm.H"
67

@@ -11,6 +12,7 @@
1112
+ turbulence->divDevReff(U) // momentum flux
1213
+ fvc::div(Rwall)
1314
- fCoriolis // Coriolis force
15+
- fSponge // Sponge layer damping
1416
- SourceU // mesoscale source terms
1517
);
1618

‎applications/solvers/incompressible/windEnergy/ABLSolver/computeSpongeForce.H

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
// Compute the sponge layer damping force
2+
if (spongeLayerType == "Rayleigh")
3+
{
4+
fSponge = spongeLayerViscosity * (spongeLayerReferenceVelocity - U);
5+
}
6+
else if (spongeLayerType == "viscous")
7+
{
8+
fSponge = fvc::laplacian(spongeLayerViscosity,U);
9+
}

‎applications/solvers/incompressible/windEnergy/ABLSolver/createFields.H

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -78,6 +78,83 @@
7878
dimensionedVector("fCoriolis",dimensionSet(0, 1, -2, 0, 0, 0, 0),vector::zero)
7979
);
8080

81+
// Create sponge layer viscosity field
82+
Info << "Creating sponge layer viscosity field, spongeLayerViscosity..." << endl;
83+
label lengthDimension;
84+
if (spongeLayerType == "Rayleigh")
85+
{
86+
lengthDimension = 0;
87+
}
88+
else
89+
{
90+
lengthDimension = 2;
91+
}
92+
93+
volScalarField spongeLayerViscosity
94+
(
95+
IOobject
96+
(
97+
"spongeLayerViscosity",
98+
runTime.timeName(),
99+
mesh,
100+
IOobject::NO_READ,
101+
IOobject::AUTO_WRITE
102+
),
103+
mesh,
104+
dimensionedScalar("spongeLayerViscosity",dimensionSet(0, lengthDimension, -1, 0, 0, 0, 0),spongeLayerViscosityTop)
105+
);
106+
107+
forAll(mesh.cells(),cellI)
108+
{
109+
scalar z = mesh.C()[cellI][upIndex];
110+
spongeLayerViscosity[cellI] *= 0.5 *
111+
(
112+
1.0 - Foam::cos
113+
(
114+
Foam::constant::mathematical::pi *
115+
max( (z - spongeLayerBaseHeight)/(spongeLayerTopHeight - spongeLayerBaseHeight) , 0.0 )
116+
)
117+
118+
);
119+
}
120+
121+
forAll(spongeLayerViscosity.boundaryField(),i)
122+
{
123+
if ( !mesh.boundary()[i].coupled() )
124+
{
125+
forAll(spongeLayerViscosity.boundaryField()[i],j)
126+
{
127+
scalar z = mesh.boundary()[i].Cf()[j].z();
128+
spongeLayerViscosity.boundaryField()[i][j] *= 0.5 *
129+
(
130+
1.0 - Foam::cos
131+
(
132+
Foam::constant::mathematical::pi *
133+
max( (z - spongeLayerBaseHeight)/(spongeLayerTopHeight - spongeLayerBaseHeight) , 0.0 )
134+
)
135+
136+
);
137+
}
138+
}
139+
}
140+
141+
// Create sponge layer force vector
142+
Info << "Creating sponge layer force vector, fSponge..." << endl;
143+
volVectorField fSponge
144+
(
145+
IOobject
146+
(
147+
"fSponge",
148+
runTime.timeName(),
149+
mesh,
150+
IOobject::READ_IF_PRESENT,
151+
IOobject::NO_WRITE
152+
),
153+
mesh,
154+
dimensionedVector("fSponge",dimensionSet(0, 1, -2, 0, 0, 0, 0),vector::zero)
155+
);
156+
157+
81158
// Create and calculate the Boussinesq density field for computing buoyancy forces
82159
Info << "Creating kinematic (Boussinesq) density field, rhok..." << endl;
83160
volScalarField rhok

‎applications/solvers/incompressible/windEnergy/ABLSolver/pEqn.H

Lines changed: 12 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,11 +23,22 @@
2323
+ rAUf*fvc::ddtCorr(U, phi)
2424
);
2525

26-
// This part balances global mass flux.
26+
// This part balances global mass flux. It does it in a temporary field, and then
27+
// applies the correction indirectly by setting the pressure gradient to be used in
28+
// the fixedFluxPressure boundary condition on p_rgh, or directly if the zeroGradient
29+
// boundary condition on p_rgh is used.
2730
surfaceScalarField phiFixedFlux = phiHbyA;
2831

2932
adjustPhi(phiFixedFlux, U, p_rgh);
3033

34+
forAll(p_rgh.boundaryField(), patchi)
35+
{
36+
if (isA<zeroGradientFvPatchScalarField>(p_rgh.boundaryField()[patchi]))
37+
{
38+
phiHbyA.boundaryField()[patchi] = phiFixedFlux.boundaryField()[patchi];
39+
}
40+
}
41+
3142
phiHbyA += phig;
3243

3344
// Update the fixedFluxPressure BCs to ensure flux consistency

‎applications/solvers/incompressible/windEnergy/ABLSolver/readABLProperties.H

Lines changed: 45 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -197,8 +197,6 @@
197197

198198

199199

200-
201-
202200
// PROPERTIES CONCERNING THE WAY IN WHICH PERTURBATION PRESSURE IS DEFINED
203201

204202
// Options for defining the perturbation pressure:
@@ -243,8 +241,51 @@
243241

244242

245243

246-
247-
244+
// PROPERTIES CONCERNING SPONGE LAYER
245+
246+
// Specify the type of sponge layer to use. The
247+
// possible types are "Rayleigh", "viscous" or "none".
248+
// - The "Rayleigh" type means that the damping term is computed as nu*(u_ref-u)
249+
// The viscosity coefficient nu has dimensions of 1/s
250+
// - The "viscous" type means that the damping term is computed as nu * Lapl(u)
251+
// The viscosity coefficient nu has dimensions of m**2/s
252+
// - The "none" type means no damping is added
253+
word spongeLayerType(ABLProperties.lookupOrDefault<word>("spongeLayerType","none"));
254+
255+
// Sponge layer base height
256+
scalar spongeLayerBaseHeight(ABLProperties.lookupOrDefault<scalar>("spongeLayerBaseHeight",0.0));
257+
258+
// Sponge layer top height
259+
scalar spongeLayerTopHeight(ABLProperties.lookupOrDefault<scalar>("spongeLayerTopHeight",10000.0));
260+
261+
// Sponge layer viscosity at the top boundary
262+
scalar spongeLayerViscosityTop(ABLProperties.lookupOrDefault<scalar>("spongeLayerViscosityTop",0.0));
263+
264+
265+
// Create sponge layer reference velocity
266+
scalar spongeLayerUx(ABLProperties.lookupOrDefault<scalar>("spongeLayerUx",0.0));
267+
scalar spongeLayerUy(ABLProperties.lookupOrDefault<scalar>("spongeLayerUy",0.0));
268+
vector Uref_;
269+
Uref_.x() = spongeLayerUx;
270+
Uref_.y() = spongeLayerUy;
271+
Uref_.z() = 0.0;
272+
uniformDimensionedVectorField spongeLayerReferenceVelocity
273+
(
274+
IOobject
275+
(
276+
"spongeLayerReferenceVelocity",
277+
runTime.constant(),
278+
mesh,
279+
IOobject::NO_READ,
280+
IOobject::NO_WRITE
281+
),
282+
dimensionedVector("spongeLayerReferenceVelocity",dimensionSet(0, 1, -1, 0, 0, 0, 0),Uref_)
283+
);
284+
285+
if (spongeLayerType == "Rayleigh")
286+
{
287+
Info << spongeLayerReferenceVelocity << endl;
288+
}
248289

249290

250291
// PROPERTIES CONCERNING GATHERING STATISTICS

‎applications/solvers/incompressible/windEnergy/ABLTerrainSolver/UEqn.H

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,8 @@
11
// Solve the momentum equation
2+
23

34
#include "computeCoriolisForce.H"
5+
#include "computeSpongeForce.H"
46

57
#include "computeBuoyancyTerm.H"
68

@@ -11,6 +13,7 @@
1113
+ turbulence->divDevReff(U) // momentum flux
1214
+ fvc::div(Rwall)
1315
- fCoriolis // Coriolis force
16+
- fSponge // Sponge layer damping
1417
+ gradP // driving pressure gradient
1518
);
1619

‎applications/solvers/incompressible/windEnergy/ABLTerrainSolver/computeSpongeForce.H

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
// Compute the sponge layer damping force
2+
if (spongeLayerType == "Rayleigh")
3+
{
4+
fSponge = spongeLayerViscosity * (spongeLayerReferenceVelocity - U);
5+
}
6+
else if (spongeLayerType == "viscous")
7+
{
8+
fSponge = fvc::laplacian(spongeLayerViscosity,U);
9+
}

‎applications/solvers/incompressible/windEnergy/ABLTerrainSolver/createFields.H

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -139,6 +139,83 @@
139139
dimensionedVector("fCoriolis",dimensionSet(0, 1, -2, 0, 0, 0, 0),vector::zero)
140140
);
141141

142+
143+
// Create sponge layer viscosity field
144+
Info << "Creating sponge layer viscosity field, spongeLayerViscosity..." << endl;
145+
label lengthDimension;
146+
if (spongeLayerType == "Rayleigh")
147+
{
148+
lengthDimension = 0;
149+
}
150+
else
151+
{
152+
lengthDimension = 2;
153+
}
154+
155+
volScalarField spongeLayerViscosity
156+
(
157+
IOobject
158+
(
159+
"spongeLayerViscosity",
160+
runTime.timeName(),
161+
mesh,
162+
IOobject::NO_READ,
163+
IOobject::AUTO_WRITE
164+
),
165+
mesh,
166+
dimensionedScalar("spongeLayerViscosity",dimensionSet(0, lengthDimension, -1, 0, 0, 0, 0),spongeLayerViscosityTop)
167+
);
168+
169+
forAll(mesh.cells(),cellI)
170+
{
171+
scalar z = mesh.C()[cellI][upIndex];
172+
spongeLayerViscosity[cellI] *= 0.5 *
173+
(
174+
1.0 - Foam::cos
175+
(
176+
Foam::constant::mathematical::pi *
177+
max( (z - spongeLayerBaseHeight)/(spongeLayerTopHeight - spongeLayerBaseHeight) , 0.0 )
178+
)
179+
180+
);
181+
}
182+
183+
forAll(spongeLayerViscosity.boundaryField(),i)
184+
{
185+
if ( !mesh.boundary()[i].coupled() )
186+
{
187+
forAll(spongeLayerViscosity.boundaryField()[i],j)
188+
{
189+
scalar z = mesh.boundary()[i].Cf()[j].z();
190+
spongeLayerViscosity.boundaryField()[i][j] *= 0.5 *
191+
(
192+
1.0 - Foam::cos
193+
(
194+
Foam::constant::mathematical::pi *
195+
max( (z - spongeLayerBaseHeight)/(spongeLayerTopHeight - spongeLayerBaseHeight) , 0.0 )
196+
)
197+
198+
);
199+
}
200+
}
201+
}
202+
203+
// Create sponge layer force vector
204+
Info << "Creating sponge layer force vector, fSponge..." << endl;
205+
volVectorField fSponge
206+
(
207+
IOobject
208+
(
209+
"fSponge",
210+
runTime.timeName(),
211+
mesh,
212+
IOobject::READ_IF_PRESENT,
213+
IOobject::NO_WRITE
214+
),
215+
mesh,
216+
dimensionedVector("fSponge",dimensionSet(0, 1, -2, 0, 0, 0, 0),vector::zero)
217+
);
218+
142219
// Create and calculate the Boussinesq density field for computing buoyancy forces
143220
Info << "Creating kinematic (Boussinesq) density field, rhok..." << endl;
144221
volScalarField rhok

‎applications/solvers/incompressible/windEnergy/ABLTerrainSolver/pEqn.H

Lines changed: 12 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,11 +23,22 @@
2323
+ rAUf*fvc::ddtCorr(U, phi)
2424
);
2525

26-
// This part balances global mass flux.
26+
// This part balances global mass flux. It does it in a temporary field, and then
27+
// applies the correction indirectly by setting the pressure gradient to be used in
28+
// the fixedFluxPressure boundary condition on p_rgh, or directly if the zeroGradient
29+
// boundary condition on p_rgh is used.
2730
surfaceScalarField phiFixedFlux = phiHbyA;
2831

2932
adjustPhi(phiFixedFlux, U, p_rgh);
3033

34+
forAll(p_rgh.boundaryField(), patchi)
35+
{
36+
if (isA<zeroGradientFvPatchScalarField>(p_rgh.boundaryField()[patchi]))
37+
{
38+
phiHbyA.boundaryField()[patchi] = phiFixedFlux.boundaryField()[patchi];
39+
}
40+
}
41+
3142
phiHbyA += phig;
3243

3344
// Update the fixedFluxPressure BCs to ensure flux consistency

‎applications/solvers/incompressible/windEnergy/ABLTerrainSolver/readABLProperties.H

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -102,6 +102,51 @@
102102
Info << Omega << endl;
103103

104104

105+
// PROPERTIES CONCERNING SPONGE LAYER
106+
107+
// Specify the type of sponge layer to use. The
108+
// possible types are "Rayleigh", "viscous" or "none".
109+
// - The "Rayleigh" type means that the damping term is computed as nu*(u_ref-u)
110+
// The viscosity coefficient nu has dimensions of 1/s
111+
// - The "viscous" type means that the damping term is computed as nu * Lapl(u)
112+
// The viscosity coefficient nu has dimensions of m**2/s
113+
// - The "none" type means no damping is added
114+
word spongeLayerType(ABLProperties.lookupOrDefault<word>("spongeLayerType","none"));
115+
116+
// Sponge layer base height
117+
scalar spongeLayerBaseHeight(ABLProperties.lookupOrDefault<scalar>("spongeLayerBaseHeight",0.0));
118+
119+
// Sponge layer top height
120+
scalar spongeLayerTopHeight(ABLProperties.lookupOrDefault<scalar>("spongeLayerTopHeight",10000.0));
121+
122+
// Sponge layer viscosity at the top boundary
123+
scalar spongeLayerViscosityTop(ABLProperties.lookupOrDefault<scalar>("spongeLayerViscosityTop",0.0));
124+
125+
126+
// Create sponge layer reference velocity
127+
scalar spongeLayerUx(ABLProperties.lookupOrDefault<scalar>("spongeLayerUx",0.0));
128+
scalar spongeLayerUy(ABLProperties.lookupOrDefault<scalar>("spongeLayerUy",0.0));
129+
vector Uref_;
130+
Uref_.x() = spongeLayerUx;
131+
Uref_.y() = spongeLayerUy;
132+
Uref_.z() = 0.0;
133+
uniformDimensionedVectorField spongeLayerReferenceVelocity
134+
(
135+
IOobject
136+
(
137+
"spongeLayerReferenceVelocity",
138+
runTime.constant(),
139+
mesh,
140+
IOobject::NO_READ,
141+
IOobject::NO_WRITE
142+
),
143+
dimensionedVector("spongeLayerReferenceVelocity",dimensionSet(0, 1, -1, 0, 0, 0, 0),Uref_)
144+
);
145+
146+
if (spongeLayerType == "Rayleigh")
147+
{
148+
Info << spongeLayerReferenceVelocity << endl;
149+
}
105150

106151

107152

0 commit comments

Comments
 (0)