1120 lines
28 KiB
C
1120 lines
28 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | www.openfoam.com
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
Copyright (C) 2007-2021 PCOpt/NTUA
|
|
Copyright (C) 2013-2021 FOSS GP
|
|
Copyright (C) 2019-2023 OpenCFD Ltd.
|
|
-------------------------------------------------------------------------------
|
|
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 "adjointSpalartAllmaras.H"
|
|
#include "addToRunTimeSelectionTable.H"
|
|
#include "wallDist.H"
|
|
#include "wallFvPatch.H"
|
|
#include "nutUSpaldingWallFunctionFvPatchScalarField.H"
|
|
#include "boundaryAdjointContribution.H"
|
|
#include "coupledFvPatch.H"
|
|
#include "ATCModel.H"
|
|
|
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
|
|
namespace Foam
|
|
{
|
|
namespace incompressibleAdjoint
|
|
{
|
|
namespace adjointRASModels
|
|
{
|
|
|
|
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
|
|
|
defineTypeNameAndDebug(adjointSpalartAllmaras, 0);
|
|
addToRunTimeSelectionTable
|
|
(
|
|
adjointRASModel,
|
|
adjointSpalartAllmaras,
|
|
dictionary
|
|
);
|
|
|
|
|
|
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
|
|
|
|
// * * * * * * * * * * * * Primal Spalart - Allmaras * * * * * * * * * * * * //
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::chi() const
|
|
{
|
|
return nuTilda()/nu();
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::fv1(const volScalarField& chi) const
|
|
{
|
|
const volScalarField chi3(pow3(chi));
|
|
return chi3/(chi3 + pow3(Cv1_));
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::fv2
|
|
(
|
|
const volScalarField& chi,
|
|
const volScalarField& fv1
|
|
) const
|
|
{
|
|
return 1.0 - chi/(1.0 + chi*fv1);
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::Stilda
|
|
(
|
|
const volScalarField& chi,
|
|
const volScalarField& fv1
|
|
) const
|
|
{
|
|
volScalarField Omega(::sqrt(2.0)*mag(skew(gradU_)));
|
|
|
|
return
|
|
(
|
|
max
|
|
(
|
|
Omega
|
|
+ fv2(chi, fv1)*nuTilda()/sqr(kappa_*y_),
|
|
Cs_*Omega
|
|
)
|
|
);
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::r
|
|
(
|
|
const volScalarField& Stilda
|
|
) const
|
|
{
|
|
auto tr =
|
|
tmp<volScalarField>::New
|
|
(
|
|
min
|
|
(
|
|
nuTilda()/(max(Stilda, minStilda_)*sqr(kappa_*y_)),
|
|
scalar(10)
|
|
)
|
|
);
|
|
tr.ref().boundaryFieldRef() == Zero;
|
|
|
|
return tr;
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::fw
|
|
(
|
|
const volScalarField& Stilda
|
|
) const
|
|
{
|
|
const volScalarField g(r_ + Cw2_*(pow6(r_) - r_));
|
|
|
|
return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::DnuTildaEff() const
|
|
{
|
|
return
|
|
tmp<volScalarField>::New
|
|
("DnuTildaEff", (nuTilda() + this->nu())/sigmaNut_);
|
|
}
|
|
|
|
|
|
const volScalarField& adjointSpalartAllmaras::nuTilda() const
|
|
{
|
|
return primalVars_.RASModelVariables()().TMVar1();
|
|
}
|
|
|
|
|
|
const volScalarField& adjointSpalartAllmaras::nut() const
|
|
{
|
|
return primalVars_.RASModelVariables()().nutRef();
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * Adjoint Spalart - Allmaras * * * * * * * * * * * * //
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dFv1_dChi
|
|
(
|
|
const volScalarField& chi
|
|
) const
|
|
{
|
|
volScalarField chi3(pow3(chi));
|
|
|
|
return 3.0*pow3(Cv1_)*sqr(chi/(chi3 + pow3(Cv1_)));
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dFv2_dChi
|
|
(
|
|
const volScalarField& chi,
|
|
const volScalarField& fv1,
|
|
const volScalarField& dFv1dChi
|
|
) const
|
|
{
|
|
return (chi*chi*dFv1dChi - 1.)/sqr(1. + chi*fv1);
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dStilda_dOmega
|
|
(
|
|
const volScalarField& Omega,
|
|
const volScalarField& fv2
|
|
) const
|
|
{
|
|
volScalarField fieldSwitch
|
|
(
|
|
Omega + fv2*nuTilda()/sqr(kappa_*y_) - Cs_*Omega
|
|
);
|
|
|
|
return pos(fieldSwitch) + neg(fieldSwitch)*Cs_;
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dStilda_dNuTilda
|
|
(
|
|
const volScalarField& Omega,
|
|
const volScalarField& fv2,
|
|
const volScalarField& dFv2dChi
|
|
) const
|
|
{
|
|
volScalarField invDenom(1./sqr(kappa_*y_));
|
|
volScalarField fieldSwitch(Omega + fv2*nuTilda()*invDenom - Cs_*Omega);
|
|
|
|
return pos(fieldSwitch)*(dFv2dChi*nuTilda()*invDenom/nu() + fv2*invDenom);
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dStilda_dDelta
|
|
(
|
|
const volScalarField& Omega,
|
|
const volScalarField& fv2
|
|
) const
|
|
{
|
|
volScalarField aux(fv2*nuTilda()/sqr(kappa_*y_));
|
|
volScalarField fieldSwitch(Omega + aux - Cs_*Omega);
|
|
|
|
return - 2.*pos(fieldSwitch)*aux/y_;
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dr_dNuTilda
|
|
(
|
|
const volScalarField& Stilda
|
|
) const
|
|
{
|
|
tmp<volScalarField> tdrdNutilda
|
|
(
|
|
1./(max(Stilda, minStilda_)*sqr(kappa_*y_))
|
|
*(scalar(10) - r_)/(scalar(10) - r_ + SMALL)
|
|
);
|
|
tdrdNutilda.ref().boundaryFieldRef() == Zero;
|
|
|
|
return tdrdNutilda;
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dr_dStilda
|
|
(
|
|
const volScalarField& Stilda
|
|
) const
|
|
{
|
|
tmp<volScalarField> tdrdStilda
|
|
(
|
|
- nuTilda()/sqr(max(Stilda, minStilda_)*kappa_*y_)
|
|
*(scalar(10) - r_)/(scalar(10) - r_ + SMALL)
|
|
);
|
|
tdrdStilda.ref().boundaryFieldRef() == Zero;
|
|
|
|
return tdrdStilda;
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dr_dDelta
|
|
(
|
|
const volScalarField& Stilda
|
|
) const
|
|
{
|
|
tmp<volScalarField> tdrdDelta
|
|
(
|
|
- 2.*nuTilda()/(max(Stilda, minStilda_)*sqr(kappa_*y_)*y_)
|
|
*(scalar(10) - r_)/(scalar(10) - r_ + SMALL)
|
|
);
|
|
tdrdDelta.ref().boundaryFieldRef() == Zero;
|
|
|
|
return tdrdDelta;
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dfw_dr
|
|
(
|
|
const volScalarField& Stilda
|
|
) const
|
|
{
|
|
volScalarField g(r_ + Cw2_*(pow6(r_) - r_));
|
|
|
|
dimensionedScalar pow6Cw3 = pow6(Cw3_);
|
|
volScalarField pow6g(pow6(g));
|
|
|
|
return
|
|
pow6Cw3/(pow6g + pow6Cw3)
|
|
*pow((1.0 + pow6Cw3)/(pow6g + pow6Cw3), 1.0/6.0)
|
|
*(1.0 + Cw2_*(6.0*pow5(r_) - 1.0));
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dfw_dNuTilda
|
|
(
|
|
const volScalarField& Stilda,
|
|
const volScalarField& dfwdr,
|
|
const volScalarField& dStildadNuTilda
|
|
) const
|
|
{
|
|
volScalarField invDenom(1./sqr(kappa_*y_));
|
|
|
|
return
|
|
dfwdr*(dr_dNuTilda(Stilda) + dr_dStilda(Stilda)*dStildadNuTilda);
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dfw_dOmega
|
|
(
|
|
const volScalarField& Stilda,
|
|
const volScalarField& dfwdr,
|
|
const volScalarField& dStildadOmega
|
|
) const
|
|
{
|
|
return dfwdr*dr_dStilda(Stilda)*dStildadOmega;
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dfw_dDelta
|
|
(
|
|
const volScalarField& Stilda,
|
|
const volScalarField& dfwdr,
|
|
const volScalarField& dStildadDelta
|
|
) const
|
|
{
|
|
return dfwdr*(dr_dDelta(Stilda) + dr_dStilda(Stilda)*dStildadDelta);
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dD_dNuTilda
|
|
(
|
|
const volScalarField& fw,
|
|
const volScalarField& dfwdNuTilda
|
|
) const
|
|
{
|
|
return Cw1_*(nuTilda()*dfwdNuTilda + fw)/sqr(y_);
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dP_dNuTilda
|
|
(
|
|
const volScalarField& dStildadNuTilda
|
|
) const
|
|
{
|
|
return - Cb1_*dStildadNuTilda;
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::dnut_dNuTilda
|
|
(
|
|
const volScalarField& fv1,
|
|
const volScalarField& dFv1dChi
|
|
) const
|
|
{
|
|
return dFv1dChi*nuTilda()/nu() + fv1;
|
|
}
|
|
|
|
|
|
tmp<volVectorField> adjointSpalartAllmaras::conservativeMomentumSource()
|
|
{
|
|
// Store boundary field of the conservative part,
|
|
// for use in adjoint outlet boundary conditions
|
|
forAll(mesh_.boundary(), pI)
|
|
{
|
|
const fvPatch& patch = mesh_.boundary()[pI];
|
|
if(!isA<coupledFvPatch>(patch))
|
|
{
|
|
tmp<vectorField> tnf = patch.nf();
|
|
adjMomentumBCSourcePtr_()[pI] =
|
|
(tnf & momentumSourceMult_.boundaryField()[pI])
|
|
*nuaTilda().boundaryField()[pI];
|
|
}
|
|
}
|
|
|
|
return fvc::div(momentumSourceMult_*nuaTilda());
|
|
}
|
|
|
|
|
|
void adjointSpalartAllmaras::updatePrimalRelatedFields()
|
|
{
|
|
if (changedPrimalSolution_)
|
|
{
|
|
Info<< "Updating primal-based fields of the adjoint turbulence "
|
|
<< "model ..." << endl;
|
|
|
|
// Grab references
|
|
const volVectorField& U = primalVars_.U();
|
|
|
|
// Update gradient fields
|
|
gradU_ = mask_*fvc::grad(U, "gradUStilda");
|
|
gradNuTilda_ = fvc::grad(nuTilda());
|
|
|
|
const volScalarField Omega(::sqrt(2.0)*mag(skew(gradU_)));
|
|
|
|
// Primal SA fields
|
|
volScalarField chi(this->chi());
|
|
volScalarField fv1(this->fv1(chi));
|
|
volScalarField fv2(this->fv2(chi, fv1));
|
|
Stilda_ = Stilda(chi, fv1);
|
|
r_ = r(Stilda_);
|
|
fw_ = this->fw(Stilda_);
|
|
|
|
// Derivatives of primal fields wrt to nuTilda
|
|
volScalarField dFv1_dChi(this->dFv1_dChi(chi));
|
|
volScalarField dFv2_dChi(this->dFv2_dChi(chi, fv1, dFv1_dChi));
|
|
volScalarField dStilda_dNuTilda
|
|
(this->dStilda_dNuTilda(Omega, fv2, dFv2_dChi));
|
|
volScalarField dfw_dr(this->dfw_dr(Stilda_));
|
|
volScalarField dfw_dNuTilda
|
|
(this->dfw_dNuTilda(Stilda_, dfw_dr, dStilda_dNuTilda));
|
|
|
|
// Fields to be used in the nuaTilda equation
|
|
symmAdjointProductionU_ =
|
|
symm(mask_*fvc::grad(U, "adjointProductionU"));
|
|
|
|
productionDestructionSource_ =
|
|
nuTilda()
|
|
*(
|
|
dD_dNuTilda(fw_, dfw_dNuTilda)
|
|
+ dP_dNuTilda(dStilda_dNuTilda)
|
|
);
|
|
|
|
Cdnut_ = dnut_dNuTilda(fv1, dFv1_dChi);
|
|
|
|
// Constant multiplier in the adjoint momentum source term
|
|
volScalarField dStilda_dOmega(this->dStilda_dOmega(Omega, fv2));
|
|
volScalarField dfw_dOmega
|
|
(this->dfw_dOmega(Stilda_, dfw_dr, dStilda_dOmega));
|
|
|
|
momentumSourceMult_ =
|
|
2.*skew(gradU_)
|
|
/(Omega + dimensionedScalar("SMALL", Omega.dimensions(), SMALL))
|
|
*(
|
|
- Cb1_*nuTilda()*dStilda_dOmega
|
|
+ Cw1_*sqr(nuTilda()/y_)*dfw_dOmega
|
|
);
|
|
|
|
// Set changedPrimalSolution_ to false to avoid recomputing these
|
|
// fields unless the primal has changed
|
|
changedPrimalSolution_ = false;
|
|
}
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::allocateMask()
|
|
{
|
|
tmp<volScalarField> mask;
|
|
if (limitAdjointProduction_)
|
|
{
|
|
mask = ATCModel::createLimiter(mesh_, coeffDict_);
|
|
}
|
|
else
|
|
{
|
|
mask = tmp<volScalarField>
|
|
(
|
|
new volScalarField
|
|
(
|
|
IOobject
|
|
(
|
|
"unitMask",
|
|
mesh_.time().timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh_,
|
|
dimensionedScalar("unit", dimless, scalar(1))
|
|
)
|
|
);
|
|
}
|
|
|
|
return mask;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
|
|
|
adjointSpalartAllmaras::adjointSpalartAllmaras
|
|
(
|
|
incompressibleVars& primalVars,
|
|
incompressibleAdjointMeanFlowVars& adjointVars,
|
|
objectiveManager& objManager,
|
|
const word& adjointTurbulenceModelName,
|
|
const word& modelName
|
|
)
|
|
:
|
|
adjointRASModel
|
|
(
|
|
modelName,
|
|
primalVars,
|
|
adjointVars,
|
|
objManager,
|
|
adjointTurbulenceModelName
|
|
),
|
|
|
|
sigmaNut_
|
|
(
|
|
dimensioned<scalar>::getOrAddToDict
|
|
(
|
|
"sigmaNut",
|
|
this->coeffDict_,
|
|
0.66666
|
|
)
|
|
),
|
|
kappa_
|
|
(
|
|
dimensioned<scalar>::getOrAddToDict
|
|
(
|
|
"kappa",
|
|
this->coeffDict_,
|
|
0.41
|
|
)
|
|
),
|
|
|
|
Cb1_
|
|
(
|
|
dimensioned<scalar>::getOrAddToDict
|
|
(
|
|
"Cb1",
|
|
this->coeffDict_,
|
|
0.1355
|
|
)
|
|
),
|
|
Cb2_
|
|
(
|
|
dimensioned<scalar>::getOrAddToDict
|
|
(
|
|
"Cb2",
|
|
this->coeffDict_,
|
|
0.622
|
|
)
|
|
),
|
|
Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
|
|
Cw2_
|
|
(
|
|
dimensioned<scalar>::getOrAddToDict
|
|
(
|
|
"Cw2",
|
|
this->coeffDict_,
|
|
0.3
|
|
)
|
|
),
|
|
Cw3_
|
|
(
|
|
dimensioned<scalar>::getOrAddToDict
|
|
(
|
|
"Cw3",
|
|
this->coeffDict_,
|
|
2.0
|
|
)
|
|
),
|
|
Cv1_
|
|
(
|
|
dimensioned<scalar>::getOrAddToDict
|
|
(
|
|
"Cv1",
|
|
this->coeffDict_,
|
|
7.1
|
|
)
|
|
),
|
|
Cs_
|
|
(
|
|
dimensioned<scalar>::getOrAddToDict
|
|
(
|
|
"Cs",
|
|
this->coeffDict_,
|
|
0.3
|
|
)
|
|
),
|
|
|
|
limitAdjointProduction_
|
|
(
|
|
coeffDict_.getOrDefault("limitAdjointProduction", true)
|
|
),
|
|
|
|
y_(primalVars_.RASModelVariables()().d()),
|
|
|
|
mask_(allocateMask()),
|
|
|
|
symmAdjointProductionU_
|
|
(
|
|
IOobject
|
|
(
|
|
"symmAdjointProductionU",
|
|
runTime_.timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh_,
|
|
dimensionedSymmTensor(dimless/dimTime, Zero)
|
|
),
|
|
|
|
productionDestructionSource_
|
|
(
|
|
IOobject
|
|
(
|
|
"productionDestructionSource",
|
|
runTime_.timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh_,
|
|
dimensionedScalar(dimless/dimTime, Zero)
|
|
),
|
|
|
|
Stilda_
|
|
(
|
|
IOobject
|
|
(
|
|
"Stilda",
|
|
runTime_.timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh_,
|
|
dimensionedScalar(dimless/dimTime, Zero)
|
|
),
|
|
|
|
r_
|
|
(
|
|
IOobject
|
|
(
|
|
"r",
|
|
runTime_.timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh_,
|
|
dimensionedScalar(dimless, Zero)
|
|
),
|
|
|
|
fw_
|
|
(
|
|
IOobject
|
|
(
|
|
"fw",
|
|
runTime_.timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh_,
|
|
dimensionedScalar(dimless, Zero)
|
|
),
|
|
|
|
Cdnut_
|
|
(
|
|
IOobject
|
|
(
|
|
"Cdnut",
|
|
runTime_.timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh_,
|
|
dimensionedScalar(dimless, Zero)
|
|
),
|
|
|
|
momentumSourceMult_
|
|
(
|
|
IOobject
|
|
(
|
|
"momentumSourceMult",
|
|
runTime_.timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh_,
|
|
dimensionedTensor(sqr(dimLength)/dimTime, Zero)
|
|
),
|
|
|
|
gradU_(fvc::grad(primalVars.U())),
|
|
gradNuTilda_(fvc::grad(nuTilda())),
|
|
minStilda_("SMALL", Stilda_.dimensions(), SMALL)
|
|
{
|
|
adjointTMVariablesBaseNames_.setSize(1);
|
|
adjointTMVariablesBaseNames_[0] = "nuaTilda";
|
|
// Read nuaTilda field and reset pointer to the first
|
|
// adjoint turbulence model variable
|
|
variablesSet::setField
|
|
(
|
|
adjointTMVariable1Ptr_,
|
|
mesh_,
|
|
"nuaTilda",
|
|
adjointVars.solverName(),
|
|
adjointVars.useSolverNameForFields()
|
|
);
|
|
|
|
setMeanFields();
|
|
|
|
// Set the includeDistance to true, to allow for the automatic solution
|
|
// of the adjoint eikonal equation when computing sensitivities
|
|
includeDistance_ = true;
|
|
|
|
// Update the primal related fields here so that functions computing
|
|
// sensitivities have the updated fields in case of continuation
|
|
updatePrimalRelatedFields();
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
tmp<volSymmTensorField> adjointSpalartAllmaras::devReff() const
|
|
{
|
|
const volVectorField& Ua = adjointVars_.UaInst();
|
|
return devReff(Ua);
|
|
}
|
|
|
|
|
|
tmp<volSymmTensorField> adjointSpalartAllmaras::devReff
|
|
(
|
|
const volVectorField& U
|
|
) const
|
|
{
|
|
return
|
|
tmp<volSymmTensorField>::New
|
|
(
|
|
IOobject
|
|
(
|
|
"devRhoReff",
|
|
runTime_.timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
-nuEff()*devTwoSymm(fvc::grad(U))
|
|
);
|
|
}
|
|
|
|
|
|
tmp<fvVectorMatrix> adjointSpalartAllmaras::divDevReff(volVectorField& Ua) const
|
|
{
|
|
tmp<volScalarField> tnuEff(nuEff());
|
|
const volScalarField& nuEff = tnuEff();
|
|
|
|
return
|
|
(
|
|
- fvm::laplacian(nuEff, Ua)
|
|
- fvc::div(nuEff*dev(fvc::grad(Ua)().T()))
|
|
);
|
|
}
|
|
|
|
|
|
tmp<volVectorField> adjointSpalartAllmaras::adjointMeanFlowSource()
|
|
{
|
|
addProfiling
|
|
(adjointSpalartAllmaras, "adjointSpalartAllmaras::addMomentumSource");
|
|
// cm formulation
|
|
//return (- nuTilda()*fvc::grad(nuaTilda() - conservativeMomentumSource());
|
|
|
|
// ncm formulation
|
|
return (nuaTilda()*gradNuTilda_ - conservativeMomentumSource());
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::nutJacobianTMVar1() const
|
|
{
|
|
volScalarField chi(this->chi());
|
|
volScalarField fv1(this->fv1(chi));
|
|
volScalarField dFv1_dChi(this->dFv1_dChi(chi));
|
|
|
|
return dnut_dNuTilda(fv1, dFv1_dChi);
|
|
}
|
|
|
|
|
|
tmp<scalarField> adjointSpalartAllmaras::diffusionCoeffVar1(label patchI) const
|
|
{
|
|
tmp<scalarField> tdiffCoeff
|
|
(
|
|
new scalarField(mesh_.boundary()[patchI].size(), Zero)
|
|
);
|
|
|
|
scalarField& diffCoeff = tdiffCoeff.ref();
|
|
|
|
diffCoeff =
|
|
(nuTilda().boundaryField()[patchI] + nu()().boundaryField()[patchI])
|
|
/sigmaNut_.value();
|
|
|
|
return tdiffCoeff;
|
|
}
|
|
|
|
|
|
const boundaryVectorField&
|
|
adjointSpalartAllmaras::adjointMomentumBCSource() const
|
|
{
|
|
// Computed in conservativeMomentumSource
|
|
return adjMomentumBCSourcePtr_();
|
|
}
|
|
|
|
|
|
const boundaryVectorField& adjointSpalartAllmaras::wallShapeSensitivities()
|
|
{
|
|
boundaryVectorField& wallShapeSens = wallShapeSensitivitiesPtr_();
|
|
|
|
forAll(mesh_.boundary(), patchI)
|
|
{
|
|
const fvPatch& patch = mesh_.boundary()[patchI];
|
|
|
|
tmp<vectorField> tnf = patch.nf();
|
|
if (isA<wallFvPatch>(patch) && patch.size() != 0)
|
|
{
|
|
wallShapeSens[patchI] =
|
|
- nuaTilda().boundaryField()[patchI].snGrad()
|
|
*diffusionCoeffVar1(patchI)
|
|
*nuTilda().boundaryField()[patchI].snGrad()*tnf;
|
|
}
|
|
}
|
|
|
|
return wallShapeSens;
|
|
}
|
|
|
|
|
|
const boundaryVectorField& adjointSpalartAllmaras::wallFloCoSensitivities()
|
|
{
|
|
boundaryVectorField& wallFloCoSens = wallFloCoSensitivitiesPtr_();
|
|
|
|
forAll(mesh_.boundary(), patchI)
|
|
{
|
|
tmp<vectorField> tnf = mesh_.boundary()[patchI].nf();
|
|
|
|
wallFloCoSens[patchI] =
|
|
nuaTilda().boundaryField()[patchI]
|
|
*nuTilda().boundaryField()[patchI]*tnf;
|
|
}
|
|
|
|
return wallFloCoSens;
|
|
}
|
|
|
|
|
|
tmp<volScalarField> adjointSpalartAllmaras::distanceSensitivities()
|
|
{
|
|
const volVectorField& U = primalVars_.U();
|
|
const volVectorField& Ua = adjointVars_.Ua();
|
|
|
|
// Primal SA fields
|
|
volScalarField chi(this->chi());
|
|
volScalarField fv1(this->fv1(chi));
|
|
volScalarField fv2(this->fv2(chi, fv1));
|
|
volScalarField Omega(::sqrt(2.0)*mag(gradU_));
|
|
|
|
// Derivatives of primal fields wrt to nuTilda
|
|
volScalarField dFv1_dChi(this->dFv1_dChi(chi));
|
|
volScalarField dFv2_dChi(this->dFv2_dChi(chi, fv1, dFv1_dChi));
|
|
volScalarField dStilda_dDelta(this->dStilda_dDelta(Omega, fv2));
|
|
volScalarField dfw_dr(this->dfw_dr(Stilda_));
|
|
volScalarField dfw_dDelta
|
|
(this->dfw_dDelta(Stilda_, dfw_dr, dStilda_dDelta));
|
|
|
|
auto tadjointEikonalSource =
|
|
tmp<volScalarField>::New
|
|
(
|
|
"adjointEikonalSource" + type(),
|
|
(
|
|
- Cb1_*nuTilda()*dStilda_dDelta
|
|
+ Cw1_*sqr(nuTilda()/y_)*(dfw_dDelta - 2.*fw_/y_)
|
|
)*nuaTilda()
|
|
);
|
|
volScalarField& adjointEikonalSource = tadjointEikonalSource.ref();
|
|
|
|
// if wall functions are used, add appropriate source terms
|
|
typedef nutUSpaldingWallFunctionFvPatchScalarField
|
|
SAwallFunctionPatchField;
|
|
|
|
const volScalarField::Boundary& nutBoundary = nut().boundaryField();
|
|
const scalarField& V = mesh_.V().field();
|
|
|
|
tmp<volScalarField> tnuEff = nuEff();
|
|
const volScalarField& nuEff = tnuEff();
|
|
|
|
forAll(nutBoundary, patchi)
|
|
{
|
|
const fvPatch& patch = mesh_.boundary()[patchi];
|
|
if
|
|
(
|
|
isA<SAwallFunctionPatchField>(nutBoundary[patchi])
|
|
&& patch.size() != 0
|
|
)
|
|
{
|
|
const scalar kappa_(0.41);
|
|
const scalar E_(9.8);
|
|
tmp<vectorField> tnf = patch.nf();
|
|
const vectorField& nf = tnf.ref();
|
|
const scalarField& magSf = patch.magSf();
|
|
|
|
const fvPatchVectorField& Up = U.boundaryField()[patchi];
|
|
const fvPatchVectorField& Uap = Ua.boundaryField()[patchi];
|
|
const vectorField Uc(Up.patchInternalField());
|
|
const vectorField Uc_t(Uc - (Uc & nf)*nf);
|
|
|
|
// By convention, tf has the direction of the tangent
|
|
// PRIMAL velocity at the first cell off the wall
|
|
const vectorField tf(Uc_t/mag(Uc_t));
|
|
|
|
const scalarField nuw(nuEff.boundaryField()[patchi]);
|
|
const scalarField nu(this->nu()().boundaryField()[patchi]);
|
|
const fvPatchScalarField& yC = y()[patchi];
|
|
|
|
const scalarField magGradU(mag(Up.snGrad()));
|
|
|
|
// Note: What happens in separation?? sign change needed
|
|
const scalarField vtau(sqrt(nuw*magGradU));
|
|
|
|
// Note: mag for positive uPlus
|
|
const scalarField uPlus(mag(Uc)/vtau);
|
|
|
|
const scalarField yPlus(yC*vtau/nu);
|
|
const scalarField kUu(min(kappa_*uPlus, scalar(50)));
|
|
const scalarField auxA
|
|
((kappa_/E_)*(exp(kUu) - 1 - kUu - 0.5*kUu*kUu));
|
|
const scalarField Cwf_d(sqr(vtau)/nu/(yPlus+uPlus*(1 + auxA)));
|
|
|
|
// Tangential components are according to tf
|
|
autoPtr<boundaryAdjointContribution> boundaryContrPtr
|
|
(
|
|
boundaryAdjointContribution::New
|
|
(
|
|
"objectiveManager" + objectiveManager_.adjointSolverName(),
|
|
objectiveManager_.adjointSolverName(),
|
|
"incompressible",
|
|
patch
|
|
)
|
|
);
|
|
tmp<vectorField> tsource(boundaryContrPtr->normalVelocitySource());
|
|
|
|
const scalarField rt(tsource() & tf);
|
|
const scalarField Uap_t(Uap & tf);
|
|
|
|
const labelList& faceCells = patch.faceCells();
|
|
forAll(faceCells, faceI)
|
|
{
|
|
label cellI = faceCells[faceI];
|
|
adjointEikonalSource[cellI] -=
|
|
2.*( rt[faceI] + Uap_t[faceI] )
|
|
*vtau[faceI]*Cwf_d[faceI]*magSf[faceI]
|
|
/V[cellI]; // Divide with cell volume since the term
|
|
// will be used as a source term in the
|
|
// adjoint eikonal equation
|
|
}
|
|
}
|
|
}
|
|
|
|
return tadjointEikonalSource;
|
|
}
|
|
|
|
|
|
tmp<volTensorField> adjointSpalartAllmaras::FISensitivityTerm()
|
|
{
|
|
const volVectorField& U = primalVars_.U();
|
|
|
|
tmp<volTensorField> tgradU = fvc::grad(U);
|
|
const volTensorField& gradU = tgradU.cref();
|
|
tmp<volVectorField> tgradNuTilda = fvc::grad(nuTilda());
|
|
volVectorField& gradNuTilda = tgradNuTilda.ref();
|
|
tmp<volVectorField> tgradNuaTilda = fvc::grad(nuaTilda());
|
|
volVectorField::Boundary& gradNuaTildab =
|
|
tgradNuaTilda.ref().boundaryFieldRef();
|
|
|
|
// Explicitly correct the boundary gradient to get rid of
|
|
// the tangential component
|
|
forAll(mesh_.boundary(), patchI)
|
|
{
|
|
const fvPatch& patch = mesh_.boundary()[patchI];
|
|
if (isA<wallFvPatch>(patch))
|
|
{
|
|
tmp<vectorField> tnf = patch.nf();
|
|
const vectorField& nf = tnf();
|
|
// gradU:: can cause problems in zeroGradient patches for U
|
|
// and zero fixedValue for nuTilda.
|
|
// S becomes 0 and is used as a denominator in G
|
|
//gradU.boundaryField()[patchI] =
|
|
// nf * U_.boundaryField()[patchI].snGrad();
|
|
gradNuTilda.boundaryFieldRef()[patchI] =
|
|
nf*nuTilda().boundaryField()[patchI].snGrad();
|
|
gradNuaTildab[patchI] =
|
|
nf*nuaTilda().boundaryField()[patchI].snGrad();
|
|
}
|
|
}
|
|
|
|
// delta vorticity
|
|
volScalarField Omega(::sqrt(2.0)*mag(skew(gradU)));
|
|
volTensorField deltaOmega
|
|
(
|
|
(
|
|
(gradU & gradU)().T() //jk
|
|
- (gradU & gradU.T()) //symmetric
|
|
)
|
|
/(Omega + dimensionedScalar("SMALL", Omega.dimensions(), SMALL))
|
|
);
|
|
tgradU.clear();
|
|
|
|
volScalarField chi(this->chi());
|
|
volScalarField fv1(this->fv1(chi));
|
|
volScalarField fv2(this->fv2(chi, fv1));
|
|
|
|
volScalarField dfw_dr(this->dfw_dr(Stilda_));
|
|
volScalarField dStilda_dOmega(this->dStilda_dOmega(Omega, fv2));
|
|
volScalarField dfw_dOmega
|
|
(this->dfw_dOmega(Stilda_, dfw_dr, dStilda_dOmega));
|
|
|
|
return
|
|
tmp<volTensorField>::New
|
|
(
|
|
"volSensTerm",
|
|
// jk, cm formulation for the TM model convection
|
|
- (nuaTilda()*(U*gradNuTilda))
|
|
// jk, symmetric in theory
|
|
+ nuaTilda()*fvc::grad(DnuTildaEff()*gradNuTilda)().T()
|
|
// jk
|
|
- DnuTildaEff()*(tgradNuaTilda*gradNuTilda)
|
|
// symmetric
|
|
+ 2.*nuaTilda()*Cb2_/sigmaNut_*(gradNuTilda*gradNuTilda)
|
|
+ (
|
|
- Cb1_*nuTilda()*dStilda_dOmega
|
|
+ Cw1_*sqr(nuTilda()/y_)*dfw_dOmega
|
|
)
|
|
*nuaTilda()*deltaOmega // jk
|
|
);
|
|
}
|
|
|
|
|
|
void adjointSpalartAllmaras::nullify()
|
|
{
|
|
variablesSet::nullifyField(nuaTilda());
|
|
}
|
|
|
|
|
|
void adjointSpalartAllmaras::correct()
|
|
{
|
|
addProfiling
|
|
(adjointSpalartAllmaras, "adjointSpalartAllmaras::correct");
|
|
if (!adjointTurbulence_)
|
|
{
|
|
return;
|
|
}
|
|
|
|
adjointTurbulenceModel::correct();
|
|
|
|
updatePrimalRelatedFields();
|
|
|
|
const surfaceScalarField& phi = primalVars_.phi();
|
|
const volVectorField& Ua = adjointVars_.UaInst();
|
|
|
|
volScalarField gradNua(gradNuTilda_ & fvc::grad(nuaTilda()));
|
|
volScalarField gradUaR
|
|
(
|
|
2.0*fvc::grad(Ua) && symmAdjointProductionU_
|
|
);
|
|
|
|
dimensionedScalar oneOverSigmaNut = 1./sigmaNut_;
|
|
|
|
nuaTilda().storePrevIter();
|
|
|
|
tmp<fvScalarMatrix> nuaTildaEqn
|
|
(
|
|
fvm::ddt(nuaTilda())
|
|
+ fvm::div(-phi, nuaTilda())
|
|
- fvm::laplacian(DnuTildaEff(), nuaTilda())
|
|
// Note: Susp
|
|
+ fvm::SuSp(productionDestructionSource_, nuaTilda())
|
|
+ fvc::laplacian(2.0*Cb2_*oneOverSigmaNut*nuaTilda(), nuTilda())
|
|
+ gradNua*oneOverSigmaNut
|
|
==
|
|
// always a negative contribution to the lhs. No Sp used!
|
|
Cb1_*Stilda_*nuaTilda()
|
|
//always a positive contribution to the lhs. no need for SuSp
|
|
- fvm::Sp(Cw1_*fw_*nuTilda()/sqr(y_), nuaTilda())
|
|
- Cdnut_*gradUaR
|
|
);
|
|
|
|
// Add sources from the objective functions
|
|
objectiveManager_.addTMEqn1Source(nuaTildaEqn.ref());
|
|
|
|
nuaTildaEqn.ref().relax();
|
|
solve(nuaTildaEqn);
|
|
nuaTilda().correctBoundaryConditions();
|
|
nuaTilda().relax();
|
|
|
|
if (adjointVars_.getSolverControl().printMaxMags())
|
|
{
|
|
scalar maxDeltaNuaTilda =
|
|
gMax(mag(nuaTilda() - nuaTilda().prevIter())());
|
|
dimensionedScalar maxNuaTilda = max(mag(nuaTilda()));
|
|
Info<< "Max mag of nuaTilda = " << maxNuaTilda.value() << endl;
|
|
Info<< "Max mag of delta nuaTilda = " << maxDeltaNuaTilda << endl;
|
|
}
|
|
}
|
|
|
|
|
|
bool adjointSpalartAllmaras::read()
|
|
{
|
|
if (adjointRASModel::read())
|
|
{
|
|
sigmaNut_.readIfPresent(this->coeffDict());
|
|
kappa_.readIfPresent(this->coeffDict());
|
|
|
|
Cb1_.readIfPresent(this->coeffDict());
|
|
Cb2_.readIfPresent(this->coeffDict());
|
|
Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
|
|
Cw2_.readIfPresent(this->coeffDict());
|
|
Cw3_.readIfPresent(this->coeffDict());
|
|
Cv1_.readIfPresent(this->coeffDict());
|
|
Cs_.readIfPresent(this->coeffDict());
|
|
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
|
|
} // End namespace adjointRASModels
|
|
} // End namespace incompressibleAdjoint
|
|
} // End namespace Foam
|
|
|
|
// ************************************************************************* //
|