Liquefied gas simulation (heat driven flows)
 

Hi Everyone,

I would like to consult with the community the problem I am facing. I am new to the CFD so probably majority of my questions/problems is due to my short experience, so please excuse me.

I am intereseted in convection flows of a liquid gas, i.e. I would like to compute the U field for liquid argon in a cryostat. The cryostat is not perfect, it has some inward heat flux, uniformly distributed over its surface. To compensate the flux there is a cooling system, which basically compensates the heat (the same flux for equilibrium). The cooling system is placed at the top of the cryostat:

-----
| |
#~~# <- cooling surface and liquid level
|~~|
|~~|
\---/

I do not mind so much about the volume of gas above. I expect that including evapouration effects would be difficult...

Until now I've prepared the model in GMSH, and exported it to OF. The model was prepared in mm dimensions (3D), so I scaled it using known utility down/up to metres (this was a problem at the beginning - huge divergence of the solution).
The model consists of physical surfaces like cryostat walls and a cooling surface, and a volume of liquid gas (and a few construction details inside the cryostat described later). The walls are a source of positive heat flux (inwards), the cooling surface has a negative heat flux.
I've set up the solution using buoyantBoussinesqSimpleFoam solver (based on iglooWithFridges), and here come the questions:
1) I chose the solver as my problem is incompressible one. Is it a correct choice?
2) I set up fields: alphat (uniform 0), epsilon (uniform 2000 to avoid convergence errors right at the beginning), k (uniform 0.1), kappat (uniform 0, Prt 0.85), nut (uniform 0), p (uniform 0), p_rgh (uniform 1E5), T (uniform 87.3) and U (uniform 0). I have doubts about "p" field and I do not know exactly where should I put information about rho of the fluid (is that nut field?). The boundary conditions for T are the following: for the heating surface (almost all cryostat walls): type is fixedGradient, and a value of 212, for the cooler value of -2800 (its surface is smaller). The numbers come from the computation using Fourier's law: the heat flux is around 25 W/m2 for the cryostat, and -330 W/m2 for the cooler, thermal conductivity is 0.12 W/mK. If I run the solver, it diverges quite fast, i.e. the U field is silly after a 100 steps (1E19 m/s etc.), and the temperature of the liquid is uncomparable with the real device (I have a cryostat and I measure the temperatures inside very precisely).
3) Transport properties are calculated by myself, and therefore there could be a mistake: nu 0.2695, beta 3E-3, Pr 2.234, Prt (guessed) 0.8. Did I miss something?
4) The mesh is quite accurate, down to a few mm cells, but there are some physical surfaces inside the cryostat denoting metal elements immersed in the liquid volume. For those I set the constant BC, i.e. temperature of 87.3 K as I do not know how to tell the software that it should compute the temperature from the surrounding liquid (I assume metal is conducting enough to equilibrate with the liquid fast). Is it OK, and how to tell the OF to compute the temperature for the "inner" surfaces?
5) Regarding convergence, I've already chaged relaxation factors to improve it a little bit. Time step is 1 s; recently I used 0.01 s but with no success (I tried to adjust the problem time-scale). Should I change anything here?

I think the problem is quite easy to be solved with OF without an extra effort in implementing own solvers etc. The problem is that I've never reached the steady state conditions, the U field is just unbelievable. If I lower the heat flux, let's say 1000 times then the problem doesn't diverge, but this is not a physical situation.
Any comments would be much appreciated!

Regards,
Krzysztof

Are you certain the heat flux balances out perfectly? Otherwise it will never reach a steady state.

I am not, and I am aware of that. Some of the energy is accumulated in the kinetics, also the integral heat flux is dependant on the real area of the surfaces, computed form the mesh. As I do not have the exact values of the areas I just approximate the fluxes (in reality I have this gaseous "buffer volume" which by evapouration/liquefication compensates the heat flux perfectly - also by changing the pressure above the liquid volume).
I thought of running the simulation up to the point in which the temperature distribution reflects what I have in reality. I would consider that point as a "steady state". That's why I chose so small time step at the beginning. (From my experience in reality the liquid in cryostat equilibrates within 12 - 24 h after a considerable change in cooling/heating flux).

buoyantBoussinesqSimpleFoam is a steady state solver, there is no time derivative term in the governing equations. Thus your notion of time doesn't really match up here. As an intermediate step, you could test using a fixedValue temperature instead of a cooling heat flux, and check if you get something that makes sense.

- Anton

Concerning the time - this I know, just the time-step seems to me to be relevant. I'll give a try to fixed temperature BC's. If I know the temperature profile along the z-axis in the cryostat, how could it be implemented? Simply saying in T BC definition file "calculated" and giving something like "87.3 + (z - 0.5)*0.01" would work? (In reality I have only point readouts of the temperature, but the functional approximation is good enough).

Regards,
Krzysztof

OK, so I've implemented BC with fixed temperature, and there are still problems:

I am using mesh created in gmsh, then during the conversion step (gmshToFoam) I get a warning:

--> FOAM Warning :
From function polyMesh::polyMesh(... construct from shapes...)
in file meshes/polyMesh/polyMeshFromShapeMesh.C at line 619
Found 38306 undefined faces in mesh; adding to default patch.

Assuimg this is not fatal, I run the simulation and the U field seems to be computed OK, but in one certain point of the volume the U is really huge (for sure it is a computational error, maybe the mesh is not OK?). This point is located next to the wall. I inspect this by running Stream Lines plugin in paraFoam. The stream lines end up and begin in the same point (black and white hole at the sime time?). But I set up U field to fixed (0, 0, 0) value for all the BC, so I have really no idea of the source of the problem.

Have you ever encountered such behaviour?

Best wishes,
Krzysztof

Based on the issue with what seems to be some kind of convergence of U, I think the mesh, while still computable, might be giving the solver a problem. Try running
on the case and see if it tells you that the mesh is alright. If it is, then I don't know what else to suggest for the moment. Is there a blockMeshDict file you could post? Or does it not create one when you convert from GMSH to OF?

Cheers,
~Ed

Hi,

I run the checkMesh utility and the result is OK, except:

But this still doesn't look severe. I don't have dict file as the mesh is created by gmsh utility. So I only do gmshToFoam conversion and scaling. I've also noticed that during mesh optimization in gmsh I have a few "bad quality" tets, but..? Might it be the source of the problem?

I've just found that I can do

and these non-orthogonal faces correspond to places where U diverges. I have some sharp corners in the geometry, for which I understand why this situation occures, but there are also smooth surfaces, on which in a few points I found such nonorthogonality.

Have you any idea how to improve the mesh after it's being created (this takes really some time...)?

Hmm. The only thing I could even think to suggest is that you check your Courant number for the smallest cells in your simulation (as reported by checkMesh). The finite volume method doesn't require faces to be orthogonal, so that shouldn't be it, I don't think.