[JayDeeUU]

I'd be interested to hear suggestions (and criticisms) in relation to the following impasse or paradox or spurious intention - whichever it may be. For several months I've been attempting to get a sense of the likely buoyancy-driven circulation patterns in the Snoplus neutrino detector: a closed container that has axial symmetry about the vertical axis, the spherical part of the vessel having radius 6 m and the "neck" having radius 0.75 m (see Figure of mesh). In view of the symmetry of the domain and forcing, in the latest efforts I've been treating the flow as 2D-axisymmetric, with the domain being the required 5 degree wedge. I'm not entirely happy with the mesh (as to resolution of boundary layers), but I don't believe that's central to the discussion here: checkMesh finds it OK.

The fluid in the vessel is a scintillator of known density, thermal expansion coefficient, specific heat, viscosity and thermal diffusivity. An important aspect of the flow regime is that (experimentally) the fluid is kept in a state of stable thermal stratification (about 0.2 K/m). Accordingly, any motion resulting from thermal inhomogeneity (or non-stationarity) of the wall temperature is expected to be weak; and therefore it seems essential to configure an initial state (for time-dependent simulations) that, absent any imposed forcing, would not evolve. This is one of the reasons for choosing to use the incompressible buoyantBoussinesqPimpleFoam ("bBPF") solver that is available with OFv2006, for here the density is treated as constant except in the body force term (gravity) - and indeed the fluid density is not even accessible to (because not needed by) the solver. I did make an excursion into the possibility of using the compressible solver buoyantPimpleFoam with the 'eqnOfState Boussinesq' option, but despite a prolonged effort to tediously manicure the initial profiles (of T, p, p_rgh) for hydrostatic self-consistency, I never achieved an initial (motionless) state that didn't break out into exuberant motion under the action of the solver.

Now, using bBPF (or, for that matter, bBSF - if one wanted steady-state simulations) I aim to resolve the influence of the stable "background" thermal stratification, and, a perturbation temperature field T amenable to being used to disturb an initially motionless state. Let's suppose there's a "background" component of the total temperature field whose vertical gradient is temporally and spatially uniform. To represent its influence, one may modify TEqn.h by adding (into the right hand side of the T-equation) an advection term , where is the vertical velocity ("UZ" in OpenFoam) and is a specified constant, accessed from constant/transportProperties (in that same file, I set TRef=0). In effect, if the added term alone were acting, then the temperature perturbation T for a parcel rising (descending) along the stable background gradient would tend to decrease (increase), attenuating its buoyancy.

The conundrum is this: even with the perturbation temperature T held fixed at T=0 on all wall patches, and with or without the background stratification, the OF solution develops a field of motion that, in terms of its magnitude, is not qualitatively different from what arises if (say) one invoked a mildly inhomogeneous wall perturbation temperature (see especially Figure 5). How is one to recognise or discriminate a weakly-forced buoyancy-driven motion, if an undriven solution develops comparably energetic motion? Readers experienced with OpenFoam may wonder how I could have expected life to be that simple. I have sometimes solved the discretized RANS equations for turbulent micrometeorological flows, in which normally one knows (and represents) the proper form of the inflow profiles (of mean alongstream velocity U=U(z), turbulent kinetic energy and so forth, z being the vertical coordinate) based on Monin-Obukhov theory: in simulations of "disturbed" flows a first criterion for validity of the code is that the inflow profiles should be preserved if the code is run without whatever feature it is that would (otherwise) cause the disturbance (e.g. a windbreak, or an undulation of terrain). At the upwind boundary, which lies in the region of undisturbed motion, the vertical velocity component W=0 (by requirement). The imposed inflow profiles should be equivalent to a 1D solution of the problem, and if they are, then the vertical velocity should remain negligible (apart from the influence of roundoff errors) everywhere in the flow domain. Furthermore it is simple to judge what could be termed a negligible disturbance to the flow - say, |W| everywhere smaller than max|U|/10^6.

So perhaps there's a paradox intrinsic to what I've been trying to do. I don't have a criterion for "weak circulation", nor for "weak disturbance". Setting that possibility aside, I also have a secondary conundrum. This latter is not fundamental, but relates, very likely, to my inexperience with OpenFoam - although others have alluded to similar difficulty (https://develop.openfoam.com/Develop...s/-/issues/946). I would like to "drive" this flow by imposing a cooling heat flux to the vessel walls, over (say) the top metre of the "neck" (I have a distinct patches for the sidewall of that section, and its "lid"). Unfortunately however, although "buoyantBoussinesqPimpleFoam -listScalarBCs" returns a list that includes "type externalWallHeatFluxTemperature" I haven't been able to make it work. As an alternative I've considered "type fixedGradient", but here the problem is that I don't know how to supply the (fluid side) thermal conductivity (laminar plus turbulent) needed in order to specify the gradient. To be able to impose a heat flux as boundary condition seems, a priori, a basic and essential option, so it is surprising to find that it is not automatically available.


Images, in order of left to right. (1) the mesh. (2) "undriven" solution at t=5 hrs with no background temperature gradient and T=0 on all wall patches. (3) "undriven" solution at t=5 hrs with background temperature gradient 0.2K/m and T=0 on all wall patches. (4) "driven" solution at t=5 hrs with background temperature gradient and T=-0.1 on walls of the uppermost 1 m of the vessel (T=0 on other walls). (5) Profiles (at t=5 hrs) of vertical velocity along a vertical line at radius 0.74 m, i.e. from base (z=-3.6 m) to top (z=12.75 m) of the vessel, 1 cm inside the neck.


Additional information added 11 March 2021 to clarify the above:

I used the term "solver noise" without much consideration. Perhaps a better term would be "false flow"? What I meant to identify by the term was this: whatever "solution" (i.e. velocity field and all other fields dynamically connected to it) emerges spontaneously when "the solver" (as a package, including all the discretionary choices, such as the time step, discretization schemes, tolerances, relaxation factors etc.) is applied to a given mesh with boundary and (static) initial conditions that should prevent the development of motion.

I should, therefore, have specified in more detail the various choices I was using: deltaT 0.1, adjustTimeStep yes, maxCo 0.5.

fvSchemes:

ddtSchemes
 {
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    
    div(phi,U)      Gauss upwind;
    div(phi,e)      Gauss upwind;
    div(phi,h)      Gauss upwind;
    
    div(phi,T)      Gauss upwind; 
    div(phi,C)      Gauss upwind;

    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;

    div(phiv,p)     Gauss upwind;
    div(phi,K)      Gauss upwind;
    
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
    
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fvSolution:

solvers
{
    "rho.*"
    {
        solver          diagonal;
    }

    p_rgh
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-8;
        relTol          0.01;
    }

    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

    "(U|e|k|epsilon|T|C)"
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0.1;
     }