MaxwellianThermal.C 3.66 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
6
     \\/     M anipulation  |
OpenFOAM bot's avatar
OpenFOAM bot committed
7
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2011-2017 OpenFOAM Foundation
9 10 11 12
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

13 14 15 16
    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.
17 18 19 20 21 22 23

    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
24
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
25 26 27 28

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

#include "MaxwellianThermal.H"
29 30 31
#include "constants.H"

using namespace Foam::constant;
32 33 34

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

Henry's avatar
Henry committed
35
template<class CloudType>
36 37 38 39 40 41
Foam::MaxwellianThermal<CloudType>::MaxwellianThermal
(
    const dictionary& dict,
    CloudType& cloud
)
:
42
    WallInteractionModel<CloudType>(cloud)
43 44 45 46 47
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Henry's avatar
Henry committed
48
template<class CloudType>
49 50 51 52 53 54
Foam::MaxwellianThermal<CloudType>::~MaxwellianThermal()
{}


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

Henry's avatar
Henry committed
55
template<class CloudType>
56 57
void Foam::MaxwellianThermal<CloudType>::correct
(
58
    typename CloudType::parcelType& p
59
)
60
{
61 62 63
    vector& U = p.U();

    scalar& Ei = p.Ei();
64

65 66
    label typeId = p.typeId();

67 68 69
    const label wppIndex = p.patch();

    const polyPatch& wpp = p.mesh().boundaryMesh()[wppIndex];
70

71
    label wppLocalFace = wpp.whichFace(p.face());
72

73
    const vector nw = p.normal();
74

75
    // Normal velocity magnitude
76
    scalar U_dot_nw = U & nw;
77 78

    // Wall tangential velocity (flow direction)
79
    vector Ut = U - U_dot_nw*nw;
80

81
    CloudType& cloud(this->owner());
82

83
    Random& rndGen = cloud.rndGen();
84 85 86 87 88 89 90 91

    while (mag(Ut) < SMALL)
    {
        // If the incident velocity is parallel to the face normal, no
        // tangential direction can be chosen.  Add a perturbation to the
        // incoming velocity and recalculate.
        U = vector
        (
92 93 94
            U.x()*(0.8 + 0.2*rndGen.sample01<scalar>()),
            U.y()*(0.8 + 0.2*rndGen.sample01<scalar>()),
            U.z()*(0.8 + 0.2*rndGen.sample01<scalar>())
95 96
        );

97
        U_dot_nw = U & nw;
98

99
        Ut = U - U_dot_nw*nw;
100 101 102 103 104 105
    }

    // Wall tangential unit vector
    vector tw1 = Ut/mag(Ut);

    // Other tangential unit vector
andy's avatar
andy committed
106
    vector tw2 = nw^tw1;
107

108
    scalar T = cloud.boundaryT().boundaryField()[wppIndex][wppLocalFace];
109

110 111
    scalar mass = cloud.constProps(typeId).mass();

112
    direction iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
113

andy's avatar
andy committed
114
    U =
115
        sqrt(physicoChemical::k.value()*T/mass)
andy's avatar
andy committed
116
       *(
117 118 119
            rndGen.GaussNormal<scalar>()*tw1
          + rndGen.GaussNormal<scalar>()*tw2
          - sqrt(-2.0*log(max(1 - rndGen.sample01<scalar>(), VSMALL)))*nw
andy's avatar
andy committed
120
        );
121

122
    U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace];
123 124

    Ei = cloud.equipartitionInternalEnergy(T, iDof);
125 126 127 128
}


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