[JayDeeUU]

I'd be grateful for any suggestions in relation to the following OpenFoam task, which is not covered by any of the tutorials relating to chtMultiRegionFoam. I expect what I'm attempting has been done by others, but I haven't chanced to find any forum entry or external tutorial that wholly covers it. Please note that the setup here is for 2D-axisymmetric calculations.

Attached Figures (1,2) show (in paraFoam) two (meshed) fluid regions, that had been generated using Gmsh. One, a sphere (radius 6 m) with an access neck, represents the Snoplus neutrino detector ("detector"), and is filled with a scintillator fluid (Gmsh physical volume name and OpenFoam region name "scintillator"). The other region represents a volume of water in the "cavity" surrounding the detector (Gmsh physical volume name and OpenFoam region name "water"). The underlying interest is to compute buoyancy-driven circulation of the scintillator fluid, whereas some earlier simulations have modelled only the detector volume. Coupling to the water volume will allow to use a measured profile of temperature to drive the scintillator circulation.

The domain shown by the figures was created by rotating/extruding two plane surfaces, each of which was defined by a loop in Gmsh: the loop that (when extruded) defines the outer edge of the detector forms (also) part of the loop defining an edge whose extrusion defines the cavity (this is evident from the figures). In reality detector and cavity (i.e. the eventual OpenFoam fluid regions "scintillator" and "water") are separated by a 5.5 cm thick plexiglass wall [note: I originally mis-stated this dimension as 2 cm, consistent with the four layers of 0.5 cm seen in files below], but I'm hoping to be able to represent this (conducting) barrier without actually meshing it - somewhat along the lines of a 2011 forum post https://www.cfd-online.com/Forums/openfoam-solving/80114-two-sided-wall-heat-transfer-bc-no-separate-solid-mesh-no-heat-transfer-coeff.html . Also relevant is the 2016 tutorial by Dr. Johann Turnow, https://www.foamacademy.com/wp-conte...out_public.pdf

For what follows, it will be useful to state the names (assigned in Gmsh) for various "physical surfaces", names which hopefully are self-explanatory. There are four physical surfaces whose corresponding OpenFoam patches are of the "wedge" type used in 2D-axisymmetric calculations, namely westface, eastface, westface_cavity and eastface_cavity. The outer boundary of the detector (scintillator) is defined by three faces: lid, sidewall_neck, sidewall_shell. And finally the boundary of the cavity (water) volme is defined partially by these latter three, and completed by surfaces lid_cavity, sidewall_cavity and base_cavity.

Having copied the mesh file (domain.msh) to a case folder, and having set up regionProperties

regions
(
    fluid       (water scintillator)
    solid       ()                     
);

with two fluid regions but no solid regions, I issued

gmshToFoam domain.msh
splitMeshRegions -cellZones -overwrite

The result of these commands was that in constant/polyMesh/boundary there are seven patches:

7
(
    westface
    {
        type            wedge;
        physicalType    wedge;
        nFaces          10984;
        startFace       20578;
    }
    westface_cavity
    {
        type            wedge;
        physicalType    wedge;
        nFaces          2791;
        startFace       31562;
    }
    eastface
    {
        type            wedge;
        physicalType    wedge;
        nFaces          10984;
        startFace       34353;
    }
    base_cavity
    {
        type            wall;
        physicalType    wall;
        nFaces          4;
        startFace       45337;
    }
    sidewall_cavity
    {
        type            wall;
        physicalType    wall;
        nFaces          10;
        startFace       45341;
    }
    lid_cavity
    {
        type            wall;
        physicalType    wall;
        nFaces          4;
        startFace       45351;
    }
    eastface_cavity
    {
        type            wedge;
        physicalType    wedge;
        nFaces          2791;
        startFace       45355;
    }
)

which together bound the union of the scintillator and water regions. In constant/water/polyMesh/boundary there are six patches,

6
(
    westface_cavity
    {
        type            wedge;
        physicalType    wedge;
        nFaces          2791;
        startFace       4018;
    }
    base_cavity
    {
        type            wall;
        physicalType    wall;
        nFaces          4;
        startFace       6809;
    }
    sidewall_cavity
    {
        type            wall;
        physicalType    wall;
        nFaces          10;
        startFace       6813;
    }
    lid_cavity
    {
        type            wall;
        physicalType    wall;
        nFaces          4;
        startFace       6823;
    }
    eastface_cavity
    {
        type            wedge;
        physicalType    wedge;
        nFaces          2791;
        startFace       6827;
    }
    water_to_scintillator
    {
        type            mappedWall;
        inGroups        1(wall);
        nFaces          317;
        startFace       9618;
        sampleMode      nearestPatchFace;
        sampleRegion    scintillator;
        samplePatch     scintillator_to_water;
    }
)

while in constant/scintillator/polyMesh/boundary there are only three patches:

3
(
    westface
    {
        type            wedge;
        physicalType    wedge;
        nFaces          10984;
        startFace       16243;
    }
    eastface
    {
        type            wedge;
        physicalType    wedge;
        nFaces          10984;
        startFace       27227;
    }
    scintillator_to_water
    {
        type            mappedWall;
        inGroups        1(wall);
        nFaces          317;
        startFace       38211;
        sampleMode      nearestPatchFace;
        sampleRegion    water;
        samplePatch     water_to_scintillator;
    }
)

Here notice that three Gmsh physical surfaces, "lid, sidewall_neck" and "sidewall_shell", do not appear - having been compounded to make patches scintillator_to_water and water_to_scintillator.

I set up the boundaries in 0/scintillator/T as:

boundaryField
{
    westface
    {
        type            wedge;
    }
    eastface
    {
        type            wedge;
    }    
    scintillator_to_water
    {
        type compressible::turbulentTemperatureCoupledBaffleMixed;
        Tnbr T;
        //kappa fluidThermo;
        kappaMethod fluidThermo;
        kappaName none;
        value uniform 0;
        thicknessLayers (0.5e-2 0.5e-2 0.5e-2 0.5e-2); // four layers of .5 cm
        kappaLayers (0.19 0.19 0.19 0.19);
    }

For present purposes specifying the thermophysicalProperties and other elements of the simulation is probably not warranted. The outcome is that when I run chtMultiRegionFoam I get warnings to the effect that "Wedge patch 'eastface' is not planar" (etc.), then (after various messages):

Creating finite volume options from "system/fvOptions"

Region: water Courant Number mean: 0 max: 0
Region: scintillator Courant Number mean: 0 max: 0
--> FOAM Warning : 
    From bool Foam::functionObjects::vtkWrite::readSelection(const Foam::dictionary&)
    in file vtkWrite/vtkWriteUpdate.C at line 252
    No mesh regions selected for function object walls
Region: water Courant Number mean: 0 max: 0
Region: scintillator Courant Number mean: 0 max: 0
deltaT = 0.001200048
Time = 0.00120005

Solving for fluid region water
#0  Foam::error::printStack(Foam::Ostream&) at ??:?

My best guess is that I am not correctly handling the (unresolved, unmeshed) conductive wall (the patches scintillator_to_water and water_to_scintillator). I'd appreciate any suggestions. Thanks, John.

[dlahaye]

Did you forget the attachments?

[JayDeeUU]

Thanks Domenico, yes I did forget.

[dlahaye]

1/ What is 0/water/T (just to make sure)?

2/ Is chtMultiRegionFoam confused because geometry/mesh definition assumes zero layer thickness whereas scintillator_to_water defines defines a non-zero thickness? It it valuable to try a run in which zero thickness in scintillator_to_water is assumed (as it should not take much)?

3/ I would like to understand the error message "wedge not planar" better than I currently do. I tried to extend a similar case (Sandia Flame D tutorial) with a solid domain and got the "wedge not planar" error. Possibly one should simplify the geometry is a first attempt.

[JayDeeUU]

Thanks for your suggestions Domenico. I'll respond one-by-one:


1) 0/water/T was:

boundaryField
{
    westface_cavity
    {
        type            wedge;
    }
    eastface_cavity
    {
        type            wedge;
    }    
    water_to_scintillator
    {
        type compressible::turbulentTemperatureCoupledBaffleMixed;
        Tnbr T;
        //kappa fluidThermo;
        kappaMethod fluidThermo;
        kappaName none;
        value uniform 0;
        thicknessLayers (0.5e-2 0.5e-2 0.5e-2 0.5e-2); // four layers of .5 cm
        kappaLayers (0.19 0.19 0.19 0.19);
    }      
    base_cavity
    {
        type            fixedValue;
        value           uniform 290;
    }
    sidewall_cavity
    {
        type            fixedValue;
        value           uniform 290;
    }
    lid_cavity
    {
        type            fixedValue;
        value           uniform 290;
    }    
}

2) Following your suggestion, I amended the prescription for the (implicit) boundary patch (in both 0/scintillator/T and 0/water/T) to:

type compressible::turbulentTemperatureCoupledBaffleMixed;
Tnbr T;
//kappa fluidThermo;
kappaMethod fluidThermo;
kappaName none;
value uniform 0;
thicknessLayers (0.0); 
kappaLayers (1.e10);

When I ran the code again with this alteration again I got many warnings about non-planar wedge patches. Below I show only the final warning and the subsequent crash:

--> FOAM Warning : 
    From virtual void Foam::wedgePolyPatch::calcGeometry(Foam::PstreamBuffers&)
    in file meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C at line 72
    Wedge patch 'eastface' is not planar.
At local face at (4.6547833 0.20323223 3.5138368) the normal (-0.043619485 0.99904822 -1.811081e-09) differs from the average normal (-0.043619388 0.99904822 -2.125639e-11) by 9.4842994e-15
Either correct the patch or split it into planar parts
*** Reading fluid mesh thermophysical properties for region water

    Adding to thermoFluid

Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState Boussinesq;
    specie          specie;
    energy          sensibleInternalEnergy;
}

#0  Foam::error::printStack(Foam::Ostream&) at ??:?

3) Yes maybe I shall have to do some experimenting with a simpler domain. But right now I'm in the process of going in the opposite direction and making it more complicated... by resolving the plexiglass shell (5.5 cm thick - I mis-stated the thickness earlier) in Gmsh. I'll add to this post when I have some progress, positive or negative.


Thanks again!

[dlahaye]

Thx for your input.

The runtime error in your post #5 seems to be different than in your post #1. Is this correct? Is post #5 any progress over post #1? Is the solver in post #5 getting beyond a first solve on one of the solid domains?

[JayDeeUU]

Yes, you're correct about there being a difference in outcome. With the original 4 layers the solver does report that it is "Solving for fluid region water", but not so in the altered case, that is, the solver in post #5 is not getting beyond a first solve on any of the regions.

I've re-run both setups, and I'll attach the log files. I inserted into the bottom of the log file the text that appeared in the terminal... and had to cut out a vast volume of the warnings (to reduce file size).

[JayDeeUU]

I've made some progress on this. I created the (wedge) mesh anew in Gmsh, this time resolving the plexiglass wall that separates scintillator from water. Then, starting with a new case folder that was empty but for

 /constant/regionProperties
/system/controlDict
/system/fvSolution
 /system/fvSchemes

and the .msh file, I ran "gmshToFoam" and then "splitMeshRegions -cellZones -overwrite" . I kept track of files that were newly created by these commands, and the naming of patches (some flowing from choices of names in Gmsh, some automatically generated). There were no problems.



I changed the patch names as necessary (e.g. for the bounding side planes of the wedge, type "wedge" is obligatory), copied back in the thermophysicalProperties files and so forth, and amended the 0/region/U,T etc. files. The program still crashed, with messages reporting the arising of absurd temperatures. Proceeding by trial and error I altered the "equationOfState" choice for the fluids from "Boussinesq" to "rhoConst" - the program runs only provided neither fluid is treated using the Boussinesq option. I've attached an image of the UZ field at t=100 s. Here the T field was initialized as being isothermal, so there is no reason for motion to develop, and indeed practically it does not. In the case shown I set nNonOrthogonalCorrectors to 5 (whereas if I set it to 1, the motion field is about 5 orders of magnitude more vigorous). One other point to note is that with this new mesh, resolving the plexiglass layer, the solver is not giving the warnings that had occurred with the original mesh ("Wedge patch 'eastface' is not planar, etc).

Going ahead the rhoConst option will not be suitable, so I will have to plug on and try to get the "Boussinesq" option to work, at least for the scintillator (it has been fine during months of work with other solvers). I'll pursue a more suitable formulation for the water (whose coefft. of thermal expansion of course varies strongly with temperature). In past work I've coded in a spatially-uniform "background" thermal stratification by including a term in the TEqn.H (where w is the vertical velocity UZ and is a positive constant, e.g. 0.2 [K/m]), and interpreting the T variable as a perturbation (acting on the motion through the buoyancy term). Then I can initialize with T isothermal, and force the flow by imposing a cool zone (in T) on the boundary. I'm hoping I may be able to do this here, too.

Finally, I've also revisited the case that initialized this post and discussion, which differs in that the solid zone was not meshed. I confess that I discovered the "fvSchemes" and "fvSolution" files for both regions, that I'd used in the attempts that triggered this entire post, were merely template files. With that oversight corrected, and with the "equationOfState" option "rhoConst" for both fluids, the solver now does run - albeit still with the warnings about non-planar faces. If I choose equationOfState as "Boussinesq" for either fluid, the solver crashes with huge negative temperature.

[dlahaye]

This seems like wonderful progress!

Possibly you can try to bound rho (especially in early stage of convergence) to get Boussinesq to converge. Alternatively, you can try to provide an initial guess generated using rhoConst to Boussinesq.

I would like to gain an understand of what caused the non-planar faces warning.

[JayDeeUU]

Thanks for the suggestions, which I'm going to try (to get Boussinesq to converge). I'll post the outcome. Regarding the cause of the non-planar faces warnings that arise when one runs the solver on the mesh (case "snoplus25") that does not mesh the plexiglass, here is an instance:

Create time

Create fluid mesh for region water for time = 0

--> FOAM Warning : 
    From virtual void Foam::wedgePolyPatch::calcGeometry(Foam::PstreamBuffers&)
    in file meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C at line 72
    Wedge patch 'westface_cavity' is not planar.
At local face at (6.1629055 -0.26907827 0.47202264) the normal (-0.043619405 -0.99904822 -4.4252931e-08) differs from the average normal (-0.043619387 -0.99904822 2.0481131e-10) by 2.3293888e-15
Either correct the patch or split it into planar parts
.
.
.
Next line was numbered 48516... i.e. screeds of these warnings

--> FOAM Warning : 
    From virtual void Foam::wedgePolyPatch::calcGeometry(Foam::PstreamBuffers&)
    in file meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C at line 72
    Wedge patch 'eastface' is not planar.
At local face at (4.6547833 0.20323223 3.5138368) the normal (-0.043619485 0.99904822 -1.811081e-09) differs from the average normal (-0.043619388 0.99904822 -2.125639e-11) by 9.4842994e-15
Either correct the patch or split it into planar parts
*** Reading fluid mesh thermophysical properties for region water

    Adding to thermoFluid
.
.

Recall this is happening with a "wedge" type domain (azimuthal span ) for 2D-axisymmetric simulations, and it was produced in Gmsh. There's no way (so far as I can see) to inspect the mesh in Blender or MeshLab, because the mesh output format is .msh. The gmshToFoam command gives what is needed by OpenFoam, and the OF command "checkMesh -allTopology -allGeometry" indicates that this mesh fails two checks (see attached file "checkMesh_snoplus25.txt"). One can look at the mesh using paraFoam. What is odd is that in some views paraFoam displays spurious blotches (for want of a better word) - see attached image.

Comparing the above diagnostics with those pertaining to a similar mesh (case "snoplus25a") that DOES resolve (mesh) the plexiglass boundary and that DOES NOT generate the "non-planar'' warnings when the solver is run, checkMesh still flags one issue (see atttached file checkMesh_snoplus25a.txt), but, no matter how this domain is viewed in paraFoam, there are none of the oddities seen in case snoplus25.


checkMesh_snoplus25.txt

Mesh_snoplus25_inParaFoam_internalMesh_PLUS_waterInternalMesh_scintInternalMesh_vtkBlockColors_n.png

checkMesh_snoplus25a.txt

Mesh_snoplus25a_inParaFoam_internalMesh_plusAll_plexiToWaterPatches_vtkBlockColors.png

[JayDeeUU]

I have attached here the .geo file (renamed with extension .txt) corresponding to the bad mesh of case snoplus25 (plexiglass not meshed). It can be used to reproduce Gmsh's .msh file which can then be used with gmshToFoam.

[dlahaye]

Interesting.

I checked that the "not planar" error is given by the function Foam::symmetryPlanePolyPatch::calcGeometry with implementation here https://github.com/OpenFOAM/OpenFOAM...anePolyPatch.C

We should be able to loop over all the cells constituting the wedge patch, print the cell normal and see by how much they differ.

[JayDeeUU]

In case this should help anyone encountering the difficulty described above (chtMultiRegionFoam crashing when fluid density is treated using the Boussinesq approximation), I have eventually stumbled on the recommendation (e.g. at http://www.wolfdynamics.com/training...F8/all-2x2.pdf ) to set "dpdt off" in the fluid thermophysicalProperties files - where, by the way, I have chosen:

thermo hConst; 
equationOfState Boussinesq; 
energy sensibleEnthalpy;

This worked for me. It's the spirit of the Boussinesq approximation that absolute pressure is unimportant, only its gradients matter, and, density is a unique function of temperature.

Although I continue to have problems getting a wedge-type mesh (for 2D axisymmetric calculations) that passes "checkMesh -allTopology -allGeometry", one of many attempts does pass, and is adequate to at least progress with the project.

[dkeogh]

Hi J.D.,

Did you figure this problem out? I am currently running computations on an axisymmetric mesh which I then split and use in chtmultiregionfoam. When I use the mesh in the axisymmetric computations before splitting by region everything works fine. Once I split the mesh though I get the same "not planar" error message. It is possible for me to use a three-dimensional mesh but it would be much more computationally demanding...

[JayDeeUU]

Hi, I don't have any new insights on that problem of the "not planar" messages, and my reply here is liable to be a bit unfocused, sorry. At the moment I'm not using chtMultiRegionFoam (I'm just focused on just one, central region), but I'll probably go back to it eventually.

One recent discovery (for me) is that after using "gmsgToFoam" one needs to set up the constant/polyMesh/boundary file before running checkMesh - on that score, in the past I must on some (but probably not all) occasions have wrongly concluded I had a bad mesh.

I think I also experienced that same (or a similar) message about surfaces not being planar when I used the "symmetryPlane" boundary condition on the bounding walls of a single-quadrant (rather than 5 degree wedge) mesh (also generated using Gmsh). I switched to using the "symmetry" condition instead, after reading the following forum discussion: symmetry and symmetryPlane boundary condition difference in OpenFOAM

[dkeogh]

Hi J.D.,

I managed to solve the problem although I have been using fluent to generate my mesh. I find that the problem seems to go away when I execute splitMeshRegions -cellZones -overwrite with a write precision of 20. I found that I had to change the write precision in controlDict before splitting the mesh. I did not have to change the writeFormat to binary as has been suggested in other posts. Also, I was able to keep the wedge boundary condition.

[JayDeeUU]

Thanks Declan. Glad you've found a way around the problem.