From March to May 2023, CFD Direct completely replaced the liquid film functionality in the development version of OpenFOAM (OpenFOAM-dev). The new film functionality conserves mass, unlike its predecessor which was non-conservative and consequently inaccurate and unreliable for many problems. It is implemented with the solver module framework, which enables coupling to other regions with gas flows, multiphase flows, particle clouds, solids etc, including calculations of conjugate heat transfer (CHT). It is consistent with the rest of OpenFOAM so uses existing modelling libraries, e.g. thermodynamics, transport, etc. The new implementation contains 50% of the code lines of the original one, despite being more functional. It is therefore cheaper and easier to maintain, while being more extensible and robust.
The Original Implementation
The original film implementation was released with OpenFOAM v2.0.0 in 2011. It included a surfaceFilmModels library which included sub-models, e.g. for thermodynamics and transport, which were specific to films. Film calculations used a mesh, extruded one cell thick from a boundary surface using a extrudeToRegionMesh utility which accompanied the release.
A transport equation for the film thickness δ (deltain the code) described the film flow. The equation in volume integral form, i.e. starting ∫ [∂(ρδ)/∂t] dV, calculates the change in (cell mass × film height) per unit time. This is an extensive property, i.e. it is dependent on cell size, so the solution is meaningless unless the cell size is uniform. In that case, the equation represents mass conservation of the film (by dividing by the uniform cell height). However if the boundary surface is not flat, the extruded mesh will contain cells of non-uniform height, e.g. the height of a face will be typically greater than at the cell centre. The δ equation does not therefore conserve mass in general.
Making Film Conservative
The new film and isothermalFilm modules solve the flow of compressible liquid films, with the film module also calculating thermal energy. The isothermalFilm module contains the flow calculation (inherited by film). Instead of characterising the film by its height δ, it uses a volume fraction α to represent the volume of the cell occupied by the film, like phase fraction in the volume of fluid (VoF) method. The transport equation for α represents mass conservation of the film, irrespective of cell size, so works reliably for curved boundary surfaces. The approach is similar to the VoF method, but for a liquid layer whose height is not resolved. The film solution is even valid when α > 1. Note that the user configures a δ field for the film which the module converts to α for solving.
The film and cell dimensions can differ considerably. The cell dimensions are inappropriate for surface normal gradient calculations at film wall and surface boundaries. To overcome this problem, there are special filmWall and filmSurface patch types, and equivalent conditions mappedFilmWall and mappedFilmSurface for mapping between regions, e.g. for heat transfer (see $FOAM_APP/modules/isothermalFilm/patches). These conditions generate the patch deltas based on half the film height, rather than the cell dimensions to yield suitable surface normal gradient calculations.
Flow Coupling
A film occupies a region of mesh which can be coupled to other fluid and solid regions. Boundary conditions then calculate the interactions between regions. The filmSurfaceVelocity condition can be applied to the film at a boundary with a surrounding (bulk) fluid to evaluate its surface velocity from the shear imposed by that fluid. The mappedValue condition is applied on the fluid side to transfer the film motion to the fluid. The mappedFilmPressure condition maps the pressure from the fluid region to the film. These conditions enable each region to drive flow in the other, improving on the original implementation in which the fluid could drive the film, but not vice versa.
The film can include the surface tension and contact angle effects at a solid wall using the filmContactAngle boundary condition. The condition is applied to the δ field (which is converted to α) and accesses the same contact angle models as VoF or Euler-Euler multiphase simulations in OpenFOAM. A directionally-dependent contact angle makes rivulet flows chaotic without the need for random perturbation, which the original implementation required. The chaotic behaviour is shown (below, centre) in the rivuletBox tutorial example which couples three regions:
- a vertical aluminium panel (below, right), represented by the solid module;
- a film of water introduced at the top of panel (yellow region, below centre) with the film module;
- a thin air box (left) using the fluid module.
Multi-Region Modelling
The foamMultiRun application provides multi-region solutions using solver modules. This multi-region approach promotes code consistency, e.g. it uses standard boundary conditions to couple the regions. In the rivuletBox example heat transfer is simply managed through the rationalised coupledTemperature condition for temperature. The image above, right, shows the temperature in panel, initially at 330K, which is cooled by the film entering at 300K. While the heat transfer worked correctly for films without any bespoke code, the original implementation incorrectly calculated a different heat flux from a bulk fluid to film and vice versa.
The inherent consistency means that very little additional code is required beyond the core film functionality. Code reuse extends to large, existing modelling libraries including momentum transport, thermal transport and thermodynamics. This is in contrast with the original implementation which duplicated models, e.g. the thixotropic non-Newtonian model. Inconsistency makes the code harder to extend and inevitably less robust and more difficult to use.
Other Multiphase Coupling
Further coupling is available between films and other multiphase modelling in OpenFOAM. Firstly, the filmVoFTransfer and VoFFilmTransfer fvModels can be applied to a film or VoF region, respectively, to transfer a film to the VoF phase when the film exceeds a specified thickness and, vice versa, when a VoF film becomes thin. Particles can transfer into a film by specifying the cloudFilmTransfer surfaceFilmModel in the fluid cloudProperties file and applying the filmCloudTransfer fvModel to the film. The filmCloudTransfer fvModel includes an optional ejection model can transfer fluid from a film to particles (parcels) by dripping from an inverted surface or curvature separation.
Conservation and Consistency
CFD Direct maintains OpenFOAM to meet the critical needs of users which are: availability; usability; extensibility and robustness. Code redesign is one aspect of maintenance which aims to correct flaws and improve code structure (targeting robustness and extensibility, respectively). Since 2020, we have redesigned and/or created the following components of OpenFOAM.
- Solver modules: general CFD solution for single and multiple regions.
- Dynamic meshes: mesh motion, topology change (e.g. cell refinement), parallel decomposition and mesh-to-mesh mapping.
- Non-conformal coupling: connecting domain regions, especially with rotating geometry.
- Models and constraints: general interface for models (e.g. particles, thermal radiation, heat sources) and numerical constraints.
- Transport modelling: consistent framework for momentum and thermophysical transport.
Our work adheres to two key principles: 1) implementing conservative numerics for accurate, robust solution; 2) demanding consistency between components of OpenFOAM. The film redesign demonstrates this clearly. The new implementation is conservative when the original one was not. We used existing structures, e.g. a standard mesh, finite volume numerics, boundary conditions and model libraries, within a generalised multi-region framework, to implement this functionality on a reasonable budget in two months. The new implementation contains 50% of the code lines of the old one, making it easier to maintain. We get it, right.