Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

This commit is contained in:
mattijs
2010-09-30 14:13:54 +01:00
410 changed files with 32063 additions and 35332 deletions

209
README
View File

@ -1,209 +0,0 @@
# -*- mode: org; -*-
#
#+TITLE: OpenFOAM README for version 1.6
#+AUTHOR: OpenCFD Ltd.
#+DATE: April 2010
#+LINK: http://www.opencfd.co.uk
#+OPTIONS: author:nil ^:{}
* Copyright
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. See the file COPYING in this directory, for a description of the
GNU General Public License terms under which you can copy the files.
* System requirements
OpenFOAM is developed and tested on Linux, but should work with other POSIX
systems. To check your system setup, execute the foamSystemCheck script in
the bin/ directory of the OpenFOAM installation. If no problems are reported,
proceed to "3. Installation"; otherwise contact your system administrator.
If the user wishes to run OpenFOAM in 32/64-bit mode they should consult the
section "Running OpenFOAM in 32-bit mode".
*** Qt (from http://trolltech.com/products/qt)
The ParaView 3.7.0 visualisation package requires Qt to be installed on the
system. ParaView's producers state that ParaView is only officially
supported on Qt version 4.6.x. However, we have found in limited tests that
ParaView works satisfactorily with Qt than 4.5.x. To
check whether Qt4 is installed, and the version, type:
+ qmake --version
Both 32-bit and 64-bit version of ParaView were compiled with Qt-4.4.3 (with
openSUSE-11.1). If the user finds that a ParaView binary fails to run, then
it is almost certainly due to a conflict in compiled and installed Qt
versions and they will need to consult the section below on "Compiling
ParaView and the PV3FoamReader module."
The default versions of Qt used by some GNU/Linux releases are as follows.
+ ubuntu-7.10: Version 4.3.2
+ ubuntu-8.04: Version 4.3.4
+ ubuntu-9.04: Version 4.5.0
+ openSUSE-10.2: Version 4.2.1 - too old
+ openSUSE-10.3: Version 4.3.1
+ openSUSE-11.0: Version 4.4.0
+ openSUSE-11.1: Version 4.4.3
+ openSUSE-11.2: Version 4.5.3
Compilation and running of ParaView has been successful using the libraries
downloaded in the "libqt4-dev" package on ubuntu.
If you don't have an appropriate version of Qt installed you can download
the sources e.g.:
http://get.qt.nokia.com/qt/source/qt-everywhere-opensource-src-4.6.2.tar.gz
and compile and install in /usr/local or some other location that does not
conflict with the pre-installed version.
* Installation
Download and unpack the files in the $HOME/OpenFOAM directory as described in:
http://www.OpenFOAM.org/download.html
The environment variable settings are contained in files in an etc/ directory
in the OpenFOAM release. e.g. in
+ $HOME/OpenFOAM/OpenFOAM-1.6/etc/
1) EITHER, if running bash or ksh (if in doubt type 'echo $SHELL'), source the
etc/bashrc file by adding the following line to the end of your
$HOME/.bashrc file:
+ . $HOME/OpenFOAM/OpenFOAM-1.6/etc/bashrc
Then update the environment variables by sourcing the $HOME/.bashrc file by
typing in the terminal:
+ . $HOME/.bashrc
2) OR, if running tcsh or csh, source the etc/cshrc file by adding the
following line to the end of your $HOME/.cshrc file:
+ source $HOME/OpenFOAM/OpenFOAM-1.6/etc/cshrc
Then update the environment variables by sourcing the $HOME/.cshrc file by
typing in the terminal:
+ source $HOME/.cshrc
*** Installation in alternative locations
OpenFOAM may also be installed in alternative locations. However, the
installation directory should be network available (e.g., NFS) if parallel
calculations are planned.
The environment variable 'FOAM_INST_DIR' can be used to find and source the
appropriate resource file. Here is a bash/ksh/sh example:
+ export FOAM_INST_DIR=/data/app/OpenFOAM
+ foamDotFile=$FOAM_INST_DIR/OpenFOAM-1.6/etc/bashrc
+ [ -f $foamDotFile ] && . $foamDotFile
and a csh/tcsh example:
+ setenv FOAM_INST_DIR /data/app/OpenFOAM
+ foamDotFile=$FOAM_INST_DIR/OpenFOAM-1.6/etc/cshrc
+ if ( -f $foamDotFile ) source $foamDotFile
The value set in '$FOAM_INST_DIR' will be used to locate the remaining parts
of the OpenFOAM installation.
* Building from Sources (Optional)
If you cannot find an appropriate binary pack for your platform, you can build
the complete OpenFOAM from the source-pack. You will first need to compile or
obtain a recent version of gcc (we recommend gcc-4.4.?) for your platform,
which may be obtained from http://gcc.gnu.org/.
Install the compiler in
$WM_THIRD_PARTY_DIR/platforms/$WM_ARCH$WM_COMPILER_ARCH/gcc-<GCC_VERSION>
and change the gcc version number in $WM_PROJECT_DIR/etc/settings.sh and
$WM_PROJECT_DIR/etc/settings.csh appropriately and finally update the
environment variables as in section 3.
Now go to the top-level source directory $WM_PROJECT_DIR and execute the
top-level build script './Allwmake'. In principle this will build everything,
but if problems occur with the build order it may be necessary to update the
environment variables and re-execute './Allwmake'.
If you experience difficulties with building the source-pack, or your platform
is not currently supported, please contact <enquiries@OpenCFD.co.uk> to
negotiate a support contract and we will do the port and maintain it for
future releases.
* Testing the installation
To check your installation setup, execute the 'foamInstallationTest' script
(in the bin/ directory of the OpenFOAM installation). If no problems are
reported, proceed to getting started with OpenFOAM; otherwise, go back and
check you have installed the software correctly and/or contact your system
administrator.
* Getting Started
Create a project directory within the $HOME/OpenFOAM directory named
<USER>-1.6 (e.g. 'chris-1.6' for user chris and OpenFOAM version 1.6)
and create a directory named 'run' within it, e.g. by typing:
+ mkdir -p $FOAM_RUN/run
Copy the 'tutorial' examples directory in the OpenFOAM distribution to the
'run' directory. If the OpenFOAM environment variables are set correctly,
then the following command will be correct:
+ cp -r $WM_PROJECT_DIR/tutorials $FOAM_RUN
Run the first example case of incompressible laminar flow in a cavity:
+ cd $FOAM_RUN/tutorials/incompressible/icoFoam/cavity
+ blockMesh
+ icoFoam
+ paraFoam
Refer to the OpenFOAM User Guide at http://www.OpenFOAM.org/doc/user.html for
more information.
* Compiling Paraview 3.7.0 and the PV3FoamReader module
If there are problems encountered with ParaView, then it may be necessary to
compile ParaView from sources. The compilation
is a fairly simple process using the makeParaView script
(found in ThirdParty directory), which has worked in our tests with other
packages supplied in the ThirdParty directory, namely cmake-2.8.0 and
gcc-4.4.3. Execute the following:
+ cd $WM_THIRD_PARTY_DIR
+ rm -rf paraview-3.7.0/platforms
+ rm -rf platforms/*/paraview-3.7.0
+ ./makeParaView
The PV3blockMeshReader and the PV3FoamReader ParaView plugins are compiled
as usual for OpenFOAM utilities:
+ cd $FOAM_UTILITIES/postProcessing/graphics/PV3Readers/
+ ./Allwclean
+ ./Allwmake
*** Compiling Paraview with a local version of Qt
If the user still encounters problems with ParaView, it may relate to the
version of Qt, in which case, it is recommended that the user first
downloads a supported version of Qt /e.g./ 4.5.3 as described in the section
on "Qt". The user should unpack the source pack in the $WM_THIRD_PARTY_DIR.
Then the user can build Qt by executing from within $WM_THIRD_PARTY_DIR:
+ ./makeQt
The user should then compile ParaView using the local version of Qt by
executing makeParaView with the -qmake option, giving the full path of the
newly built qmake as an argument:
+ ./makeParaView -qmake <path_to_qmake>
The user must then recompile the PV3blockMeshReader and the
PV3FoamReader plugins as usual (see above).
* Documentation
http://www.OpenFOAM.org/doc
* Help
http://www.OpenFOAM.org http://www.OpenFOAM.org/discussion.html
* Reporting Bugs in OpenFOAM
http://www.OpenFOAM.org/bugs.html
* Running OpenFOAM in 32-bit mode on 64-bit machines
Linux users with a 64-bit machine may install either the OpenFOAM 32-bit
version (linux) or the OpenFOAM 64-bit version (linux64), or both. The 64-bit
is the default mode on a 64-bit machine. To use an installed 32-bit version,
the user must set the environment variable WM_ARCH_OPTION to 32 before
sourcing the etc/bashrc (or etc/cshrc) file.

View File

@ -1,374 +0,0 @@
# -*- mode: org; -*-
#
#+TITLE: OpenFOAM release notes for version 1.6
#+AUTHOR: OpenCFD Ltd.
#+DATE: July 2009
#+LINK: http://www.opencfd.co.uk
#+OPTIONS: author:nil ^:{}
* Overview
OpenFOAM-1.6 is a significant upgrade to version 1.5 in ways that are
outlined below. This release passes all our standard tests and the
tutorials have been broadly checked. If there are any bugs, please report
them using the instructions set out here:
http://www.OpenFOAM.org/bugs.html.
* GNU/Linux version
The 32bit and 64bit binary packs of the OpenFOAM release were compiled on
a machine running openSUSE GNU/Linux version 11.1 and also tested on
Ubuntu 9. We recommend that users run OpenFOAM on one of these, or on a
similarly recent version of GNU/Linux. This release has also been
successfully compiled and tested on older GNU/Linux releases, but this
requires the installation of Qt 4.3.? (the sources for which are supplied
with OpenFOAM-1.6, see README) for ParaView-3 to run.
* C++ Compiler version
+ Release compiled with GCC 4.3.3.
+ Built-in support for the Intel C++ 10.? compiler (untested).
+ The choice of the compiler is controlled by the setting of the
~$WM_COMPILER~ and ~$WM_COMPILER_ARCH~ environment variables in the
/OpenFOAM-1.6/etc/bashrc/ (or /cshrc/) file.
+ The location of the compiler installation is controlled by the
~$compilerInstall~ environment variable in the
/OpenFOAM-1.6/etc/settings.sh/ (or /settings.csh/) file.
* Library developments
*** Core library
***** Dictionary improvements/changes
+ Dictionaries can use words (unquoted) or regular expressions (quoted)
for their keywords. When searching, an exact match has priority over a
regular expression match. Multiple regular expressions are matched in
reverse order.
+ The *new* =#includeIfPresent= directive is similar to the =#include=
directive, but does not generate an error if the file does not exist.
+ The default =#inputMode= is now '=merge=', which corresponds to the most
general usage. The =#inputMode warn= corresponds to the previous default
behaviour.
+ The *new* =#inputMode protect= can be used to conditionally merge
default values into existing dictionaries.
+ *New* =digest()= method to calculate and return the SHA1 message digest.
***** Regular Expressions
The addition of regular expressions marks a major improvement in
usability.
+ *New* =regExp= class provides support for accessing POSIX extended
regular expresssions from within OpenFOAM.
+ *New* =wordRe= class can contain a =word= or a =regExp= .
+ *New* =stringListOps= to search string lists based on regular
expressions, =wordRe= or =wordReList=.
+ =Istream= and =Ostream= now retain backslashes when reading/writing
strings.
***** Convenience changes
+ =IOobject= has a *new* constructor for creating an =IOobject= from a
single-path specification (eg, see =blockMesh -dict= option).
+ =argList= has *new* convenience methods for accessing options more
directly: =option()=, =optionFound()=, =optionLookup()=, =optionRead()=,
=optionReadIfPresent()=.
+ The *new* =readList(Istream&)= can read a bracket-delimited list or
handle a single value as a list of size 1. This can be a useful
convenience when processing command-line options.
+ Export *new* environment variable =FOAM_CASENAME= that contains the
name part of the =FOAM_CASE= environment variable.
*** Turbulence modelling
+ Major development of turbulence model libraries to give extra flexibility
at the solver level. For solvers that can support either RAS/LES
computations, the selection is made in the
/constant/turbulenceProperties/, by setting the =simulationType= keyword
to:
- =laminar=,
- =RASModel=,
- =LESModel=.
+ Depending on the selection, the model is the instantiated from /constant//
- /RASProperties/,
- /LESProperties/.
***** RAS wall functions
Wall functions are now run-time selectable per patch for RAS.
+ Velocity:
- Apply to turbulent viscosities =nut= or =mut=,
- Apply to =k=, =Q=, =R=,
- Apply to =epsilon=, =omega=.
+ Temperature:
- Apply to turbulent thermal diffusivity, =alphat= (compressible only).
+ To apply wall functions:
- To recapture the functionality of previous OpenFOAM versions (v1.5 and
earlier) assign:
- for velocity:
- =nut=: =nutWallFunction=,
- =mut=: =muWallFunction=,
- =epsilon=: =epsilonWallFunction=,
- =omega=: =omegaWallFunction=,
- =k=, =q=, =R=: =kqRWallFunction=.
- for temperature:
- =alphat=: =alphatWallFunction=.
- New =alphaSgsJayatillekeWallFunction= thermal wall function for
compressible LES.
***** *New* LES turbulence models
+ Spalart-Allmaras DDES.
+ Spalart-Allmaras IDDES.
***** Upgrading:
+ *New* utility - =applyWallFunctionBoundaryConditions=.
+ Solvers will automatically update existing cases.
- New fields created based on the presence of the =nut/mut= field.
- Boundary conditions include scoping, i.e compressibility:: for
compressible solvers.
- Modified fields will be backed-up to /<field>.old/.
+ NOTE:
- Fields are only updated for those fields associated with the current
turbulence model selection, i.e. if fields exist for use with other
models, they will not be updated.
- The new specification is not backwards compatible.
*** Thermo-physical Models
+ Old compressibility-based thermo package renamed
=basicThermo= \rightarrow =basicPsiThermo=.
+ *New* =basicRhoThermo= thermo package.
- Additional density field stored.
- General form - can be used for other types of media, e.g. liquids.
- Additional polynomial-based thermodynamics:
- Equation of state: =icoPolynomial=,
- Transport: =polynomialTransport=,
- Thermo: =hPolynomialThermo=.
+ Removed earlier hard-coding of gas thermophysics for chemistry modelling:
- =reactingMixture= now templated on thermo package,
- =chemistryModel= now templated on thermo package,
- =chemistrySolver= now templated on thermo package.
+ *New* =fvDOM= radition model
- finite volume, discrete ordinates method.
+ *New* (reinstated) =eThermo= thermodynamics package
- internal energy-based thermodynamics.
*** Lagrangian
***** Intermediate
+ Overhaul of the underlying framework.
+ Reacting now split into reacting and reacting multiphase.
+ New structure for variable composition.
+ Many new sub-models, including:
- Injection
- =PatchInjection= - injection local to patch face cells,
- =FieldActivatedInjection= - injection based on satisfying external
criterion,
- LookupTableInjection - explicity define injection locations and all
parcel properties.
- Post-processing
- patch post-processing - collect data for parcels impacting user,
defined patches.
- Patch interaction
- generalised behaviour for parcel interaction with patch.
- Phase change
- liquid evaporation.
***** Coal combustion
+ *New* library - extension of reacting-multiphase functionality.
- Surface reaction/combustion models.
*** Discrete methods
+ *New* library offering DSMC simulation functionality - see =dsmcFoam=
below.
+ Significant development of the libraries offering molecular dynamics
simulation functionality - see =mdFoam= and =mdEquilibrationFoam= below.
*** Numerics
+ *new* polynomial-fit higher-order interpolation schemes:
- =biLinearFit=
- =linearFit=
- =quadraticLinearFit=
- =quadraticFit=
- =linearPureUpwindFit=
- =quadraticLinearPureUpwindFit=
- =quadraticLinearUpwindFit=
- =quadraticUpwindFit=
- =cubicUpwindFit=
+ *new* polynomial-fit higher-order Sn-Grad: =quadraticFitSnGrad=.
*** *New* surfMesh library
Provides a more efficient storage mechanism than possible with =triSurface=
without restrictions on the shape of the face (templated parameter).
+ =MeshedSurface= class - with zero or more contiguous =surfZones= .
+ =UnsortedMeshedSurface= class - unordered surface zones (as per
=triSurface=).
+ =surfMesh= class - for reading/writing in native OpenFOAM format.
* Solvers
*** Solver restructuring
The upgrade to the turbulence models means that the simulation type, i.e.
laminar, RAS or LES can be selected at run time. This has allowed a reduction
in the number of solvers, simplifying the overall code structure
+ Solvers which support laminar, RAS and LES:
- =turbFoam=, =oodles= \rightarrow =pisoFoam=.
- =turbDyMFoam= \rightarrow =pimpleDyMFoam=.
- =rhoTurbFoam=, =coodles= \rightarrow =rhoPisoFoam=.
- =xoodles= \rightarrow absorbed into =XiFoam=.
- =buoyantFoam=, =lesBuoyantFoam= \rightarrow =buoyantPisoFoam=.
- =interFoam=, =rasInterFoam=, =lesInterFoam= \rightarrow =interFoam=.
- =lesCavitatingFoam=, =rasCavitatingFoam= \rightarrow =cavitatingFoam=.
+ Solvers which support LES only:
- =channelOodles= \rightarrow =channelFoam= (LES).
+ =pd= replaced by static pressure =p=. All solvers in which buoyancy affects
might be strong have been converted from using =pd= to =p= with improved
numerics to give equally good accuracy and stability. This change is
prompted by the need to remove the confusion surrounding the meaning and
purpose of =pd=.
+ =g= (acceleration due to gravity) is now a *new*
=uniformDimensionedVectorField= which has the behaviour of a field, is
registered to an =objectRegistry=, but stores only a single value. Thus
=g= and other =UniformDimensionedFields= can be created and looked-up
elsewhere in the code, /e.g./ in =fvPatchFields=.
*** Solver control improvements
Now uses consistent dictionary entries for the solver controls.
+ This Allows dictionary substitutions and regular expressions in
/system/fvSolution/.
+ The old solver control syntax is still supported (warning emitted), but
the *new* =foamUpgradeFvSolution= utility can be used to convert
/system/fvSolution/ to the new format.
*** *New* Solvers
+ =buoyantBoussinesqSimpleFoam= Steady state heat transfer solver using a
Boussinesq approximation for buoyancy, with laminar, RAS or LES turbulence
modelling.
+ =buoyantBoussinesqPisoFoam= Transient heat transfer solver using a
Boussinesq approximation for buoyancy, with laminar, RAS or LES turbulence
modelling.
+ =coalChemistryFoam= Transient, reacting lagrangian solver, employing a coal
cloud and a thermo cloud, with chemistry, and laminar, RAS or LES turbulence
modelling.
+ =porousExplicitSourceReactingParcelFoam= Transient, reacting lagrangian
solver, employing a single phase reacting cloud, with porous media, explicit
mass sources, and laminar, RAS or LES turbulence modelling.
+ =rhoReactingFoam= Density-based thermodynamics variant of the reactingFoam
solver, i.e. now applicable to liquid systems.
+ =dsmcFoam= DSMC (Direct Simulation Monte-Carlo) solver for rarefied gas
dynamics simulations, able to simulate mixtures of an arbitrary number of
gas species. The variable hard sphere collision model with Larsen-Borgnakke
internal energy redistribution (see "Molecular Gas Dynamics and the Direct
Simulation of Gas Flows" G.A. Bird, 1994) is available; other run-time
selectable collision models can be easily added.
*** Updated solvers
+ =mdFoam= Molecular Dynamics (MD) solver able to simulate a mixture of an
arbitrary number of mono-atomic and small, rigid polyatomic (i.e. H2O, N2)
molecular species, with 6 degree of freedom motion, in complex geometries. A
molecule of any species can be built by specifying its sites of mass and
charge. All molecules interact with short-range dispersion forces and
pairwise electrostatic interactions using methods described in: Fennell and
Gezelter, J. Chem. Phys. 124, 234104 (2006).
+ =mdEquilibrationFoam= Similar to mdFoam, but employs velocity scaling to
adjust the simulation temperature to a target value. Useful to equilibrate a
case before simulation.
+ =chtMultiRegionFoam= New boundary condition allows independent decomposition
of coupled regions without any constraint on the decomposition.
* Boundary conditions
+ Improved set of direct mapped boundary conditions.
+ =buoyantPressureFvPatchScalarField=, the *new* buoyancy pressure boundary
condition now supports =p= and =pd= for backward compatibility.
+ =uniformDensityHydrostaticPressure= is an additional pressure boundary
condition to aid the transition from =pd= to =p= as it behaves similarly to
specifying a uniform =pd= at an outlet for example.
+ =activeBaffleVelocity= dynamically combines cyclic and wall patches so that
the flow through the patch can be controlled /e.g./ by pressure drop.
+ =rotatingWallVelocity= specifies a rotating velocity, given the rotational
speed, origin and axis.
* Utilities
*** Improvements
+ =blockMesh= has a *new* =-dict= option for specifying an alternative
dictionary for the block mesh description. The '=convertToMeters=' entry
is now optional, and the alternative '=scale=' entry can be used for
less typing.
+ =foamToEnsight= has a *new* =-noPatches= option to suppress generation
of patches.
+ =foamToEnsightParts= has *new* =-noMesh= and =-index= options that can
be useful when post-processing results incrementally.
+ =snappyHexMesh= has lower memory footprint. New distributed triangulated
surface type for meshing surfaces with extremely large triangle count.
Now supports multi-region meshing of arbitrarily complex regions.
*** *New* utilities
+ =particleTracks= - generate particle tracks for lagrangian calculations.
+ =dsmcInitialise= - preprocessing utility to create initial configurations
of DSMC particles in a geometry.
+ =surfaceRedistributePar= - preprocessing utility to create distributed
triangulated surface.
*** *New* foamCalc functions
+ =interpolate= performs fvc::interpolate(<field>).
+ =randomise= randomises a <field> by a given perturbation.
+ =addSubtract= provides simple add/subtract field functionality.
*** Usage
+ =timeSelector= can now combine =-time ranges= and =-latestTime= options.
For example, -time '0.01:0.09' -latestTime vs. -time '0.01:'.
More reliable behaviour for cases missing /constant// or /0// directories.
When the =-noZero= option is enabled, =-latestTime= will not select the
=0/= directory unless the =-zeroTime= option is given.
This helps avoid ill effects caused by accidentally using the
/0// directory in certain utilities (eg, =reconstructPar=).
+ =-region= option added to more utilities.
*** Improvements to Paraview reader module
+ =PV3FoamReader= added mesh region handling. The region name is parsed
from the filename. Eg, /case{region}.OpenFOAM/.
+ =paraFoam= with a *new* =-region= option for specifying an alternative
region. A *new* =-touch= option to generate the /.OpenFOAM/ file only.
Only creates (and removes) /.OpenFOAM/ files if they didn't already
exist, which is useful in connection with the =-touch= option.
* Post-processing
+ Sampling on iso-surfaces, interpolated or non-interpolated.
+ Sampling on surface defined by distance to surface (=distanceSurface=).
+ Cutting planes for arbitrary meshes.
+ Output to any surface geometry format supported by the =surfMesh= library.
*** Function objects
***** Improvements for function objects and time-looping
+ The =functionObjectList= retains the order of the =functionObject=
order, which allows a chaining of operations. It is thus internally more
efficient when /system/controlDict/ uses =functions {..}= instead of
=functions (..)=, but both forms are supported.
+ The =functionObject= now has an additional =end()= method that is called
when =Time::loop()= or =Time::run()= determine that the time-loop exits.
Accordingly, one of these two idioms should be used in solver code:
1. =while (runTime.loop() { ... }=,
2. =while (runTime.run()) { runTime++; ... }=.
+ *New* =functionObjectList= now tracks the SHA1 message digest of the
sub-directories. This avoids reloading a =functionObject= when
something unrelated in /system/controlDict/ changed.
***** *New* function objects:
+ =systemCall= - executes a list of system instructions.
+ =fieldMinMax= - computes the min/max of a <field>.
+ =staticPressure= - converts kinematic pressure to static pressure.
+ =dsmcFields= - calculates intensive fields (velocity and temperature)
from averaged extensive fields (i.e. momentum and energy).
***** Usage
+ Improved output control: =timeStep= or =outputTime=.
* Tutorial restructuring
to reflect solver application structure.
* Third-party Software
+ =gcc= upgraded to version 4.3.3.
+ =OpenMPI= upgraded to version 1.3.3.
+ =ParaView= upgraded to version 3.6.1.
+ =Scotch= *new* decomposition method: \\
Scotch (http://gforge.inria.fr/projects/scotch/) is a general multi-level
decomposition method originating from the ScAlApplix project (Inria). It is
a framework for general recursive partitioning methods and a such comparable
to Metis but with a permissive licence.
The corresponding decomposition method (in =decomposeParDict=) is
=scotch=. An optional =strategy= string can be supplied to change the
decomposition methods; initial testing shows the default strategy producing
decompositions comparable in quality to Metis.

View File

@ -0,0 +1,3 @@
buoyantSimpleFoam.C
EXE = $(FOAM_APPBIN)/buoyantSimpleFoam

View File

@ -0,0 +1,13 @@
EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lmeshTools \
-lbasicThermophysicalModels \
-lspecie \
-lcompressibleRASModels \
-lfiniteVolume

View File

@ -0,0 +1,24 @@
// Solve the Momentum equation
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
+ turbulence->divDevRhoReff(U)
);
UEqn().relax();
solve
(
UEqn()
==
rho*g
- fvc::grad(p)
/*
fvc::reconstruct
(
fvc::interpolate(rho)*(g & mesh.Sf())
- fvc::snGrad(p)*mesh.magSf()
)
*/
);

View File

@ -0,0 +1,83 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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/>.
Application
buoyantSimpleFoam
Description
Steady-state solver for buoyant, turbulent flow of compressible fluids
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "basicPsiThermo.H"
#include "RASModel.H"
#include "fixedGradientFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "createFields.H"
#include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "readSIMPLEControls.H"
p.storePrevIter();
rho.storePrevIter();
// Pressure-velocity SIMPLE corrector
{
#include "UEqn.H"
#include "hEqn.H"
#include "pEqn.H"
}
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,69 @@
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<basicPsiThermo> pThermo
(
basicPsiThermo::New(mesh)
);
basicPsiThermo& thermo = pThermo();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
thermo.rho()
);
volScalarField& p = thermo.p();
volScalarField& h = thermo.h();
const volScalarField& psi = thermo.psi();
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "compressibleCreatePhi.H"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::RASModel> turbulence
(
compressible::RASModel::New
(
rho,
U,
phi,
thermo
)
);
thermo.correct();
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
mesh.solutionDict().subDict("SIMPLE"),
pRefCell,
pRefValue
);
dimensionedScalar initialMass = fvc::domainIntegrate(rho);

View File

@ -0,0 +1,17 @@
{
fvScalarMatrix hEqn
(
fvm::div(phi, h)
- fvm::Sp(fvc::div(phi), h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
- p*fvc::div(phi/fvc::interpolate(rho))
);
hEqn.relax();
hEqn.solve();
thermo.correct();
}

View File

@ -0,0 +1,57 @@
{
rho = thermo.rho();
volScalarField rUA = 1.0/UEqn().A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
U = rUA*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p);
surfaceScalarField buoyancyPhi =
rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
phi += buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rhorUAf, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (nonOrth == nNonOrthCorr)
{
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
}
// Calculate the conservative fluxes
phi -= pEqn.flux();
// Explicitly relax pressure for momentum corrector
p.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U += rUA*(rho*g - fvc::grad(p));
//U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rhorUAf);
U.correctBoundaryConditions();
}
}
#include "continuityErrs.H"
rho = thermo.rho();
rho.relax();
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
<< endl;
}

View File

@ -7,8 +7,8 @@ EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/barotropicCompressibilityModel/lnInclude
EXE_LIBS = \
-lincompressibleTurbulenceModel \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume \

View File

@ -20,3 +20,5 @@
{
solve(UEqn == -fvc::grad(p));
}
Info<< "max(U) " << max(mag(U)).value() << endl;

View File

@ -25,7 +25,8 @@ Application
cavitatingFoam
Description
Transient cavitation code based on the barotropic equation of state.
Transient cavitation code based on the homogeneous equilibrium model
from which the compressibility of the liquid/vapour "mixture" is obtained.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.

View File

@ -52,9 +52,32 @@
}
}
Info<< "max-min p: " << max(p).value()
Info<< "Predicted p max-min : " << max(p).value()
<< " " << min(p).value() << endl;
rho == max
(
psi*p
+ (1.0 - gamma)*rhol0
+ ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
rhoMin
);
#include "gammaPsi.H"
p =
(
rho
- (1.0 - gamma)*rhol0
- ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
)/psi;
p.correctBoundaryConditions();
Info<< "Phase-change corrected p max-min : " << max(p).value()
<< " " << min(p).value() << endl;
// Correct velocity
U = HbyA - rUA*fvc::grad(p);
@ -70,17 +93,4 @@
U.correctBoundaryConditions();
Info<< "max(U) " << max(mag(U)).value() << endl;
rho == max
(
psi*p
+ (1.0 - gamma)*rhol0
+ ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
rhoMin
);
Info<< "max-min rho: " << max(rho).value()
<< " " << min(rho).value() << endl;
#include "gammaPsi.H"
}

View File

@ -1,86 +0,0 @@
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phir = phic*interface.nHatf();
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
{
volScalarField::DimensionedInternalField Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
);
volScalarField::DimensionedInternalField Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
// Divergence term is handled explicitly to be
// consistent with the explicit transport solution
divU*min(alpha1, scalar(1))
);
forAll(dgdt, celli)
{
if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
{
Sp[celli] -= dgdt[celli]*alpha1[celli];
Su[celli] += dgdt[celli]*alpha1[celli];
}
else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
{
Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
}
}
surfaceScalarField phiAlpha1 =
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
);
MULES::explicitSolve
(
geometricOneField(),
alpha1,
phi,
phiAlpha1,
Sp,
Su,
1,
0
);
surfaceScalarField rho1f = fvc::interpolate(rho1);
surfaceScalarField rho2f = fvc::interpolate(rho2);
rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f;
alpha2 = scalar(1) - alpha1;
}
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(alpha1).value()
<< " Min(alpha2) = " << min(alpha2).value()
<< endl;
}

View File

@ -1,144 +0,0 @@
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field alpha1\n" << endl;
volScalarField alpha1
(
IOobject
(
"alpha1",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Calculating field alpha1\n" << endl;
volScalarField alpha2("alpha2", scalar(1) - alpha1);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
Info<< "Reading transportProperties\n" << endl;
twoPhaseMixture twoPhaseProperties(U, phi);
dimensionedScalar rho10
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase1Name()
).lookup("rho0")
);
dimensionedScalar rho20
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase2Name()
).lookup("rho0")
);
dimensionedScalar psi1
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase1Name()
).lookup("psi")
);
dimensionedScalar psi2
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase2Name()
).lookup("psi")
);
dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
volScalarField rho1 = rho10 + psi1*p;
volScalarField rho2 = rho20 + psi2*p;
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
alpha1*rho1 + alpha2*rho2
);
// Mass flux
// Initialisation does not matter because rhoPhi is reset after the
// alpha1 solution before it is used in the U equation.
surfaceScalarField rhoPhi
(
IOobject
(
"rho*phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
fvc::interpolate(rho)*phi
);
volScalarField dgdt =
pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001));
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
);
wordList pcorrTypes
(
p.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
forAll(p.boundaryField(), i)
{
if (p.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}

View File

@ -1,95 +0,0 @@
{
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA);
tmp<fvScalarMatrix> pEqnComp;
if (transonic)
{
pEqnComp =
(fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
}
else
{
pEqnComp =
(fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
}
U = rUA*UEqn.H();
surfaceScalarField phiU
(
"phiU",
(fvc::interpolate(U) & mesh.Sf())
);
phi = phiU +
(
fvc::interpolate(interface.sigmaK())
*fvc::snGrad(alpha1)*mesh.magSf()
+ fvc::interpolate(rho)*(g & mesh.Sf())
)*rUAf;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqnIncomp
(
fvc::div(phi)
- fvm::laplacian(rUAf, p)
);
if
(
oCorr == nOuterCorr-1
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
{
solve
(
(
max(alpha1, scalar(0))*(psi1/rho1)
+ max(alpha2, scalar(0))*(psi2/rho2)
)
*pEqnComp()
+ pEqnIncomp,
mesh.solver(p.name() + "Final")
);
}
else
{
solve
(
(
max(alpha1, scalar(0))*(psi1/rho1)
+ max(alpha2, scalar(0))*(psi2/rho2)
)
*pEqnComp()
+ pEqnIncomp
);
}
if (nonOrth == nNonOrthCorr)
{
dgdt =
(pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
*(pEqnComp & p);
phi += pEqnIncomp.flux();
}
}
U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
U.correctBoundaryConditions();
p.max(pMin);
rho1 = rho10 + psi1*p;
rho2 = rho20 + psi2*p;
Info<< "max(U) " << max(mag(U)).value() << endl;
Info<< "min(p) " << min(p).value() << endl;
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
}

View File

@ -0,0 +1,8 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wclean
wclean compressibleInterDyMFoam
# ----------------------------------------------------------------- end-of-file

View File

@ -0,0 +1,8 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wmake
wmake compressibleInterDyMFoam
# ----------------------------------------------------------------- end-of-file

View File

@ -6,7 +6,7 @@ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-linterfaceProperties \
-ltwoPhaseInterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \

View File

@ -24,10 +24,10 @@
==
fvc::reconstruct
(
fvc::interpolate(rho)*(g & mesh.Sf())
+ (
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- fvc::snGrad(p)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);

View File

@ -1,5 +1,5 @@
EXE_INC = \
-I../interFoam \
-I.. \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
@ -10,7 +10,7 @@ EXE_INC = \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude
EXE_LIBS = \
-linterfaceProperties \
-ltwoPhaseInterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
@ -18,5 +18,4 @@ EXE_LIBS = \
-lfiniteVolume \
-ldynamicMesh \
-lmeshTools \
-ldynamicFvMesh \
-ltopoChangerFvMesh
-ldynamicFvMesh

View File

@ -27,7 +27,7 @@
!(++alphaSubCycle).end();
)
{
# include "alphaEqns.H"
#include "alphaEqns.H"
rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
@ -35,7 +35,7 @@
}
else
{
# include "alphaEqns.H"
#include "alphaEqns.H"
}
if (oCorr == 0)

View File

@ -22,7 +22,7 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
compressibleLesInterFoam
compressibleInterDyMFoam
Description
Solver for 2 compressible, isothermal immiscible fluids using a VOF
@ -30,9 +30,10 @@ Description
with optional mesh motion and mesh topology changes including adaptive
re-meshing.
The momentum and other fluid properties are of the "mixture" and a
single momentum equation is solved. Turbulence modelling is generic,
i.e. laminar, RAS or LES may be selected.
The momentum and other fluid properties are of the "mixture" and a single
momentum equation is solved.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
\*---------------------------------------------------------------------------*/
@ -55,6 +56,7 @@ int main(int argc, char *argv[])
#include "readControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createPcorrTypes.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
@ -89,6 +91,9 @@ int main(int argc, char *argv[])
Info<< "Execution time for mesh.update() = "
<< runTime.elapsedCpuTime() - timeBeforeMeshUpdate
<< " s" << endl;
gh = g & mesh.C();
ghf = g & mesh.Cf();
}
if (mesh.changing() && correctPhi)

View File

@ -42,7 +42,7 @@
adjustPhi(phi, U, pcorr);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pcorrEqn
(

View File

@ -0,0 +1,13 @@
wordList pcorrTypes
(
p_rgh.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
for (label i=0; i<p_rgh.boundaryField().size(); i++)
{
if (p_rgh.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}

View File

@ -0,0 +1,96 @@
{
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA);
tmp<fvScalarMatrix> p_rghEqnComp;
if (transonic)
{
p_rghEqnComp =
(
fvm::ddt(p_rgh)
+ fvm::div(phi, p_rgh)
- fvm::Sp(fvc::div(phi), p_rgh)
);
}
else
{
p_rghEqnComp =
(
fvm::ddt(p_rgh)
+ fvc::div(phi, p_rgh)
- fvc::Sp(fvc::div(phi), p_rgh)
);
}
U = rUA*UEqn.H();
surfaceScalarField phiU
(
"phiU",
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
);
phi = phiU +
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf();
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix p_rghEqnIncomp
(
fvc::div(phi)
- fvm::laplacian(rUAf, p_rgh)
);
solve
(
(
max(alpha1, scalar(0))*(psi1/rho1)
+ max(alpha2, scalar(0))*(psi2/rho2)
)
*p_rghEqnComp()
+ p_rghEqnIncomp,
mesh.solver
(
p_rgh.select
(
oCorr == nOuterCorr-1
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
);
if (nonOrth == nNonOrthCorr)
{
dgdt =
(pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
*(p_rghEqnComp & p_rgh);
phi += p_rghEqnIncomp.flux();
}
}
U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
U.correctBoundaryConditions();
p = max
(
(p_rgh + gh*(alpha1*rho10 + alpha2*rho20))
/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
pMin
);
rho1 = rho10 + psi1*p;
rho2 = rho20 + psi2*p;
Info<< "max(U) " << max(mag(U)).value() << endl;
Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
}

View File

@ -1,5 +1,5 @@
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "readPISOControls.H"
#include "readTimeControls.H"
label nAlphaCorr
(
@ -19,9 +19,14 @@
<< exit(FatalError);
}
const bool correctPhi =
piso.lookupOrDefault("correctPhi", true);
const bool checkMeshCourantNo =
piso.lookupOrDefault("checkMeshCourantNo", false);
bool correctPhi = true;
if (piso.found("correctPhi"))
{
correctPhi = Switch(piso.lookup("correctPhi"));
}
bool checkMeshCourantNo = false;
if (piso.found("checkMeshCourantNo"))
{
checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
}

View File

@ -1,9 +1,9 @@
Info<< "Reading field p\n" << endl;
volScalarField p
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p",
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
@ -83,6 +83,28 @@
dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
max
(
(p_rgh + gh*(alpha1*rho10 + alpha2*rho20))
/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
pMin
)
);
volScalarField rho1 = rho10 + psi1*p;
volScalarField rho2 = rho20 + psi2*p;

View File

@ -2,17 +2,25 @@
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA);
tmp<fvScalarMatrix> pEqnComp;
tmp<fvScalarMatrix> p_rghEqnComp;
if (transonic)
{
pEqnComp =
(fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
p_rghEqnComp =
(
fvm::ddt(p_rgh)
+ fvm::div(phi, p_rgh)
- fvm::Sp(fvc::div(phi), p_rgh)
);
}
else
{
pEqnComp =
(fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
p_rghEqnComp =
(
fvm::ddt(p_rgh)
+ fvc::div(phi, p_rgh)
- fvc::Sp(fvc::div(phi), p_rgh)
);
}
@ -21,22 +29,22 @@
surfaceScalarField phiU
(
"phiU",
(fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi)
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
);
phi = phiU +
(
fvc::interpolate(interface.sigmaK())
*fvc::snGrad(alpha1)*mesh.magSf()
+ fvc::interpolate(rho)*(g & mesh.Sf())
)*rUAf;
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqnIncomp
fvScalarMatrix p_rghEqnIncomp
(
fvc::div(phi)
- fvm::laplacian(rUAf, p)
- fvm::laplacian(rUAf, p_rgh)
);
solve
@ -45,27 +53,41 @@
max(alpha1, scalar(0))*(psi1/rho1)
+ max(alpha2, scalar(0))*(psi2/rho2)
)
*pEqnComp()
+ pEqnIncomp
*p_rghEqnComp()
+ p_rghEqnIncomp,
mesh.solver
(
p_rgh.select
(
oCorr == nOuterCorr-1
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
);
if (nonOrth == nNonOrthCorr)
{
dgdt =
(pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
*(pEqnComp & p);
phi += pEqnIncomp.flux();
*(p_rghEqnComp & p_rgh);
phi += p_rghEqnIncomp.flux();
}
}
U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
U.correctBoundaryConditions();
p.max(pMin);
p = max
(
(p_rgh + gh*(alpha1*rho10 + alpha2*rho20))
/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
pMin
);
rho1 = rho10 + psi1*p;
rho2 = rho20 + psi2*p;
Info<< "max(U) " << max(mag(U)).value() << endl;
Info<< "min(p) " << min(p).value() << endl;
Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
}

View File

@ -1,62 +0,0 @@
{
if (mesh.changing())
{
forAll(U.boundaryField(), patchi)
{
if (U.boundaryField()[patchi].fixesValue())
{
U.boundaryField()[patchi].initEvaluate();
}
}
forAll(U.boundaryField(), patchi)
{
if (U.boundaryField()[patchi].fixesValue())
{
U.boundaryField()[patchi].evaluate();
phi.boundaryField()[patchi] =
U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi];
}
}
}
#include "continuityErrs.H"
volScalarField pcorr
(
IOobject
(
"pcorr",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("pcorr", p.dimensions(), 0.0),
pcorrTypes
);
dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
adjustPhi(phi, U, pcorr);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
);
pcorrEqn.setReference(pRefCell, pRefValue);
pcorrEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phi -= pcorrEqn.flux();
}
}
#include "continuityErrs.H"
}

View File

@ -1,112 +0,0 @@
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field alpha1\n" << endl;
volScalarField alpha1
(
IOobject
(
"alpha1",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
Info<< "Reading transportProperties\n" << endl;
twoPhaseMixture twoPhaseProperties(U, phi);
const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + (scalar(1) - alpha1)*rho2,
alpha1.boundaryField().types()
);
rho.oldTime();
// Mass flux
// Initialisation does not matter because rhoPhi is reset after the
// alpha1 solution before it is used in the U equation.
surfaceScalarField rhoPhi
(
IOobject
(
"rho*phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
rho1*phi
);
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
);
wordList pcorrTypes
(
p.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
forAll(p.boundaryField(), i)
{
if (p.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);

View File

@ -1,52 +0,0 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
U = rAU*UEqn.H();
surfaceScalarField phiU("phiU", (fvc::interpolate(U) & mesh.Sf()));
if (p.needReference())
{
fvc::makeRelative(phiU, U);
adjustPhi(phiU, U, p);
fvc::makeAbsolute(phiU, U);
}
phi = phiU +
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+ fvc::interpolate(rho)*(g & mesh.Sf())
)*rAUf;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rAUf, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pEqn.solve(mesh.solver(p.name()));
}
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
U.correctBoundaryConditions();
#include "continuityErrs.H"
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
}

View File

@ -1,9 +0,0 @@
# include "readTimeControls.H"
# include "readPISOControls.H"
const bool correctPhi =
piso.lookupOrDefault("correctPhi", true);
const bool checkMeshCourantNo =
piso.lookupOrDefault("checkMeshCourantNo", false);

View File

@ -0,0 +1,10 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wclean
wclean interDyMFoam
wclean MRFInterFoam
wclean porousInterFoam
# ----------------------------------------------------------------- end-of-file

View File

@ -0,0 +1,10 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wmake
wmake interDyMFoam
wmake MRFInterFoam
wmake porousInterFoam
# ----------------------------------------------------------------- end-of-file

View File

@ -51,7 +51,6 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
@ -70,6 +69,7 @@ int main(int argc, char *argv[])
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
runTime++;
@ -79,6 +79,7 @@ int main(int argc, char *argv[])
twoPhaseProperties.correct();
#include "alphaEqnSubCycle.H"
#include "zonePhaseVolumes.H"
#include "UEqn.H"
@ -88,8 +89,6 @@ int main(int argc, char *argv[])
#include "pEqn.H"
}
#include "continuityErrs.H"
turbulence->correct();
runTime.write();

View File

@ -0,0 +1,3 @@
MRFInterFoam.C
EXE = $(FOAM_APPBIN)/MRFInterFoam

View File

@ -1,5 +1,5 @@
EXE_INC = \
-I$(FOAM_SOLVERS)/multiphase/interFoam \
-I.. \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
@ -7,8 +7,9 @@ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-linterfaceProperties \
-ltwoPhaseInterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume

View File

@ -25,10 +25,10 @@
==
fvc::reconstruct
(
fvc::interpolate(rho)*(g & mesh.Sf())
+ (
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- fvc::snGrad(p)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);

View File

@ -0,0 +1,62 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
U = rAU*UEqn.H();
surfaceScalarField phiU
(
"phiU",
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phi)
);
mrfZones.relativeFlux(phiU);
adjustPhi(phiU, U, p_rgh);
phi = phiU +
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf();
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
p_rghEqn.solve
(
mesh.solver
(
p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
)
);
if (nonOrth == nNonOrthCorr)
{
phi -= p_rghEqn.flux();
}
}
U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
U.correctBoundaryConditions();
#include "continuityErrs.H"
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
}

View File

@ -0,0 +1,21 @@
{
const scalarField& V = mesh.V();
forAll(mesh.cellZones(), czi)
{
const labelList& cellLabels = mesh.cellZones()[czi];
scalar phaseVolume = 0;
forAll(cellLabels, cli)
{
label celli = cellLabels[cli];
phaseVolume += alpha1[celli]*V[celli];
}
reduce(phaseVolume, sumOp<scalar>());
Info<< "Phase volume in zone " << mesh.cellZones()[czi].name()
<< " = " << phaseVolume*1e6 << " ml " << endl;
}
}

View File

@ -6,7 +6,7 @@ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-linterfaceProperties \
-ltwoPhaseInterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \

View File

@ -1,6 +0,0 @@
surfaceScalarField alpha1f = fvc::interpolate(alpha1);
surfaceScalarField UBlendingFactor
(
"UBlendingFactor",
sqrt(max(min(4*alpha1f*(1.0 - alpha1f), 1.0), 0.0))
);

View File

@ -24,10 +24,10 @@
==
fvc::reconstruct
(
fvc::interpolate(rho)*(g & mesh.Sf())
+ (
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- fvc::snGrad(p)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);

View File

@ -0,0 +1,58 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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/>.
Global
CourantNo
Description
Calculates and outputs the mean and maximum Courant Numbers.
\*---------------------------------------------------------------------------*/
scalar maxAlphaCo
(
readScalar(runTime.controlDict().lookup("maxAlphaCo"))
);
scalar alphaCoNum = 0.0;
scalar meanAlphaCoNum = 0.0;
if (mesh.nInternalFaces())
{
surfaceScalarField alphaf = fvc::interpolate(alpha1);
surfaceScalarField SfUfbyDelta =
pos(alphaf - 0.01)*pos(0.99 - alphaf)
*mesh.surfaceInterpolation::deltaCoeffs()*mag(phi);
alphaCoNum = max(SfUfbyDelta/mesh.magSf())
.value()*runTime.deltaT().value();
meanAlphaCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
.value()*runTime.deltaT().value();
}
Info<< "Interface Courant Number mean: " << meanAlphaCoNum
<< " max: " << alphaCoNum << endl;
// ************************************************************************* //

View File

@ -3,13 +3,13 @@
wordList pcorrTypes
(
p.boundaryField().size(),
p_rgh.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
forAll(p.boundaryField(), i)
forAll (p_rgh.boundaryField(), i)
{
if (p.boundaryField()[i].fixesValue())
if (p_rgh.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
@ -26,11 +26,11 @@
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("pcorr", p.dimensions(), 0.0),
dimensionedScalar("pcorr", p_rgh.dimensions(), 0.0),
pcorrTypes
);
dimensionedScalar rUAf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
adjustPhi(phi, U, pcorr);
@ -38,7 +38,7 @@
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
);
pcorrEqn.setReference(pRefCell, pRefValue);

View File

@ -1,9 +1,9 @@
Info<< "Reading field p\n" << endl;
volScalarField p
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p",
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
@ -40,7 +40,7 @@
mesh
);
# include "createPhi.H"
#include "createPhi.H"
Info<< "Reading transportProperties\n" << endl;
@ -83,11 +83,6 @@
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties);
@ -97,3 +92,54 @@
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
);
#include "readGravitationalAcceleration.H"
/*
dimensionedVector g0(g);
// Read the data file and initialise the interpolation table
interpolationTable<vector> timeSeriesAcceleration
(
runTime.path()/runTime.caseConstant()/"acceleration.dat"
);
*/
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
mesh.solutionDict().subDict("PISO"),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}

View File

@ -1,4 +1,5 @@
EXE_INC = \
-I.. \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
@ -9,7 +10,7 @@ EXE_INC = \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude
EXE_LIBS = \
-linterfaceProperties \
-ltwoPhaseInterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
@ -19,4 +20,3 @@ EXE_LIBS = \
-lmeshTools \
-ldynamicFvMesh \
-ltopoChangerFvMesh

View File

@ -47,7 +47,6 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "readGravitationalAcceleration.H"
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
@ -62,6 +61,7 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readControls.H"
#include "alphaCourantNo.H"
#include "CourantNo.H"
// Make the fluxes absolute
@ -83,6 +83,9 @@ int main(int argc, char *argv[])
Info<< "Execution time for mesh.update() = "
<< runTime.elapsedCpuTime() - timeBeforeMeshUpdate
<< " s" << endl;
gh = g & mesh.C();
ghf = g & mesh.Cf();
}
if (mesh.changing() && correctPhi)

View File

@ -0,0 +1,64 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
U = rAU*UEqn.H();
surfaceScalarField phiU("phiU", (fvc::interpolate(U) & mesh.Sf()));
if (p_rgh.needReference())
{
fvc::makeRelative(phiU, U);
adjustPhi(phiU, U, p_rgh);
fvc::makeAbsolute(phiU, U);
}
phi = phiU +
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf();
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
p_rghEqn.solve
(
mesh.solver
(
p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
)
);
if (nonOrth == nNonOrthCorr)
{
phi -= p_rghEqn.flux();
}
}
U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
U.correctBoundaryConditions();
#include "continuityErrs.H"
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
}

View File

@ -0,0 +1,14 @@
# include "readTimeControls.H"
# include "readPISOControls.H"
bool correctPhi = true;
if (piso.found("correctPhi"))
{
correctPhi = Switch(piso.lookup("correctPhi"));
}
bool checkMeshCourantNo = false;
if (piso.found("checkMeshCourantNo"))
{
checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
}

View File

@ -43,6 +43,7 @@ Description
#include "interfaceProperties.H"
#include "twoPhaseMixture.H"
#include "turbulenceModel.H"
#include "interpolationTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -51,7 +52,6 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
@ -69,6 +69,7 @@ int main(int argc, char *argv[])
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
runTime++;
@ -87,8 +88,6 @@ int main(int argc, char *argv[])
#include "pEqn.H"
}
#include "continuityErrs.H"
turbulence->correct();
runTime.write();

View File

@ -1,49 +1,61 @@
{
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA);
U = rUA*UEqn.H();
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
U = rAU*UEqn.H();
surfaceScalarField phiU
(
"phiU",
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
+ fvc::ddtPhiCorr(rAU, rho, U, phi)
);
adjustPhi(phiU, U, p);
adjustPhi(phiU, U, p_rgh);
phi = phiU +
(
fvc::interpolate(interface.sigmaK())
*fvc::snGrad(alpha1)*mesh.magSf()
+ fvc::interpolate(rho)*(g & mesh.Sf())
)*rUAf;
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rUAf, p) == fvc::div(phi)
fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pEqn.solve(mesh.solver(p.name()));
}
p_rghEqn.solve
(
mesh.solver
(
p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
)
);
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
phi -= p_rghEqn.flux();
}
}
U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
U.correctBoundaryConditions();
#include "continuityErrs.H"
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
}

View File

@ -0,0 +1,3 @@
porousInterFoam.C
EXE = $(FOAM_APPBIN)/porousInterFoam

View File

@ -0,0 +1,17 @@
EXE_INC = \
-I.. \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-ltwoPhaseInterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume \
-lmeshTools

View File

@ -5,17 +5,22 @@
+ fvc::interpolate(rho*turbulence->nut())
);
// Calculate and cache mu for the porous media
volScalarField mu(twoPhaseProperties.mu());
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
pZones.ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muEff, U)
- (fvc::grad(U) & fvc::grad(muEff))
//- fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T())))
//- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
);
UEqn.relax();
pZones.addResistance(UEqn);
if (momentumPredictor)
{
solve
@ -24,10 +29,10 @@
==
fvc::reconstruct
(
fvc::interpolate(rho)*(g & mesh.Sf())
+ (
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- fvc::snGrad(p)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);

View File

@ -0,0 +1 @@
porousZones pZones(mesh);

View File

@ -0,0 +1,108 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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/>.
Application
porousInterFoam
Description
Solver for 2 incompressible, isothermal immiscible fluids using a VOF
(volume of fluid) phase-fraction based interface capturing approach.
The momentum and other fluid properties are of the "mixture" and a single
momentum equation is solved.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
For a two-fluid approach see twoPhaseEulerFoam.
Explicit handling of porous zones is included.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "MULES.H"
#include "subCycle.H"
#include "interfaceProperties.H"
#include "twoPhaseMixture.H"
#include "turbulenceModel.H"
#include "porousZones.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createPorousZones.H"
#include "readTimeControls.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
twoPhaseProperties.correct();
#include "alphaEqnSubCycle.H"
#include "UEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,53 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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/>.
Global
setDeltaT
Description
Reset the timestep to maintain a constant maximum courant Number.
Reduction of time-step is immediate, but increase is damped to avoid
unstable oscillations.
\*---------------------------------------------------------------------------*/
if (adjustTimeStep)
{
scalar maxDeltaTFact =
min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
runTime.setDeltaT
(
min
(
deltaTFact*runTime.deltaT().value(),
maxDeltaT
)
);
Info<< "deltaT = " << runTime.deltaT().value() << endl;
}
// ************************************************************************* //

View File

@ -5,12 +5,13 @@ EXE_INC = \
-IincompressibleThreePhaseMixture \
-IthreePhaseInterfaceProperties \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/transportModels/twoPhaseInterfaceProperties/alphaContactAngle/alphaContactAngle \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/transportModels
EXE_LIBS = \
-linterfaceProperties \
-ltwoPhaseInterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \

View File

@ -0,0 +1,61 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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/>.
Global
CourantNo
Description
Calculates and outputs the mean and maximum Courant Numbers.
\*---------------------------------------------------------------------------*/
scalar maxAlphaCo
(
readScalar(runTime.controlDict().lookup("maxAlphaCo"))
);
scalar alphaCoNum = 0.0;
scalar meanAlphaCoNum = 0.0;
if (mesh.nInternalFaces())
{
surfaceScalarField alpha1f = fvc::interpolate(alpha1);
surfaceScalarField alpha2f = fvc::interpolate(alpha2);
surfaceScalarField SfUfbyDelta = max
(
pos(alpha1f - 0.01)*pos(0.99 - alpha1f),
pos(alpha2f - 0.01)*pos(0.99 - alpha2f)
)*mesh.surfaceInterpolation::deltaCoeffs()*mag(phi);
alphaCoNum = max(SfUfbyDelta/mesh.magSf())
.value()*runTime.deltaT().value();
meanAlphaCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
.value()*runTime.deltaT().value();
}
Info<< "Interface Courant Number mean: " << meanAlphaCoNum
<< " max: " << alphaCoNum << endl;
// ************************************************************************* //

View File

@ -1,9 +1,9 @@
Info<< "Reading field p\n" << endl;
volScalarField p
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p",
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
@ -73,7 +73,7 @@
mesh
);
# include "createPhi.H"
#include "createPhi.H"
threePhaseMixture threePhaseProperties(U, phi);
@ -116,11 +116,6 @@
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
// Construct interface from alpha distribution
threePhaseInterfaceProperties interface(threePhaseProperties);
@ -130,3 +125,43 @@
(
incompressible::turbulenceModel::New(U, phi, threePhaseProperties)
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
mesh.solutionDict().subDict("PISO"),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}

View File

@ -62,6 +62,7 @@ int main(int argc, char *argv[])
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
runTime++;

View File

@ -1,6 +1,6 @@
interPhaseChangeFoam.C
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.C
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixtureNew.C
phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/newPhaseChangeTwoPhaseMixture.C
phaseChangeTwoPhaseMixtures/Kunz/Kunz.C
phaseChangeTwoPhaseMixtures/Merkle/Merkle.C
phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C

View File

@ -7,7 +7,7 @@ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-linterfaceProperties \
-ltwoPhaseInterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \

View File

@ -25,10 +25,10 @@
==
fvc::reconstruct
(
fvc::interpolate(rho)*(g & mesh.Sf())
+ (
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- fvc::snGrad(p)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);

View File

@ -50,18 +50,18 @@
+ vDotcAlphal
);
// MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
// MULES::explicitSolve
// (
// geometricOneField(),
// alpha1,
// phi,
// phiAlpha,
// Sp,
// Su,
// 1,
// 0
// );
//MULES::explicitSolve
//(
// geometricOneField(),
// alpha1,
// phi,
// phiAlpha,
// Sp,
// Su,
// 1,
// 0
//);
MULES::implicitSolve
(
geometricOneField(),

View File

@ -36,12 +36,12 @@ surfaceScalarField rhoPhi
!(++alphaSubCycle).end();
)
{
# include "alphaEqn.H"
#include "alphaEqn.H"
}
}
else
{
# include "alphaEqn.H"
#include "alphaEqn.H"
}
if (nOuterCorr == 1)

View File

@ -1,54 +0,0 @@
{
#include "continuityErrs.H"
wordList pcorrTypes
(
p.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
forAll(p.boundaryField(), i)
{
if (p.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}
volScalarField pcorr
(
IOobject
(
"pcorr",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("pcorr", p.dimensions(), 0.0),
pcorrTypes
);
dimensionedScalar rUAf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
adjustPhi(phi, U, pcorr);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
);
pcorrEqn.setReference(pRefCell, pRefValue);
pcorrEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phi -= pcorrEqn.flux();
}
}
# include "continuityErrs.H"
}

View File

@ -1,9 +1,9 @@
Info<< "Reading field p\n" << endl;
volScalarField p
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p",
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
@ -40,7 +40,7 @@
mesh
);
# include "createPhi.H"
#include "createPhi.H"
Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
autoPtr<phaseChangeTwoPhaseMixture> twoPhaseProperties =
@ -65,12 +65,6 @@
);
rho.oldTime();
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties());
@ -79,3 +73,43 @@
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties())
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
mesh.solutionDict().subDict("PISO"),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}

View File

@ -59,7 +59,7 @@ int main(int argc, char *argv[])
#include "initContinuityErrs.H"
#include "createFields.H"
#include "readTimeControls.H"
#include "correctPhi.H"
#include "../interFoam/correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
@ -82,9 +82,11 @@ int main(int argc, char *argv[])
turbulence->correct();
// --- Outer-corrector loop
// --- Pressure-velocity PIMPLE corrector loop
for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
{
bool finalIter = oCorr == nOuterCorr-1;
#include "UEqn.H"
// --- PISO loop
@ -92,8 +94,6 @@ int main(int argc, char *argv[])
{
#include "pEqn.H"
}
#include "continuityErrs.H"
}
twoPhaseProperties->correct();

View File

@ -11,14 +11,13 @@
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
);
adjustPhi(phiU, U, p);
adjustPhi(phiU, U, p_rgh);
phi = phiU +
(
fvc::interpolate(interface.sigmaK())
*fvc::snGrad(alpha1)*mesh.magSf()
+ fvc::interpolate(rho)*(g & mesh.Sf())
)*rUAf;
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf();
Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
const volScalarField& vDotcP = vDotP[0]();
@ -26,29 +25,48 @@
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
fvScalarMatrix p_rghEqn
(
fvc::div(phi) - fvm::laplacian(rUAf, p)
- (vDotvP - vDotcP)*pSat + fvm::Sp(vDotvP - vDotcP, p)
fvc::div(phi) - fvm::laplacian(rUAf, p_rgh)
- (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
);
pEqn.setReference(pRefCell, pRefValue);
p_rghEqn.setReference(pRefCell, pRefValue);
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pEqn.solve(mesh.solver(p.name() + "Final"));
}
else
{
pEqn.solve(mesh.solver(p.name()));
}
p_rghEqn.solve
(
mesh.solver
(
p_rgh.select
(
finalIter
&& corr == nCorr-1
&& nonOrth == nNonOrthCorr
)
)
);
if (nonOrth == nNonOrthCorr)
{
phi += pEqn.flux();
phi += p_rghEqn.flux();
}
}
U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
U.correctBoundaryConditions();
#include "continuityErrs.H"
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
}

View File

@ -9,16 +9,16 @@ 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 Generac License as published by
the Free Software Foundation; either 2 of the License, or
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 ho it will be useful, but WITHOUT
ANY WARRANTY; without even the imarranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.he GNU General Public License
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 oNU General Public License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class

View File

@ -9,16 +9,16 @@ 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 Generac License as published by
the Free Software Foundation; either 2 of the License, or
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 ho it will be useful, but WITHOUT
ANY WARRANTY; without even the imarranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.he GNU General Public License
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 oNU General Public License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class

View File

@ -1,5 +1,5 @@
/*---------------------------------------------------------------------------*\
========Merkle= |
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
@ -9,16 +9,16 @@ 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 Generac License as published by
the Free Software Foundation; either 2 of the License, or
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 ho it will be useful, but WITHOUT
ANY WARRANTY; without even the imarranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.he GNU General Public License
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 oNU General Public License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class

View File

@ -36,27 +36,30 @@ Foam::phaseChangeTwoPhaseMixture::New
const word& alpha1Name
)
{
// get model name, but do not register the dictionary
const word mixtureType
IOdictionary transportPropertiesDict
(
IOdictionary
IOobject
(
IOobject
(
"transportProperties",
U.time().constant(),
U.db(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
).lookup("phaseChangeTwoPhaseMixture")
"transportProperties",
U.time().constant(),
U.db(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
);
Info<< "Selecting phaseChange model " << mixtureType << endl;
word phaseChangeTwoPhaseMixtureTypeName
(
transportPropertiesDict.lookup("phaseChangeTwoPhaseMixture")
);
Info<< "Selecting phaseChange model "
<< phaseChangeTwoPhaseMixtureTypeName << endl;
componentsConstructorTable::iterator cstrIter =
componentsConstructorTablePtr_->find(mixtureType);
componentsConstructorTablePtr_
->find(phaseChangeTwoPhaseMixtureTypeName);
if (cstrIter == componentsConstructorTablePtr_->end())
{
@ -64,8 +67,8 @@ Foam::phaseChangeTwoPhaseMixture::New
(
"phaseChangeTwoPhaseMixture::New"
) << "Unknown phaseChangeTwoPhaseMixture type "
<< mixtureType << nl << nl
<< "Valid phaseChangeTwoPhaseMixture types are : " << endl
<< phaseChangeTwoPhaseMixtureTypeName << endl << endl
<< "Valid phaseChangeTwoPhaseMixtures are : " << endl
<< componentsConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}

View File

@ -28,7 +28,7 @@ Description
SourceFiles
phaseChangeTwoPhaseMixture.C
phaseChangeModelNew.C
newPhaseChangeModel.C
\*---------------------------------------------------------------------------*/

View File

@ -0,0 +1,8 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wclean
wclean MRFMultiphaseInterFoam
# ----------------------------------------------------------------- end-of-file

View File

@ -0,0 +1,9 @@
#!/bin/sh
cd ${0%/*} || exit 1 # run from this directory
set -x
wmake libso multiphaseMixture
wmake
wmake MRFMultiphaseInterFoam
# ----------------------------------------------------------------- end-of-file

View File

@ -0,0 +1,99 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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/>.
Application
MRFMultiphaseInterFoam
Description
Solver for n incompressible fluids which captures the interfaces and
includes surface-tension and contact-angle effects for each phase.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "multiphaseMixture.H"
#include "turbulenceModel.H"
#include "MRFZones.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createMRFZones.H"
#include "readTimeControls.H"
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
mixture.solve();
rho = mixture.rho();
#include "zonePhaseVolumes.H"
#include "UEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,3 @@
MRFMultiphaseInterFoam.C
EXE = $(FOAM_APPBIN)/MRFMultiphaseInterFoam

View File

@ -0,0 +1,19 @@
EXE_INC = \
-I.. \
-I../../interFoam \
-I../../interFoam/MRFInterFoam \
-I../multiphaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lmultiphaseInterFoam \
-linterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume

View File

@ -0,0 +1,35 @@
surfaceScalarField muEff
(
"muEff",
mixture.muf()
+ fvc::interpolate(rho*turbulence->nut())
);
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(mixture.rhoPhi(), U)
- fvm::laplacian(muEff, U)
- (fvc::grad(U) & fvc::grad(muEff))
//- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
);
mrfZones.addCoriolis(rho, UEqn);
UEqn.relax();
if (momentumPredictor)
{
solve
(
UEqn
==
fvc::reconstruct
(
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);
}

View File

@ -0,0 +1,63 @@
{
volScalarField rAU = 1.0/UEqn.A();
surfaceScalarField rAUf = fvc::interpolate(rAU);
U = rAU*UEqn.H();
surfaceScalarField phiU
(
"phiU",
(fvc::interpolate(U) & mesh.Sf())
//+ fvc::ddtPhiCorr(rAU, rho, U, phi)
);
mrfZones.relativeFlux(phiU);
adjustPhi(phiU, U, p_rgh);
phi = phiU +
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
)*rAUf*mesh.magSf();
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
);
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
p_rghEqn.solve
(
mesh.solver
(
p_rgh.select(corr == nCorr-1 && nonOrth == nNonOrthCorr)
)
);
if (nonOrth == nNonOrthCorr)
{
phi -= p_rghEqn.flux();
}
}
U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
U.correctBoundaryConditions();
#include "continuityErrs.H"
p == p_rgh + rho*gh;
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rho*gh;
}
}

View File

@ -0,0 +1,26 @@
{
const scalarField& V = mesh.V();
forAll(mesh.cellZones(), czi)
{
const labelList& cellLabels = mesh.cellZones()[czi];
forAllConstIter(PtrDictionary<phase>, mixture.phases(), iter)
{
const volScalarField& alpha = iter();
scalar phaseVolume = 0;
forAll(cellLabels, cli)
{
label celli = cellLabels[cli];
phaseVolume += alpha[celli]*V[celli];
}
reduce(phaseVolume, sumOp<scalar>());
Info<< alpha.name()
<< " phase volume in zone " << mesh.cellZones()[czi].name()
<< " = " << phaseVolume*1e6 << " ml " << endl;
}
}
}

View File

@ -1,6 +1,3 @@
multiphaseMixture/phase/phase.C
multiphaseMixture/multiphaseAlphaContactAngle/multiphaseAlphaContactAngleFvPatchScalarField.C
multiphaseMixture/multiphaseMixture.C
multiphaseInterFoam.C
EXE = $(FOAM_APPBIN)/multiphaseInterFoam

View File

@ -1,17 +1,17 @@
EXE_INC = \
-I../interFoam \
-ImultiphaseMixture \
-ImultiphaseMixture/phase \
-ImultiphaseMixture/multiphaseAlphaContactAngle \
-ImultiphaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lmultiphaseInterFoam \
-linterfaceProperties \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume
-lfiniteVolume

View File

@ -24,10 +24,10 @@
==
fvc::reconstruct
(
fvc::interpolate(rho)*(g & mesh.Sf())
+ (
(
mixture.surfaceTensionForce()
- fvc::snGrad(p)
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);

View File

@ -0,0 +1,56 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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/>.
Global
CourantNo
Description
Calculates and outputs the mean and maximum Courant Numbers.
\*---------------------------------------------------------------------------*/
scalar maxAlphaCo
(
readScalar(runTime.controlDict().lookup("maxAlphaCo"))
);
scalar alphaCoNum = 0.0;
scalar meanAlphaCoNum = 0.0;
if (mesh.nInternalFaces())
{
surfaceScalarField SfUfbyDelta =
mixture.nearInterface()
*mesh.surfaceInterpolation::deltaCoeffs()*mag(phi);
alphaCoNum = max(SfUfbyDelta/mesh.magSf())
.value()*runTime.deltaT().value();
meanAlphaCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
.value()*runTime.deltaT().value();
}
Info<< "Interface Courant Number mean: " << meanAlphaCoNum
<< " max: " << alphaCoNum << endl;
// ************************************************************************* //

View File

@ -1,9 +1,9 @@
Info<< "Reading field p\n" << endl;
volScalarField p
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p",
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
@ -45,13 +45,48 @@
rho.oldTime();
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, mixture)
);
#include "readGravitationalAcceleration.H"
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
p_rgh + rho*gh
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
p,
p_rgh,
mesh.solutionDict().subDict("PISO"),
pRefCell,
pRefValue
);
if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
}

View File

@ -43,7 +43,6 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "readPISOControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
@ -61,13 +60,14 @@ int main(int argc, char *argv[])
#include "readPISOControls.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "alphaCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
mixture.correct();
mixture.solve();
rho = mixture.rho();
#include "UEqn.H"
@ -78,18 +78,16 @@ int main(int argc, char *argv[])
#include "pEqn.H"
}
#include "continuityErrs.H"
//turbulence->correct();
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n" << endl << endl;
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "\n end \n";
Info<< "End\n" << endl;
return 0;
}

View File

@ -0,0 +1,5 @@
phase/phase.C
alphaContactAngle/alphaContactAngleFvPatchScalarField.C
multiphaseMixture.C
LIB = $(FOAM_LIBBIN)/libmultiphaseInterFoam

View File

@ -0,0 +1,11 @@
EXE_INC = \
-IalphaContactAngle \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-linterfaceProperties \
-lincompressibleTransportModels \
-lfiniteVolume

View File

@ -23,7 +23,7 @@ License
\*---------------------------------------------------------------------------*/
#include "multiphaseAlphaContactAngleFvPatchScalarField.H"
#include "alphaContactAngleFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
@ -32,8 +32,7 @@ License
namespace Foam
{
multiphaseAlphaContactAngleFvPatchScalarField::interfaceThetaProps::
interfaceThetaProps
alphaContactAngleFvPatchScalarField::interfaceThetaProps::interfaceThetaProps
(
Istream& is
)
@ -48,10 +47,10 @@ interfaceThetaProps
Istream& operator>>
(
Istream& is,
multiphaseAlphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
)
{
is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_;
is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_;
return is;
}
@ -59,7 +58,7 @@ Istream& operator>>
Ostream& operator<<
(
Ostream& os,
const multiphaseAlphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
const alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
)
{
os << tp.theta0_ << token::SPACE
@ -73,8 +72,7 @@ Ostream& operator<<
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
multiphaseAlphaContactAngleFvPatchScalarField::
multiphaseAlphaContactAngleFvPatchScalarField
alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
@ -84,10 +82,9 @@ multiphaseAlphaContactAngleFvPatchScalarField
{}
multiphaseAlphaContactAngleFvPatchScalarField::
multiphaseAlphaContactAngleFvPatchScalarField
alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
(
const multiphaseAlphaContactAngleFvPatchScalarField& gcpsf,
const alphaContactAngleFvPatchScalarField& gcpsf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
@ -98,8 +95,7 @@ multiphaseAlphaContactAngleFvPatchScalarField
{}
multiphaseAlphaContactAngleFvPatchScalarField::
multiphaseAlphaContactAngleFvPatchScalarField
alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
@ -113,10 +109,9 @@ multiphaseAlphaContactAngleFvPatchScalarField
}
multiphaseAlphaContactAngleFvPatchScalarField::
multiphaseAlphaContactAngleFvPatchScalarField
alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
(
const multiphaseAlphaContactAngleFvPatchScalarField& gcpsf,
const alphaContactAngleFvPatchScalarField& gcpsf,
const DimensionedField<scalar, volMesh>& iF
)
:
@ -127,7 +122,7 @@ multiphaseAlphaContactAngleFvPatchScalarField
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void multiphaseAlphaContactAngleFvPatchScalarField::write(Ostream& os) const
void alphaContactAngleFvPatchScalarField::write(Ostream& os) const
{
fvPatchScalarField::write(os);
os.writeKeyword("thetaProperties")
@ -138,11 +133,7 @@ void multiphaseAlphaContactAngleFvPatchScalarField::write(Ostream& os) const
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
makePatchTypeField
(
fvPatchScalarField,
multiphaseAlphaContactAngleFvPatchScalarField
);
makePatchTypeField(fvPatchScalarField, alphaContactAngleFvPatchScalarField);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View File

@ -22,19 +22,19 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::multiphaseAlphaContactAngleFvPatchScalarField
Foam::alphaContactAngleFvPatchScalarField
Description
Contact-angle boundary condition for multi-phase interface-capturing
simulations. Used in conjuction with multiphaseMixture.
SourceFiles
multiphaseAlphaContactAngleFvPatchScalarField.C
alphaContactAngleFvPatchScalarField.C
\*---------------------------------------------------------------------------*/
#ifndef multiphaseAlphaContactAngleFvPatchScalarField_H
#define multiphaseAlphaContactAngleFvPatchScalarField_H
#ifndef alphaContactAngleFvPatchScalarField_H
#define alphaContactAngleFvPatchScalarField_H
#include "zeroGradientFvPatchFields.H"
#include "multiphaseMixture.H"
@ -45,10 +45,10 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class multiphaseAlphaContactAngleFvPatch Declaration
Class alphaContactAngleFvPatch Declaration
\*---------------------------------------------------------------------------*/
class multiphaseAlphaContactAngleFvPatchScalarField
class alphaContactAngleFvPatchScalarField
:
public zeroGradientFvPatchScalarField
{
@ -132,32 +132,31 @@ private:
public:
//- Runtime type information
TypeName("multiphaseAlphaContactAngle");
TypeName("alphaContactAngle");
// Constructors
//- Construct from patch and internal field
multiphaseAlphaContactAngleFvPatchScalarField
alphaContactAngleFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&
);
//- Construct from patch, internal field and dictionary
multiphaseAlphaContactAngleFvPatchScalarField
alphaContactAngleFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const dictionary&
);
//- Construct by mapping given
// multiphaseAlphaContactAngleFvPatchScalarField onto a new
// patch
multiphaseAlphaContactAngleFvPatchScalarField
//- Construct by mapping given alphaContactAngleFvPatchScalarField
// onto a new patch
alphaContactAngleFvPatchScalarField
(
const multiphaseAlphaContactAngleFvPatchScalarField&,
const alphaContactAngleFvPatchScalarField&,
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const fvPatchFieldMapper&
@ -168,14 +167,14 @@ public:
{
return tmp<fvPatchScalarField>
(
new multiphaseAlphaContactAngleFvPatchScalarField(*this)
new alphaContactAngleFvPatchScalarField(*this)
);
}
//- Construct as copy setting internal field reference
multiphaseAlphaContactAngleFvPatchScalarField
alphaContactAngleFvPatchScalarField
(
const multiphaseAlphaContactAngleFvPatchScalarField&,
const alphaContactAngleFvPatchScalarField&,
const DimensionedField<scalar, volMesh>&
);
@ -187,7 +186,7 @@ public:
{
return tmp<fvPatchScalarField>
(
new multiphaseAlphaContactAngleFvPatchScalarField(*this, iF)
new alphaContactAngleFvPatchScalarField(*this, iF)
);
}

View File

@ -24,11 +24,10 @@ License
\*---------------------------------------------------------------------------*/
#include "multiphaseMixture.H"
#include "multiphaseAlphaContactAngleFvPatchScalarField.H"
#include "alphaContactAngleFvPatchScalarField.H"
#include "Time.H"
#include "subCycle.H"
#include "fvCFD.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * //
@ -237,7 +236,7 @@ Foam::multiphaseMixture::surfaceTensionForce() const
}
void Foam::multiphaseMixture::correct()
void Foam::multiphaseMixture::solve()
{
forAllIter(PtrDictionary<phase>, phases_, iter)
{
@ -296,6 +295,10 @@ void Foam::multiphaseMixture::correct()
}
void Foam::multiphaseMixture::correct()
{}
Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
(
const volScalarField& alpha1,
@ -351,11 +354,10 @@ void Foam::multiphaseMixture::correctContactAngle
forAll(boundary, patchi)
{
if (isA<multiphaseAlphaContactAngleFvPatchScalarField>(gbf[patchi]))
if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
{
const multiphaseAlphaContactAngleFvPatchScalarField& acap =
refCast<const multiphaseAlphaContactAngleFvPatchScalarField>
(gbf[patchi]);
const alphaContactAngleFvPatchScalarField& acap =
refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
vectorField& nHatPatch = nHatb[patchi];
@ -363,7 +365,7 @@ void Foam::multiphaseMixture::correctContactAngle
mesh_.Sf().boundaryField()[patchi]
/mesh_.magSf().boundaryField()[patchi];
multiphaseAlphaContactAngleFvPatchScalarField::thetaPropsTable::
alphaContactAngleFvPatchScalarField::thetaPropsTable::
const_iterator tp =
acap.thetaProps().find(interfacePair(alpha1, alpha2));
@ -455,6 +457,34 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
}
Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixture::nearInterface() const
{
tmp<surfaceScalarField> tnearInt
(
new surfaceScalarField
(
IOobject
(
"nearInterface",
mesh_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar("nearInterface", dimless, 0.0)
)
);
forAllConstIter(PtrDictionary<phase>, phases_, iter)
{
surfaceScalarField alphaf = fvc::interpolate(iter());
tnearInt() = max(tnearInt(), pos(alphaf - 0.01)*pos(0.99 - alphaf));
}
return tnearInt;
}
void Foam::multiphaseMixture::solveAlphas
(
const label nAlphaCorr,
@ -466,7 +496,7 @@ void Foam::multiphaseMixture::solveAlphas
nSolves++;
word alphaScheme("div(phi,alpha)");
word alphacScheme("div(phic,alpha)");
word alphacScheme("div(phirb,alpha)");
tmp<fv::convectionScheme<scalar> > mvConvection
(

View File

@ -164,7 +164,7 @@ private:
multivariateSurfaceInterpolationScheme<scalar>::fieldTable alphaTable_;
// Private Member Functions
// Private member functions
void calcAlphas();
@ -256,6 +256,13 @@ public:
tmp<surfaceScalarField> surfaceTensionForce() const;
//- Indicator of the proximity of the interface
// Field values are 1 near and 0 away for the interface.
tmp<surfaceScalarField> nearInterface() const;
//- Solve for the mixture phase-fractions
void solve();
//- Correct the mixture properties
void correct();

Some files were not shown because too many files have changed in this diff Show More