openfoam merge requestshttps://develop.openfoam.com/Development/openfoam/-/merge_requests2021-07-16T16:39:50Zhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/475finiteArea: new models and solvers for films2021-07-16T16:39:50ZSergio FerrarisfiniteArea: new models and solvers for films- Adding `kinematicParcelFoam` solver:
- The original `thermoSurfaceFilm` sub-models were divided between `kinematicSurfaceFilm` and `thermoSurfaceFilm` in order to use the `surfaceFilm` model in a `kinematicCloud`.
- The film in...- Adding `kinematicParcelFoam` solver:
- The original `thermoSurfaceFilm` sub-models were divided between `kinematicSurfaceFilm` and `thermoSurfaceFilm` in order to use the `surfaceFilm` model in a `kinematicCloud`.
- The film interaction models are now in a `kinematicSurface` class which can be used in a kinematic cloud adding constant thermal properties (`p` and `T`) for some sub-models, e.g. `drySplashInteraction`, `wetSplashInteraction` etc.
- `pRef` and `Tref` were added to the `kinematicSurfaceFilm` as entry to the `regionFilm` when used with a kinematic cloud.
- In the finite area surface film model `Tref`, `pRef` are stored in `filmSubModel`.
- Film models:
- Only laminar turbulence is available for wall friction and gas shear stress. Wall friction models:
- quadraticProfile,
- linearProfile,
- DarcyWeisbach,
- ManningStrickle
- gas friction models: `Cf *(Ufilm - Ufluid)`. This would need to be improved using turbulent shear stress from the flow.
- `curvatureSepration` injection model from the film to Lagrangian particles
- minimum (`fThreshold`) force and minimum curvature (`minInvR1`) for separation were added to have more control on determining the points of film separation.
- To force wall-like BC from the fluid side the option `zeroWallVelocity true;`. Allows the film to be seen as zero U BC by the flow.
- forces: The only available extra forces (apart from gravity) is contact angle force: `perturbedTemperatureDependentContactAngle`v2112Sergio FerrarisSergio Ferrarishttps://develop.openfoam.com/Development/openfoam/-/merge_requests/476ENH: turbulenceFields: various improvements2021-08-04T11:37:18ZKutalmış BerçinENH: turbulenceFields: various improvements- 55f61a0037 - ENH: turbulenceFields: enable custom prefix for output fields
- 8f16527a51 - BUG: turbulenceFields: unset duplicate already-registered fields
- 86f753b7bc - DOC: turbulenceFields: improve header documentation
- 189053421c...- 55f61a0037 - ENH: turbulenceFields: enable custom prefix for output fields
- 8f16527a51 - BUG: turbulenceFields: unset duplicate already-registered fields
- 86f753b7bc - DOC: turbulenceFields: improve header documentation
- 189053421c - STYLE: turbulenceFields: apply more recent OpenFOAM-coding practices
- 4243291c52 - BUG: turbulenceFields: use omega funcs of turbulence models (fixes #2132)
#### Resolved bugs
#2132
#### Test cases
[steadyBoundaryLayer](https://tinyurl.com/5fbtx82v)v2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/477Implicit treatment of coupled boundary conditions2021-08-03T20:08:53ZSergio FerrarisImplicit treatment of coupled boundary conditions## Summary
The alternative addressing insert new internal faces on selected coupled patches. Extending the matrix addressing allows the implicit treatment of some boundary conditions. At the moment, the BCs which can be made implicit ar...## Summary
The alternative addressing insert new internal faces on selected coupled patches. Extending the matrix addressing allows the implicit treatment of some boundary conditions. At the moment, the BCs which can be made implicit are: cyclic, AMI, ACMI and mapped.
The entry `useImplicit true` is needed in the field BC. Only scalar fields can be made implicit (p, p_rgh, k, epsilon, etc). The U field **cannot** be implicit.
The top solvers are unchanged except for the cht solvers where a single he matrix is built for the multi-region case.
In the case that more than one patch is present, the user can choose which ones are treated implicitly, i.e **p** can be implicit in AMI1 and explicit in AMI2 and **k** the other way around.
Currently, there exists a limitation for parallel cases where AMI patch-pairs need to be on the same processor. For cyclics and mapped patches, the number of faces on each processor **MUST** be the same. Therefore a constraint decomposition needs to be used.
## Details of new models
The model adds a new numbering of internal/boundary coefficients and patch ID's. This is constructed at the call of the matrix.solve() - it is reset to the original addressing after solving.
The original internal/boundary coeffs are cached and used (in some cases with some extra information) to fill the coefficients on the new internal faces and/or source term on the newly formed fvMatrix.
## Tutorials
- tutorials/heatTransfer/chtMultiRegionSimpleFoam/cpuCabinet/
- tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeaterImplicit/
- tutorials/basic/laplacianFoam/implicitAMI/
- tutorials/basic/simpleFoam/implicitAMI
- tutorials/basic/chtMultiRegionFoam/2DImplicitCyclic/
## Risks
Addressing of cells on the patches are changed, this needs to be passed to all the lduInterfaces. Use alternative lduAddressing in lduMatrix.
The coupled BCs (cyclics, AMI, mapped) which made implicit don't exist at the moment of solving the matrix. Thus, any fvOption or FO referring to this BC will **Fail**
In general terms when using the implicit approach the lduaddressing changes and therefore any matrix solve will be affected in _all top level solvers and all the linear solvers_.v2112Sergio FerrarisSergio Ferrarishttps://develop.openfoam.com/Development/openfoam/-/merge_requests/478ENH: linear solvers: add variable-specific debug flags2021-08-04T09:15:33ZKutalmış BerçinENH: linear solvers: add variable-specific debug flags### Summary
* a762a311 - ENH: linear solvers: add variable-specific debug flags <Kutalmis Bercin>
* 390e7b0a - TUT: rotatingFanInRoom: perturb locationInMesh (fixes #2162) <Matej Forman>
* 7d7e012e - ENH: KirchhoffShell: simplification...### Summary
* a762a311 - ENH: linear solvers: add variable-specific debug flags <Kutalmis Bercin>
* 390e7b0a - TUT: rotatingFanInRoom: perturb locationInMesh (fixes #2162) <Matej Forman>
* 7d7e012e - ENH: KirchhoffShell: simplification of log output <Kutalmis Bercin>
Introduces a new optional keyword of label type 'log'
to linear-solver dictionaries to enable variable-specific
debug statements. For example, in fvOptions file:
solvers
{
p
{
solver GAMG;
...
log 2;
}
U
{
...
log 0;
}
}
The meanings of values of 'log' are:
log 0; <!-- no output
log 1; <!-- standard output
log 2; <!-- debug output
// values higher than 2 are expected to have no effect
### Resolved bugs (If applicable)
#2162
EP#1615 @Chiara
### Risks
This keyword does not directly affect the operations of various `DebugSwitches` and
backward compatibility has been ensured in exchange of code cleanness. The related `DebugSwitches` are:
DebugSwitches
{
SolverPerformance 0;
GAMG 0;
PCG 0;
PBiCG 0;
smoothSolver 0;
}
### Test cases
[4-ep-1615-solver-verbosity](https://tinyurl.com/4wrp4hj4)v2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/480improve handling of multiple world communication2021-08-20T14:10:33ZMark OLESENimprove handling of multiple world communication- when attempting to connect several worlds, the overlapping communicators (or other parts) cause the simulation to bloc.- when attempting to connect several worlds, the overlapping communicators (or other parts) cause the simulation to bloc.v2112Mattijs Janssens4-Mattijs@users.noreply.develop.openfoam.comMattijs Janssens4-Mattijs@users.noreply.develop.openfoam.comhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/483ENH: interpolate in cylindrical coordinates: See #21452021-09-16T08:48:17ZMattijs Janssens4-Mattijs@users.noreply.develop.openfoam.comENH: interpolate in cylindrical coordinates: See #2145This intercepts all vector/tensor AMI interpolations and
does the interpolation in cylindrical coordinates.
Use (component) 'coupled' linear solver to enable this in
the linear solver sweeps as well (however this is not very stable
sinc...This intercepts all vector/tensor AMI interpolations and
does the interpolation in cylindrical coordinates.
Use (component) 'coupled' linear solver to enable this in
the linear solver sweeps as well (however this is not very stable
since the vector solution is still in Cartesian coordinates -
only the AMI interpolation is in cylindrical)v2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/484BUG: fix yPlus predictions of nut wall functions2021-09-14T15:21:51ZKutalmış BerçinBUG: fix yPlus predictions of nut wall functions### Summary
- Corrects `y+` predictions from `nutUWallFunction`, `nutkWallFunction` and `nutLowReWallFunction`, but the internal `y+` predictions of these wall functions were kept the same in order to keep their `nut` predictions the sa...### Summary
- Corrects `y+` predictions from `nutUWallFunction`, `nutkWallFunction` and `nutLowReWallFunction`, but the internal `y+` predictions of these wall functions were kept the same in order to keep their `nut` predictions the same.
### Resolved bugs
Closes #1773, #1838, #1846, #1843
### Risks
- No changes in wall-function output.
- No changes to user input.
- New optional keyword into the `yPlus` function object: `useWallFunction`.
- Changes in the `yPlus` function object when the above nut wall functions are being used.
### Tests
- [x] Compilation: `linux64ClangDPInt32Opt`, `linux64GccDPInt32Opt`, `linux64GccSPDPInt64Debug`
- [x] Alltest: No change
- [x] [simpleFoam & rhoSimpleFoam](https://tinyurl.com/jm3n3e)
- [x] Verification: [planeChannel](https://tinyurl.com/3k3z8utc)
Verification case characteristics:
* One-dimensional smooth-wall plane channel flow, ReTau=5200
* Number of cells = 20
* nu = 0.000192827 \[m2/s\]
* The case was designed to produce => -1\*sqrt(mag(wallShearStress)) = 1
* The y1+ set (expected) = {0.05, 0.5, 1, 5, 10, 20, 30, 50, 100}
* kOmegaSST & simpleFoam
Fields of interest:
* `mag(wallShearStress)` function object correctly returns \~ 1 for all test cases.
* `yPlus` function object (min) results for `lowerWall` rounded up to 2 decimals can be seen below:
### y1+ = 0.05:
| Wall function | Pre-fix yPlus | Post-fix yPlus |
| --- | ------ |---------:|
| nutUWallFunction | **0.12** | 0.05 |
| nutkWallFunction | **0.00025** | 0.05 |
| nutUSpaldingWallFunction | 0.05 | 0.05 |
| nutUBlendedWallFunction | 0.05 | 0.05 |
| nutLowReWallFunction | 0.05 | 0.05 |
### y1+ = 0.5:
| Wall function | Pre-fix yPlus | Post-fix yPlus |
| --- | ------ |---------:|
| nutUWallFunction | **0.18** | 0.50 |
| nutkWallFunction | **0.025** | 0.50 |
| nutUSpaldingWallFunction | 0.50 | 0.50 |
| nutUBlendedWallFunction | 0.50 | 0.50 |
### y1+ = 1:
| Wall function | Pre-fix yPlus | Post-fix yPlus |
| --- | ------ |---------:|
| nutUWallFunction | **0.34** | 1.00 |
| nutkWallFunction | **0.1** | 1.00 |
| nutUSpaldingWallFunction | 1.00 | 1.00 |
| nutUBlendedWallFunction | 1.00 | 1.00 |
### y1+ = 5:
| Wall function | Pre-fix yPlus | Post-fix yPlus |
| --- | ------ |---------:|
| nutUWallFunction | **3.02** | 5.00 |
| nutkWallFunction | **2.75** | 5.00 |
| nutUSpaldingWallFunction | 5.00 | 5.00 |
| nutUBlendedWallFunction | 5.00 | 5.00 |
### y1+ = 10:
| Wall function | Pre-fix yPlus | Post-fix yPlus |
| --- | ------ |---------:|
| nutUWallFunction | **9.12** | 10.00 |
| nutkWallFunction | **8.31** | 10.00 |
| nutUSpaldingWallFunction | 10.00 | 10.00 |
| nutUBlendedWallFunction | 10.00 | 10.00 |
### y1+ = 20:
| Wall function | Pre-fix yPlus | Post-fix yPlus |
| --- | ------ |---------:|
| nutUWallFunction | 20.00 | 20.02 |
| nutkWallFunction | **19.30** | 20.01 |
| nutUSpaldingWallFunction | 20.00 | 20.00 |
| nutUBlendedWallFunction | 20.00 | 20.01 |
### y1+ = 30:
| Wall function | Pre-fix yPlus | Post-fix yPlus |
| --- | ------ |---------:|
| nutUWallFunction | 30.09 | 30.09 |
| nutkWallFunction | 29.91 | 30.12 |
| nutUSpaldingWallFunction | 30.08 | 30.08 |
| nutUBlendedWallFunction | 30.09 | 30.09 |
### y1+ = 50:
| Wall function | Pre-fix yPlus | Post-fix yPlus |
| --- | ------ |---------:|
| nutUWallFunction | 50.16 | 50.16 |
| nutkWallFunction | 50.51 | 50.19 |
| nutUSpaldingWallFunction | 50.17 | 50.17 |
| nutUBlendedWallFunction | 50.16 | 50.16 |
### y1+ = 100:
| Wall function | Pre-fix yPlus | Post-fix yPlus |
| --- | ------ |---------:|
| nutUWallFunction | 100.01 | 100.01 |
| nutkWallFunction | 100.96 | 100.01 |
| nutUSpaldingWallFunction | 100.01 | 100.01 |
| nutUBlendedWallFunction | 100.01 | 100.01 |v2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/488add point snapping to iso-surface topo algorithm (#2210)2022-12-12T16:16:46ZMark OLESENadd point snapping to iso-surface topo algorithm (#2210)- helps avoid the creation of small face cuts (near corners, edges)
that result in zero-size faces on output.- helps avoid the creation of small face cuts (near corners, edges)
that result in zero-size faces on output.v2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/493ENH: new RANS model for geophysical applications2021-11-08T18:58:08ZKutalmış BerçinENH: new RANS model for geophysical applications### Summary
- ENH: kL: new RANS model for geophysical applications
- STYLE: atmosphericModels: regroup atmospheric turbulence models
- ENH: kEpsilonLopesdaCosta: enable this RANS model for compressible applications
- TUT: atmFlatTer...### Summary
- ENH: kL: new RANS model for geophysical applications
- STYLE: atmosphericModels: regroup atmospheric turbulence models
- ENH: kEpsilonLopesdaCosta: enable this RANS model for compressible applications
- TUT: atmFlatTerrain: add an example for kL RANS model
- BUG: atmFlatTerrain: fix input settings in the plot script
### Resolved bugs (If applicable)
- A small issue in `atmFlatTerrain/precursor/plot` has been corrected.
### Details of new models (If applicable)
- For a given case, the `kL` model reduces numerical costs approximately by one-third in comparison to the `kEpsilon` model for geophysical applications.
- A set of results obtained from the `atmFlatTerrain` and `atmForestStability` tutorials (left: kEpsilon, right: kL):
![image](/uploads/79d0f1db91983844da41db3af79599b6/image.png)
![image](/uploads/e1190ba19b50ddc5a093660d7f2ee4cd/image.png)
![image](/uploads/d6999f8f9d833e90acaf0a941a428c1d/image.png)
![image](/uploads/c8fbc676998046df8f8498e935ba1f21/image.png)
![image](/uploads/03df4d3859b04af634bb2c05d7391080/image.png)
![image](/uploads/d6c34c1c9e2144096b73363e8f12297c/image.png)
### Risks
- No changes in user input.
- No changes are expected in output.
- The directory structure of `src/atmosphericModels` has been changed. Therefore, the location of the directory `kEpsilonLopesdaCosta` has also been changed. Might affect few minor things for a few developers. Not crucial, IMHO.
### Tests
* [x] Compilation (incl. submodules):
* [x] `linux64ClangDPInt32Opt` (clang11)
* [x] `linux64GccDPInt32Opt`
* [x] `linux64GccSPDPInt64Debug`
* [x] Alltest: No change in output with respect to the develop HEAD + no errorv2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/496ENH: New suite for electrostatic deposition applications2021-11-24T18:06:40ZKutalmış BerçinENH: New suite for electrostatic deposition applications### Summary
- New suite for electrostatic deposition applications
- Solver function object: `electricPotential`: Computes the steady-state equation of charge conservation to obtain the electric potential by strictly assuming a quasi-...### Summary
- New suite for electrostatic deposition applications
- Solver function object: `electricPotential`: Computes the steady-state equation of charge conservation to obtain the electric potential by strictly assuming a quasi-static electrostatic field for single-phase and multiphase applications.
- Finite-volume boundary condition: `electrostaticDeposition`: A boundary condition to calculate electric potential (`V`) on a given boundary based on film thickness (`h`) and film resistance (`R`) fields which are updated based on a given patch-normal current density field (`jn`), Coulombic efficiency and film resistivity.
- Allowing to set:
- Minimum current density for deposition onset
- Minimum accumulative specific charge for deposition onset
- Resistance due to main body and/or pretreatment layers
- Initial electric potential
### Resolved bugs
N/A
### Details of new models (If applicable)
A showcase illustrating the verification of the `electricPotential` function object. Right: dielectric-dielectric water-air; left: conductive-dielectric water-air mixture ([López-Herrera et al., 2011](https://www.doi.org/10.1016/j.jcp.2010.11.042)).
![image](/uploads/fd8508ee81de44186bcc26d60eebc0f3/image.png)
A showcase illustrating the `electrostaticDeposition` BC: pending
### Risks
- No changes in existing output.
### Constraints
- Depletion or abrasion of material due to negative current is not allowed.
- Boundary-condition updates are not allowed during outer corrections to prevent spurious accumulation of film thickness.
- `resistivity`, `jMin`, `qMin` and `Rbody` are always non-negative.
- Quasi steady
- No electrolyte/ionic transport
- Quiescent
- Sufficiently mixed
- No magnetic field
- No electrolyte and electric field interaction
- Isotropic material and electrolyte
- No film-flow interactions
- No finite-area numerics
### Remaining issues
- Micro/nano current densities/potential differences cause numerical instabilities.
- Small diffs in patch-normal current-density values of cathode faces can lead to small diffs in incremental film thickness (Delta_h) on those faces. With the very high value of 'resistivity', the small diffs in 'Delta_h' can then lead to comparable diffs in incremental electric potential differences across patch faces.
- These differences can then be accumulated throughout a simulation, causing non-zero potential difference gradients across faces of the cathode.
- Our workaround was to round the floating-point numbers of operand variables to a given number of decimal points (default: 8 decimals) - disallow micro/nano volts, so that such accumulation could be avoided.
- However, this issue may cause numerical instabilities in more complex cases.
- Therefore, the original expressions of the formulation may need to be revisited.
- Advanced electrodeposition applications using OpenFOAM can be visited as well:
- [Kauffman et al., 2020](https://doi.org/10.3390/fluids5040240)
- [López-Herrera et al., 2011](https://www.doi.org/10.1016/j.jcp.2010.11.042)
- [Roghair et al., 2013](https://ris.utwente.nl/ws/portalfiles/portal/5676495/Roghair+Paper_final_draft_v1.pdf)
- Takayama et al., 2013, Implementation of a Model for Electroplating in
OpenFOAM.
- [Kashir et al., 2019](https://www.doi.org/10.1063/1.5050793)
- [EHDIonFOAM](https://openfoamwiki.net/index.php/Extend-bazaar/solvers/multiphase/EHDIonFOAM)v2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/499Function1 objectRegistry access2022-02-25T06:53:57ZAndrew HeatherFunction1 objectRegistry access### Summary
Cross reference [EP 1657](https://exchange.openfoam.com/node/1657)
### Details of new models
Function1's
- `Function1` : The Function1 class can now be created using an `objectRegistry`, i.e. mesh or time databases, allow...### Summary
Cross reference [EP 1657](https://exchange.openfoam.com/node/1657)
### Details of new models
Function1's
- `Function1` : The Function1 class can now be created using an `objectRegistry`, i.e. mesh or time databases, allowing other registered objects to be retrieved, e.g. fields, models, function object state and results information
- change propagated across all `Function1` instances, e.g. boundary conditions, function objects, other `Function1s`
- `functionObjectValue` : set the `value` according to the result of a function object
- `sample` : set the `value` according to a field value
- `limitRange` : renamed to `inputValueMapper` and extended to include operations using function object results
Function object updates
- `sampledSets` : extended to store the set `min`, `max` and `average` values in its results dictionary
- `reference` : the reference value `refValue` is now set according to a `Function1` type. To recover the old behaviour (sampling a field value) users can employ the new `sample` Function1 type.
- registered fields are now prefixed - typically by `<function_object_name>:` - to avoid name collisions on registration, e.g. you can now have multiple `fieldAverage` function objects that act on the same field
Test case
- `$FOAM_TUTORIALS/incompressible/simpleFoam/simpleCar`
- Sample the pressure field using a `sets` function object; the result `average(p)` is created.
- The sample `average(p)` value is supplied to a `valueAverage` function object to perform time averaging. This creates the result `average(p)Mean`.
- The `reference` function object's `refValue` makes use of the new `functionObjectValue` Function1, where the value `average(p)Mean` is retrieved from the `average1` function object to determine the final reference field.
- A second `reference` function object to shows how to recover the old function object behaviour.
### Risks
Testing required
- `TimeFunction1` has been deprecated
- check that previous behaviours are maintained, e.g. for lagrangian injection models when using `engineTime`v2112Mark OLESENMark OLESENhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/500ENH: New gradient scheme iterativeGaussGrad2021-12-08T11:05:44ZKutalmış BerçinENH: New gradient scheme iterativeGaussGrad### Summary
**Aim:**
Improvement of numerical **fidelity** only via improving **gradient** computations for meshes with **non-zero face skewness**.
**(Theory) What is face skewness:**
- Level of closeness of the following two spatial...### Summary
**Aim:**
Improvement of numerical **fidelity** only via improving **gradient** computations for meshes with **non-zero face skewness**.
**(Theory) What is face skewness:**
- Level of closeness of the following two spatial points:
- Face intersection of the PN vector between cell centres of an owner cell and a neighbour cell
- Face centre of the common face of the two cells
**(Theory) Why is face skewness important:**
- Prediction fidelity
- Any deviation from face centre as the face average centre reduces the second-order accuracy (i.e. non-zero face skewness)
- Hence OpenFOAM has various skewness corrections for centre-to-face interpolations
- Numerical stability
- Reduces the diagonal dominance of the discrete Poisson operator which slows down convergence. ADI (i.e. alternating-direction implicit) type preconditioners also become less effective.
**(Theory) What is the technical challenge in face skewness:**
- Omission wherever possible due to cost concerns rather than a true technical challenge
- Gradient computations are accounted for ~20% of a typical simulation
**(Theory) Treatments in OpenFOAM:**
- Many methods to quantify skewness
- OpenFOAM has its own definition
- Skewness formulations in `checkMesh` and `snappyHexMesh` are a bit different
- Various skewness corrections available
- `skewCorrectedSnGrad` surface-normal gradient scheme
- `skewCorrected` surface interpolation scheme
- No corrections on boundary skewness
- No corrections in Gauss-Green gradient computations
- Not needed for `leastSquares` and `pointLinear` based gradient schemes
### Resolved bugs
N/A
### Methodology
**Problem definition:**
- Face skewness is not accounted during any Gauss-Green gradient computations.
- Incorporation of face skewness into gradient computations **may or may not** improve the level of prediction fidelity.
**A problem solution: Iterative Gauss gradient**
- Add an outer loop around the gradient and skew-correction vector computations
- Add an optional relaxation factor inside the gradient computation
<img src="/uploads/d8930fd27b2515bb0cfd25814927918e/image.png" width="50%" height="50%">
**Test case:**
- Very difficult to design a test case where skewness is isolated from non-orthogonality
- [Manufactured solution](https://www.openfoam.com/documentation/guides/latest/doc/guide-schemes-gradient-example.html): A two-dimensional triangle-cell domain where theoretical gradient is square-root of 2 in each co-ordinate direction
**Metrics:**
- Numerical stability
- Simulation crashes or not
- Prediction fidelity - Level of discrepancy between theoretical and numerical values:
<img src="/uploads/8eb4d34c480bf6683a8814a24938c7b2/image.png" width="50%" height="50%">
- Arithmetic average of error field
- Coefficient of variation (CoV) of error field
- Relative standard deviation: std/mean
- Measure of amount of dispersion
- Low CoV: Values are close to mean
- High CoV: Values spread out over a wider range, potentially indicating outliers with respect to the mean
**Control variable:**
- Skewness is not the control variable since changing skewness also changes non-orthogonality in general
- Therefore:
- We fixed skewness and non-orthogonality with a fixed-geometry test case
- Control variable becomes the gradient scheme itself
- We quantified the performance of each gradient scheme with the chosen metric and compared them with each other
**Mesh metrics:**
![image](/uploads/69778f9a3803614e35436dd8aac054a8/image.png)
### Results
![image](/uploads/26960502a9e3d24f21e98798ab070e4b/image.png)
![image](/uploads/fe29294c40d765384b930115275d9259/image.png)
![image](/uploads/ecb33932f656ded148e1ebd2afa2218f/image.png)
![image](/uploads/5f80d7f935a35a84b5da51a4484601fe/image.png)
![image](/uploads/11ca5223abb62764749ed71a5b01dce7/image.png)
![image](/uploads/836727a0db9861dac012a9caa9355f12/image.png)
![image](/uploads/127620e5dbc5bc140a2d888e2c1419f1/image.png)
### Discussion
**What happened?**
- New iterative Gauss scheme
- Reduced errors in internal fields in comparison to other Gauss schemes
- Reduced errors in boundary fields in comparison to least square schemes
- Increasing number of iterations monotonically reduced errors
- Application of interpolation-scheme limiters increased errors
**What do the results mean in practise? (How will it change in how we do things?)**
- The new scheme **may** be used to improve prediction fidelity for gradient computations on meshes with skewness and non-orthogonality
- More tests are needed for:
- Effects of extreme levels of face skewness
- Cost estimations
- The results **do not suggest** that the new scheme can improve numerical stability
### Risks
- No changes in existing output.
- No changes in existing user input.
### Constraints
- No finite-area
- Needs extensive tests to assume safe for single-phase, multiphase and overset finite-volume applications
- No boundary treatments for Dirichlet and Neumann conditions
### Tests
* [x] Compilation (incl. submodules):
* [x] `linux64ClangDPInt32Opt` (clang11)
* [x] `linux64GccDPInt32Opt`
* [x] `linux64GccSPDPInt64Debug`
* [x] Alltest: No change in output with respect to the develop HEAD + no errorv2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/502ENH: New schemes for non-orthogonal meshes2021-12-08T12:08:02ZKutalmış BerçinENH: New schemes for non-orthogonal meshes### Summary
**Aim:**
Improvement of numerical **stability** of **non-orthogonal** mesh simulations **only** via changes in the **diffusion operator** handling.
**(Theory) What is face non-orthogonality:**
- "The angle between the ce...### Summary
**Aim:**
Improvement of numerical **stability** of **non-orthogonal** mesh simulations **only** via changes in the **diffusion operator** handling.
**(Theory) What is face non-orthogonality:**
- "The angle between the centroid vector and the unit normal vector is the orthogonality angle." (Wimhurst, 2019)
**(Theory) Why is face non-orthogonality important:**
- Can affect the discrete diffusion operator:
- Numerical stability (hence cost)
- Level of errors
- Numerical stability
- Larger explicit treatment -> less stability
- Level of errors
- Larger corrections -> larger truncation errors
- Cannot be reduced by mesh refinement
**(Theory) What is the technical challenge in face non-orthogonality:**
- "The difficulty is to evaluate the dot product of the normal vector and the velocity gradient on the cell face." (Wimhurst, 2019)
**(Theory) Treatments in OpenFOAM:**
- Explicit correction by using an over-relaxed-type unit-normal vector decomposition
### Resolved bugs
N/A
### Methodology
**Problem definition:**
- Convergence rate
- Explicit component can become limiting factor for numerical stability
- During initial iterations, especially for steady-state runs or when starting from uniform fields in transient runs
- Yet its removal reduces accuracy
- Convergence order
- Non-orthogonality treatments are neglected at boundaries (except cyclic conditions)
- Even minor non-orthogonality on patches reduce the convergence order for cases where Dirichlet/Neumann boundary conditions are applied. (Noriega et al., 2018)
**A solution: Field relaxation**
- Field relaxation for explicit non-orthogonality components
- Under-relaxation for numerical stability
- Over-relaxation for numerical cost
- Implementation: a Laplacian scheme
- Based on an idea from `nextFoam`
- Named `relaxedNonOrthoLaplacian`
- Implementation: a surface-normal gradient scheme:
- Named `relaxedSnGrad`
**Test case:**
- The DNS study of a smooth-wall plane channel flow by (Moser et al., 1999) where the friction Reynolds number of the flow is ReTau=395.
**Metrics:**
- Numerical stability
- Simulation crashes or not
- Prediction fidelity
- Flow field predictions vs DNS statistics
**Control variable:**
<img src="/uploads/5d6921a672ee7d70d6672d48505861d5/image.png" width="50%" height="50%">
### Results
![image](/uploads/5b997d95e965d8a895b387a6d37f762d/image.png)
![image](/uploads/496cedf11adafb1a9933d0d10deca83d/image.png)
![image](/uploads/cb03fc29aece34350d7576a1c21126e6/image.png)
![image](/uploads/f3cbd51e47723f982277520b42517794/image.png)
![image](/uploads/57b268055f93e271c3766fd78be95290/image.png)
![image](/uploads/e50faee3fb0988e7f557fd5302cbf299/image.png)
### Discussion
**What happened?**
- New schemes
- Did not dampen initial residuals unlike nextFoam observations (at least for this particular case)
- Improved prediction fidelity for above 50 degrees
- Prevented simulations crashing after a certain non-orthogonality angle (above 50)
- Numerical cost
- For below-50-degree non-orthogonality cases, the runtime remained similar to that of the base
- For above-50-degree cases, the numerical cost seems to be reduced
<img src="/uploads/59a1e422f5be56ac6e631e60c7b82e22/image.png" width="50%" height="50%">
**What do the results mean in practise? (How will it change in how we do things?)**
- New schemes **may** be used
- to avoid any simulation crashes due to non-orthogonality treatments
- without affecting the fidelity of numerical predictions
- without affecting numerical cost estimations
- more tests are needed
### Risks
- No changes in existing output.
- No changes in existing user input.
### Constraints
- No finite-area numerics
- Needs extensive tests to assume safe for single-phase, multiphase and overset finite-volume applications
- No boundary treatments for Dirichlet and Neumann conditions
### Tests
* [x] Compilation (incl. submodules):
* [x] `linux64ClangDPInt32Opt` (clang11)
* [x] `linux64GccDPInt32Opt`
* [x] `linux64GccSPDPInt64Debug`
* [x] Alltest: No change in output with respect to the develop HEAD + no error
##### References
- Wimhurst, A. (2019). Mesh Non-Orthogonality. https://tinyurl.com/49b55tzw
- Wimhurst, A. (2019). Mesh Non-Orthogonality 2: The Over-Relaxed Approach https://tinyurl.com/234h7hbt
- Noriega et al. (2018). A case-study in open-source CFD code verification, Part I: Convergence rate loss diagnosis. DOI: 10.1016/j.matcom.2017.12.002v2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/503sampledSets - enable writer construction from dictionary; glTF export2021-12-06T13:14:10ZAndrew HeathersampledSets - enable writer construction from dictionary; glTF export## sampledSets - added writer construction from dictionary
The `sampledSets` `writer` is now constructed with a dictionary to enable custom options using a new `formatOptions` dictionary - see example below.
## glTF format
glTF is a c...## sampledSets - added writer construction from dictionary
The `sampledSets` `writer` is now constructed with a dictionary to enable custom options using a new `formatOptions` dictionary - see example below.
## glTF format
glTF is a common data exchange format for virtual reality software, e.g. [ESI's IC.IDO](https://www.esi-group.com/products/virtual-reality). Scenes are described using a JSON representation and the field/binary data is stored in an external file.
These changes add [glTF v2](https://www.khronos.org/registry/glTF/) as an option when writing `sampledSets` and particle tracks.
Note: glTF files can also be read in ParaView v5.9
### Additions:
- The `particleTracks` utility:
- output path has been updated to `postProcessing/lagrangian/<cloudNme>/particleTracks/`;
- has a new option to limit the number of tracks exported; and
- supports field export (previously only geometry was written).
- example:
cloud reactingCloud1;
sampleFrequency 1;
maxPositions 1000000;
//maxTracks 5; // optional limit on the number of tracks
setFormat <format>;
// Optional field specification; wildcards supported
// - if ommitted all fields are written
fields (d T);
- glTF format support added to sampledSets, e.g.
type sets;
libs (sampling);
interpolationScheme cellPointFace;
setFormat gltf;
### glTF format options
Controls are exposed via the `formatOptions` dictionary, e.g.
setFormat gltf;
formatOptions
{
...
}
Options include:
- Adding colour (RGBA):
formatOptions
{
colour yes;
}
By default, this will add colour fields named `COLOR_0`, `COLOR_1`, ... `COLOR_N` corresponding to each value field, using the default colour map (coolToWarm), anchoring the colour map limits to the field minimum and maximum.
Individual field specification is enabled using a `fieldInfo` sub-dictionary, e.g.
formatOptions
{
colour yes;
fieldInfo
{
T // Field name
{
// Optional colour map; default = coolToWarm
colourMap fire;
// Alpha channel description
alpha field; // uniform | field;
//alphaValue 0.1; // uniform alpha value
alphaField T;
normalise yes;
}
}
}
Note that when using the `field` option for alpha channel scaling, the selected field should be the same type as the field being output, i.e. scalar, vector etc. since the fields are grouped into their primitive types when output.
![Screenshot_from_2021-11-26_14-40-53](/uploads/28f5fd037e102b9887a38c08d59cf5eb/Screenshot_from_2021-11-26_14-40-53.png)
Particle tracks have some additional options:
- Static tracks can be written as a set of points along the track (lines support may be added in the future). Specification is the same as above
![Screenshot_from_2021-11-26_14-40-27](/uploads/a3f687dd4ccd9d5cca5e80e830c0fdef/Screenshot_from_2021-11-26_14-40-27.png)
- Tracks can be animated, i.e. seeding a particle at the start of the track and then updating its position as a function of time:
```
formatOptions
{
animate yes;
colour yes;
animationInfo
{
colour field;
colourMap rainbow;
colourField d;
// Optional colour map min and max limits
//min 0;
//max 0.002;
//alpha uniform;
//alphaValue 1;
alpha field;
alphaField d;
normalise yes;
}
}
```
Note that the particle colour for animated tracks remains constant; dynamic colour values are not currently supported by the format.
## Examples
See:
- $FOAM_TUTORIALS/lagrangian/reactingParcelFoam/filter/system/sample
- $FOAM_TUTORIALS/lagrangian/reactingParcelFoam/filter/constant/particleTrackProperties
- $FOAM_TUTORIALS/lagrangian/reactingParcelFoam/filter/constant/particleTrackProperties.animatev2112Mark OLESENMark OLESENhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/504Added functionality for smoothing the sensitivity derivatives2021-12-09T09:18:37ZAndrew HeatherAdded functionality for smoothing the sensitivity derivatives### Summary
Adjoint-based sensitivity maps can now be optionally smoothed using a Laplace-Beltrami operator.
When computing sensitivity maps on surface meshes generated from industrial
geometries, the outcome might appear noisy, especi...### Summary
Adjoint-based sensitivity maps can now be optionally smoothed using a Laplace-Beltrami operator.
When computing sensitivity maps on surface meshes generated from industrial
geometries, the outcome might appear noisy, especially if a volume-to-surface
approach is used for meshing, like the one utilised by snappyHexMesh. Even
though the sensitivity map is technically correct, the noisy patterns that
appear might make the extraction of useful information challenging.
This new feature facilitates the interpretation of the sensitivity map in such cases.
### Details of new models
The sensitivity map can be smoothed using a Laplace-Beltrami operator of the form
```math
\tilde{m} - \frac{\partial}{\partial x_j}\left[R^2\frac{\partial \tilde{m}}{\partial x_j}\right] = m
```
where $`\tilde{m}`$ is the smoothed sensitivity map, $`m`$ is the original
sensitivity map and $`R`$ a smoothing radius. The latter is computed as a
multiple of the average 'length' of the boundary faces, if not provided by the
user explicitly.
The above-mentioned equation is solved on the part of the surface mesh defined
by the patches on which the sensitivity map is computed, using the finiteArea
infrastructure of OpenFOAM.
If an faMesh is provided, it will be used; otherwise it will be created on the
fly based on either an faMeshDefinition dictionary in system or one constructed
internally based on the sensitivity patches.
From an optimisation point of view, this smoothing can alternatively be seen as computing the sensitivity derivatives $`\frac{\delta J}{\delta b_i}`$ of the objective function $`J`$ w.r.t. a different set of design variables $`b_i, i\in[1,N]`$, defined as
```math
x_i = x_i^{init} + \tilde{b_i} \\
\tilde{b_i} - \frac{\partial}{\partial x_j}\left[R^2\frac{\partial \tilde{b_i}}{\partial x_j}\right] = b_i
```
where $`x_i`$ are the coordinates of the updated geometry and $`\tilde{b_i}`$ a smooth displacement field. In other words,
no loss of accuracy is incurred by the smoothing; instead, sensitivities are computed w.r.t. a different set of design variables.
### Examples
An example of a noisy sensitivity map and its smoothed variants using different $`R`$ values is presented below, taken from the updated motorbike tutorial under
$FOAM_TUTORIALS/incompressible/adjointOptimisationFoam/sensitivityMaps/motorBike,
![portFront](/uploads/c7a772e62d98ea0a34414c7cd05ffa1c/portFront.png)
![port](/uploads/ec1bb55d3255e3dd403818ec21c09245/port.png)
![top](/uploads/9ddc8a60520cedba99780974f71e422f/top.png)
### Risks
No risks are foreseen. To enable the smoothing, the `smoothSensitivities` bool should be set to `true`, along with some optional entries.
```
sensitivities
{
type surfacePoints;
patches (motorBikeGroup);
includeSurfaceArea false;
smoothSensitivities true;
meanRadiusMultiplier 10; // Optional, defaults to 10. Controls the smoothing radius
iters 2000; // Optional, defaults to 500
}
```
Since the Laplace-Beltrami equation is solved using the finiteArea framework, `faSchemes` and `faSolution` should be present under `system`.
The smoothed sensitivity map is written in the `smoothedSurfaceSens + adjointSolverName + sensitivityFormulation` file, under the current time-step.
### References
The notion of sensitivity smoothing using a Laplace-Beltrami operator was used throught the seminal works of Antony Jameson. A reference can be found in
Vassberg J. C., Jameson A. (2006). Aerodynamic Shape Optimization Part I: Theoretical Background. VKI Lecture Series, Introduction to Optimization and Multidisciplinary Design, Brussels, Belgium, 8 March, 2006.v2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/505ENH: Adding dynamic-mesh motion capabilities to various solvers2021-12-16T09:05:01ZKutalmış BerçinENH: Adding dynamic-mesh motion capabilities to various solvers### Summary
With the commit, the `dynamicMeshMotion` functionalities and the following solvers can operate in combination.
- `buoyantBoussinesqPimpleFoam`
- `buoyantPimpleFoam`
- `rhoCentralFoam`
- `rhoCentralDyMFoam` has been merged...### Summary
With the commit, the `dynamicMeshMotion` functionalities and the following solvers can operate in combination.
- `buoyantBoussinesqPimpleFoam`
- `buoyantPimpleFoam`
- `rhoCentralFoam`
- `rhoCentralDyMFoam` has been merged with `rhoCentralFoam`
### Resolved bugs
N/A
### Risks
- No changes in user input.
- No changes in any output of the existing tutorials except those for `rhoCentralFoam`.
### Tests
* Compilation (incl. submodules):
* [x] `linux64ClangDPInt32Opt` (clang11)
* [x] `linux64GccDPInt32Opt`
* [x] `linux64GccSPDPInt64Debug`
* [x] Alltest: No change in output with respect to the develop HEAD + no errorv2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/508ENH: New liquid film features2021-12-09T17:02:19ZSergio FerrarisENH: New liquid film featuresAdding addition wall-function shear stress model to the film friction at the gas interface.
Updating tutorials accordinglyAdding addition wall-function shear stress model to the film friction at the gas interface.
Updating tutorials accordinglyv2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/509Additional sub-cooling heat transfer correlations for liquid H22021-12-15T10:16:50ZSergio FerrarisAdditional sub-cooling heat transfer correlations for liquid H21) Adding basic thermodynamic functions for rho and Cp at (T,p)
2) Adding heat transfer correlations for boiling film, transient and sub-cooling boiling regimes for liquid-H2.
3) Option to use a correlation-type of HTC for sub-cooling re...1) Adding basic thermodynamic functions for rho and Cp at (T,p)
2) Adding heat transfer correlations for boiling film, transient and sub-cooling boiling regimes for liquid-H2.
3) Option to use a correlation-type of HTC for sub-cooling regime instead of the standard RTI model.
4) Minor fix to the thermo for h-Polynomial type.v2112Mattijs Janssens4-Mattijs@users.noreply.develop.openfoam.comMattijs Janssens4-Mattijs@users.noreply.develop.openfoam.comhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/510ENH: interfaceOxideRate: new mass transfer model for icoReactingMultiphaseInterFoam solver2021-12-09T17:05:24ZSergio FerrarisENH: interfaceOxideRate: new mass transfer model for icoReactingMultiphaseInterFoam solver- Adding new mass transfer model based on:
Cao, L., Sun, F., Chen, T., Tang, Y., & Liao, D. (2018).
Quantitative prediction of oxide inclusion defects inside
the casting and on the walls during cast-f...- Adding new mass transfer model based on:
Cao, L., Sun, F., Chen, T., Tang, Y., & Liao, D. (2018).
Quantitative prediction of oxide inclusion defects inside
the casting and on the walls during cast-filling processes.
International Journal of Heat and Mass Transfer, 119, 614-623.
DOI:10.1016/j.ijheatmasstransfer.2017.11.127
The model calculates the oxide produced at the air-liquid interface of a VOF system.
- A new BC was created to quantify the amount of oxide mass absorbed by walls.v2112Andrew HeatherAndrew Heatherhttps://develop.openfoam.com/Development/openfoam/-/merge_requests/511Update of interIsoFoam and isoAdvection extending to work with porous media2021-12-13T17:11:35ZJohan RoenbyUpdate of interIsoFoam and isoAdvection extending to work with porous mediaExtension of isoAdvection and interIsoFoam to deal with a porous medium
occupying part of the cell volumes and described by a porosity volScalarField
taking values larger than 0 (cell full of solid) and less than or equal to one
(cell f...Extension of isoAdvection and interIsoFoam to deal with a porous medium
occupying part of the cell volumes and described by a porosity volScalarField
taking values larger than 0 (cell full of solid) and less than or equal to one
(cell full of fluid).
The current implementation deals with both the advection part of the problem
and the modification of the UEqn (p_rghEqn needs no changes).
The aim of the implementation is to make it invisible to the user who does not
work with porosity both in terms of efficiency, memory usage and code complexity.
It would be preferable if the required interIsoFaom changes could be done without
touching the solver files, e.g. implementing the UEqn modifications as an fvOption.
This seems to be impossible because we need to modify the existing terms in the
UEqn and the fvMatrix is build from right to left meaning the fvOption cannot
touch e.g. the fvm::div(rhoPhi,U)-term. The chosen solution is therefore to add an
#include "UEqnAddPorosity.H" line to the UEqn.H file AFTER the UEqn is constructed.
The activation of porosity treatment is controlled by a bool parameter called
porosityEnabled in a constant/porosityProperties dictionary. This is disabled by default.
If it is enabled the user must provide a porosity volScalarField in the 0 time dir.
The isoAdvection class has two new data members: porosityPtr_ to the porosity field,
is nullptr by default, and a bool called porosityEnabled_ (default false) read
from the constant/porosityProperties by the isoAdvection constructor.
It has been necessary to introduce "if (porosityEnabled_)" lines various places in
the code but most places this is not inside loops over cells or faces so it is
cheap. An exception is in the isoAdvection bounding functions but these are only
called for very few cells so this should not be very expensive.
At this point I am not sure that the interplay with the source terms Su and Sp.
These should probably also be divided by porosity when alpha is updated.
The interIsoFoam solver now reads the porosityEnabled bool in constant/
porosityProperties,reads the porosity field if needed and the porosity parameters
for the Darcy-Forchheimer force and added mass force from the constant/
porosityProperties dict.
I have introduced a porousAlphaCourantNo.H and porousCourantNo.H file as
replacements for alphaCourantNo.H and CourantNo.H to correct the adaptive time
stepping so it accounts for the increased flow velocity in porous cells if
porosity is enabled.
Also the "Phase-1 volume fraction" write out has been modified to take porosity
properly into account.
There are two test cases included:
1. discInConstantPorousFlow which is the same as discInConstantFlow but with a
porous zone in the middle of the domain where the disc elongates and moves faster.
2. A porousDamBreak testing porosity implementation and comparing with exp. data.
This code contribution is joined work between:
- Konstantinos Missios
- Niels Gjøl Jacobsen
- Henning Scheufler
- Johan Roenby
Johan Roenby, December 2021v2112Andrew HeatherAndrew Heather