| An Example Usage of meanVelocityForce |
|
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 |
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: dev | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object fvOptions; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // momentumSource { type meanVelocityForce; active yes; meanVelocityForceCoeffs { selectionMode all; fields (U); Ubar (0.1335 0 0); relaxation 1.0; } } // ************************************************************************* // |
- fields: [Required] Name of the velocity field
- Ubar: [Required] Desired mean velocity
- relaxation: [Optional] Relaxation factor (default = 1.0)
- selectionMode: [Required] Domain the source is applied (all/cellSet/cellZone/points)
The parameter name fields has been changed from fieldNames in the latest version as described in this page.
The meanVelocityForce fvOptions calculates a momentum source so that the volume averaged velocity \eqref{eq:magUbarAve} in the whole computational domain (all) or a part of domain specified using cellSet or cellZone reaches the desired mean velocity Ubar.
\begin{align}
\frac{\sum_{i}^{} \left(\frac{\bar{\boldsymbol{u}}}{|\bar{\boldsymbol{u}}|} \cdot \boldsymbol{u}_i \right)V_i}{\sum_{i}^{} V_i}, \tag{1} \label{eq:magUbarAve}
\end{align}
where the summation is over the cells that belong to the user-specified domain, \(\bar{\boldsymbol{u}}\) is Ubar, \(\boldsymbol{u}_i\) is the velocity in the i-th cell and \(V_i\) is the volume of the i-th cell.
| patchMeanVelocityForce |
The patchMeanVelocityForce fvOptions is also available to specify the desired mean velocity on a patch instead of the volumetric average \eqref{eq:magUbarAve}.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
momentumSource { type patchMeanVelocityForce; active yes; patchMeanVelocityForceCoeffs { selectionMode all; fields (U); Ubar (0.1335 0 0); patch inout1_half0; relaxation 1.0; } } |
| Tutorial |
| Source Code |
In the case of pimpleFoam, we can find three lines related to the fvOptions as shown in the following code (UEqn.H):
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
// Solve the Momentum equation MRF.correctBoundaryVelocity(U); tmp<fvVectorMatrix> tUEqn ( fvm::ddt(U) + fvm::div(phi, U) + MRF.DDt(U) + turbulence->divDevReff(U) == fvOptions(U) ); fvVectorMatrix& UEqn = tUEqn.ref(); UEqn.relax(); fvOptions.constrain(UEqn); if (pimple.momentumPredictor()) { solve(UEqn == -fvc::grad(p)); fvOptions.correct(U); } |
In what follows, we will take a closer look at each of three highlighted lines when the meanVelocityForce fvOptions is used. The source files of this class are in src/fvOptions/sources/derived/meanVelocityForce.
- l. 11
The addSup function is called from l.11 in UEqn.H when discretizing the momentum equation and source term Su is added to the equation.
|
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 |
void Foam::fv::meanVelocityForce::addSup ( fvMatrix<vector>& eqn, const label fieldi ) { DimensionedField<vector, volMesh> Su ( IOobject ( name_ + fieldNames_[fieldi] + "Sup", mesh_.time().timeName(), mesh_, IOobject::NO_READ, IOobject::NO_WRITE ), mesh_, dimensionedVector("zero", eqn.dimensions()/dimVolume, Zero) ); scalar gradP = gradP0_ + dGradP_; UIndirectList<vector>(Su, cells_) = flowDir_*gradP; eqn += Su; } void Foam::fv::meanVelocityForce::addSup ( const volScalarField& rho, fvMatrix<vector>& eqn, const label fieldi ) { this->addSup(eqn, fieldi); } |
- l. 17
The constrain function is called from l.17 in UEqn.H. The if loop is true in the first iteration and else loop is executed for the other iterations to initialize and update the volScalarField defined as the reciprocal of the diagonal coefficients.
|
229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 |
void Foam::fv::meanVelocityForce::constrain ( fvMatrix<vector>& eqn, const label ) { if (rAPtr_.empty()) { rAPtr_.reset ( new volScalarField ( IOobject ( name_ + ":rA", mesh_.time().timeName(), mesh_, IOobject::NO_READ, IOobject::NO_WRITE ), 1.0/eqn.A() ) ); } else { rAPtr_() = 1.0/eqn.A(); } gradP0_ += dGradP_; dGradP_ = 0.0; } |
- l. 23
As the average velocity magUbarAve \eqref{eq:magUbarAve} reaches a user-specified value Ubar, the pressure gradient increment dGradP_ needed to adjust the average velocity converges to 0.
|
148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 |
void Foam::fv::meanVelocityForce::correct(volVectorField& U) { const scalarField& rAU = rAPtr_(); // Integrate flow variables over cell set scalar rAUave = 0.0; const scalarField& cv = mesh_.V(); forAll(cells_, i) { label celli = cells_[i]; scalar volCell = cv[celli]; rAUave += rAU[celli]*volCell; } // Collect across all processors reduce(rAUave, sumOp<scalar>()); // Volume averages rAUave /= V_; scalar magUbarAve = this->magUbarAve(U); // Calculate the pressure gradient increment needed to adjust the average // flow-rate to the desired value dGradP_ = relaxation_*(mag(Ubar_) - magUbarAve)/rAUave; // Apply correction to velocity field forAll(cells_, i) { label celli = cells_[i]; U[celli] += flowDir_*rAU[celli]*dGradP_; } scalar gradP = gradP0_ + dGradP_; Info<< "Pressure gradient source: uncorrected Ubar = " << magUbarAve << ", pressure gradient = " << gradP << endl; writeProps(gradP); } |
More from my site
I saw in the code it only handle pressure incompressible case, so can you guide me to apply it for compressible solver? Thank you!