[maddalena]

Hello everybody,
I am going to simulate a conduction + convection problem in a multi region case, thus including solid regions and fluid regions. I am using OF 1.6.x, June update.

As a first step, I modified the standard solveSolid.H to include heat generation. I tested this approach on multiRegionHeater tutorial: I set all the temperatures to 300 K and waiting for the heat generation having its effect. That worked.

Than I created a simple model to understand how to proper set up the solver for a case similar to what I want to investigate. In that case, I have a slice 1 cell thick, with a cyclic boundary condition on the sides (see the attachment). There are different regions, corresponding to different materials: A to E are solids, G is fluid (laminar). The solid mesh is hexa, while the fluid is made of prisms. checkMesh did not complain:

...
Checking geometry...
    Overall domain bounding box (-0.191013 -0.1 0) (0.1910131 6.735557e-17 1.1)
    Mesh (non-empty, non-wedge) directions (1 0 1)
    Mesh (non-empty) directions (1 0 1)
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (1.818783e-16 0 8.783043e-17) OK.
    Max cell openness = 2.201851e-16 OK.
    Max aspect ratio = 57.5877 OK.
    Minumum face area = 1.710101e-05. Maximum face area = 0.001553425.  Face area magnitudes OK.
    Min volume = 1.710101e-06. Max volume = 7.959095e-06.  Total volume = 0.02111826.  Cell volumes OK.
    Mesh non-orthogonality Max: 34.24404 average: 8.211823
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.5108016 OK.

I did not change fvSchemes and fvSolution at the moment, since I think my case is still simple to be solved with standard schemes. In any case, I modified controlDict to take into account the maxCo and the maxDi.
However, all my run were not successful: different maxCo and maxDi always bring to a sudden divergence on h or Uz equation on fluid region, while the solid regions did not get a sensible temperature increase. For example: at time 7.30839 the log file for F says:

Solving for fluid region F
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 2.55267e-05, Final residual = 3.744047e-09, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.0001146642, Final residual = 5.379014e-09, No Iterations 7
DILUPBiCG:  Solving for h, Initial residual = 0.0001534736, Final residual = 4.188709e-10, No Iterations 2
Min/max T:263.3289 300.0309
GAMG:  Solving for p, Initial residual = 0.1169048, Final residual = 1.614516e-06, No Iterations 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors (F): sum local = 9.09263e-11, global = 8.529874e-14, cumulative = 1.255866e-07
GAMG:  Solving for p, Initial residual = 0.08567001, Final residual = 4.9401e-09, No Iterations 5
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors (F): sum local = 2.772171e-13, global = -3.119677e-15, cumulative = 1.255866e-07

and the following step...

Solving for fluid region F
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 2.552565e-05, Final residual = 3.75006e-09, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.0001146813, Final residual = 2.374168e+27, No Iterations 1001
DILUPBiCG:  Solving for h, Initial residual = 0.00015347, Final residual = 4.189117e-10, No Iterations 2
Min/max T:263.3284 300.0309
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.005713555, No Iterations 4
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors (F): sum local = 4.336668e+08, global = -8.60693e+07, cumulative = -8.60693e+07
#0  Foam::error::printStack(Foam::Ostream&) in "/root/OpenFOAM/OpenFOAM-1.6.x/lib/linux64GccDPOpt/libOpenFOAM.so"
#1  Foam::sigFpe::sigFpeHandler(int) in "/root/OpenFOAM/OpenFOAM-1.6.x/lib/linux64GccDPOpt/libOpenFOAM.so"
#2  ?? in "/lib/libc.so.6"
#3  Foam::GAMGSolver::scalingFactor(Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double> const&, Foam::Field<double> const&) const in "/root/OpenFOAM/OpenFOAM-1.6.x/lib/linux64GccDPOpt/libOpenFOAM.so"
#4  Foam::GAMGSolver::scalingFactor(Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double>&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, Foam::Field<double> const&, unsigned char) const in "/root/OpenFOAM/OpenFOAM-1.6.x/lib/linux64GccDPOpt/libOpenFOAM.so"
#5  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const in "/root/OpenFOAM/OpenFOAM-1.6.x/lib/linux64GccDPOpt/libOpenFOAM.so"
#6  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/root/OpenFOAM/OpenFOAM-1.6.x/lib/linux64GccDPOpt/libOpenFOAM.so"
#7  Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/root/OpenFOAM/OpenFOAM-1.6.x/lib/linux64GccDPOpt/libfiniteVolume.so"
#8  main in "/root/OpenFOAM/root-1.6.x/applications/bin/linux64GccDPOpt/chtMultiRegionHeatSourceFoam"
#9  __libc_start_main in "/lib/libc.so.6"
#10  _start in "/root/OpenFOAM/root-1.6.x/applications/bin/linux64GccDPOpt/chtMultiRegionHeatSourceFoam"

Tracing back the error, my logs/Uz_iter shows something like that:

0                2
...               3
.....             3
7.0             4
....     
7.30649    4
7.30658    5
7.30667    5
7.30676    6
7.30685    7
7.30694    661
7.30703    33
7.30712    11
7.30721    12
7.3073    18
7.30739    7
7.30748    11
7.30757    18
7.30766    33
7.30775    8
7.30785    7
7.30794    7
7.30803    7
7.30812    7
7.30821    7
7.3083    7
7.30839    7
7.30848    1001

At the last time step, I can see that the highest Uz values are located close to the slice tip. In addition, since I set T as 300 K everywhere, the

Min/max T:263.3284 300.0309

does not seems logical... Something like that happens with different maxCo or maxDi, at different simulating time.

I am thinking that all my problems may be due to the model itself:

• is the chtMultiRegionFoam able to solve a cyclic + 2D model?
• what kind of mesh is usually suggested for a cht simulation? Tetra, hexa?
• is there any general rule for the mesh definition on the solid part? Something like 30<y+<100 for a turbulent highRe flow
• And what about the fluid region if the flow is laminar?
• How can I set U as uniform (0 0 0) at all the boundaries and in the internalField? I tried that, but a floatingPointError appear at the first time step, while it disappeards if I used U uniform (1e-5 0 0).

At the same time, it seems like there is no good references on cht simulations using cfd. All the cfd books I checked give some general rules, which applies mainly in cases not involving heat transfer...
Now, I am almost out of ideas... any suggestion is really welcome.

maddalena

[maddalena]

Hello,
I think it is a good idea to post some of the answer I have found up to now.

• Has the Min/max T:263.3284 300.0309 logical if the starting temperature is 300K?

No, it is not. The problem was due to the fact that I did not used a pressure waveTransmissive boundary condition on the top boundary... Actually there was a fluid region on the A solid region... How to set up a waveTransmissive BC? See here.

• is the chtMultiRegionFoam able to solve a cyclic + 2D model?

Yes, it is. However, the solution of this kind of problem is difficult if the tip angle is small and the solid surface is "close" to the tip, since the buoyancy of the flow has not enough space to develop. If a simulation like that is required, it is better to use a larger slice, or 1/4 of the domain with symmetry bc.

• what kind of mesh is usually suggested for a cht simulation? Tetra, hexa?

The usual rules applies: better to use a hexa layer on the fluid region at the boundary with the solid part

• How can I set U as uniform (0 0 0) at all the boundaries and in the internalField? I tried that, but a floatingPointError appear at the first time step, while it disappeards if I used U uniform (1e-5 0 0).

That was a bug, now fixed on 1.7.x. See here.

I have not find a good reference for cht simulations. Still waiting for your suggestions.

Cheers,

mad

[strakakl]

Hy Maddalena,

I am simulating a conduction-convection case, similar to your test setup in post 1. It's an axissymmetric case with a solid cylinder in the middle and some kind of heater band (solid) around the cylinder. Around the cylinder and the heater band is a fluid region.
I added a heat source to the solver for the solid regions. The heat source is set to zero in the inner cylinder and nonzero in the heater band.

The case is initialized with T=300K. After a about 10s the Temperature in the fluide region ranges from about Tmin=270K to Tmax=305K, so there is some mystic cooling .

I found in your post above that you have had the same problem some time ago and that you solved it using a pressure waveTransmissive boundary condition on the top boundary. I implemented this BC in my case but it didn't solve the problem at all.

What's about the other BC's in your case. Or do you have some other suggestions how to solve this problem?

Thanks for your help in advance.

Cheers, Klaus

[maddalena]

Originally Posted by strakakl

What's about the other BC's in your case. Or do you have some other suggestions how to solve this problem?

Hi Klaus,
due to the geometry, there are not other bc to fix, since everything else is cyclic.
Unfortunately, I completely changed approach (deleted the G region) and I cannot help you more than this...
Sorry.

mad

[strakakl]

Hy Maddalena,

no problem, thanks for your answer.

Klaus

[ulli]

Hi Maddalena, hi Klaus,

I look for some experts for cyclic CHT problems: I have a periodic geometry and would like to define cyclic patches between different regions. Please see attached picture: The right side of fluid region1 should be coupled with the left side of fluid region 2. As the CHT models the meshes complete separate for each region I wonder if such a cyclic case can be generated (and solved). So I think cyclic conditions can be possible with the same region but not across different regions (?).
Any suggestions?
Thanks in advance.

Ulrich

[maddalena]

Hallo,
do the two regions have the same properties? If so, than you can merge them (using mergeMesh) even if they do not touch each other, thus they become a single region. At this point, you can create the cyclics. hope to be clear...
btw, have you had this problem when exporting mesh from an external mesh generator?


mad

[ulli]

Thanks a lot Mads,
great idea. I will check this out. I just came to the problem when I tried to prepare my mesh generator for these kinds of problems. I looked up the options in createPatch to generate cyclic patches for identical mesh faces. But there was no region option... This of course wouldn't make sense, because a mesh face can only belong to a single region. So another kind of coupling would be needed, I think. But your idea sounds great and will most likely work for my case.

Best regards

Ulrich

[maddalena]

Originally Posted by ulli

Thanks a lot Mads,
great idea. [...]But your idea sounds great and will most likely work for my case.

to be honest, this is hrv's idea...
please, cold you report if you meet any problem when exporting from a mesh generator?
thank you!

mad

[ulli]

Hi Mad,

I can confirm that all is true:
I used createPatch to generate the cyclic patches for the complete domain first. Than I used splitMeshRegion to separate the fluid from the solid regions. This works pretty well. Thanks once more for your help!


Ulrich

[math.n3rd91]

Hello,

I am relatively new to OpenFOAM still, so I apologize if I end up asking some questions that seem too simple.

I am trying to set up a simulation in chtMultiRegionFoam that will consist of a solid cylinder containing a heat source (in W/m^3), surrounded by a gap that is filled with a gas, and then that is going to be surrounded by another solid, so it'll be: solid cylinder, gas annulus, and solid annulus. Outside that, I wanted to have a fluid flowing past. The purpose is to see how the fluid cools the materials from the volumetric heat source in the center.

I realize that you all have been talking about modifying code to incorporate the heat source, but seeing as how I'm still new to most of this, I don't really understand how to go about doing any of it.

Also, if anyone knows how to get the chtMultiRegionFoam to run with cylinders and annuli that would be great! (I think I set them up correctly, but they won't mesh very well and aren't all being counted as "regions" after I run it).

Thanks,

~ Denise

[lhcamilo]

Denise,

you would have to add a source term to the Temperature/Energy equation.

The first step here would be to derive the Temperature/Energy equation that you will be solving.

I would advise you to take a look at the works of I Tiselj, R Monod & G Brillant and Kasagi. they have dealt mainly with turbulent channels channels but I guess the principle is more of less the same.

I reckon that the first steps would be to decide what kind of simulation you will be running and boggle down on defining your governing equation in the context of finite volume methods.

your proposal is doable, but there are many variables to be considered.

In order to gain experience with tinkering in openfoam I would recommend you restrict yourself to simpler geometries with fewer regions. At most 2. Once you get the hang of that, then moving on to more complex geometry would be a lot simpler.


I wish you good luck and don't hesitate to ask for help.

[math.n3rd91]

Originally Posted by lhcamilo

Denise,

you would have to add a source term to the Temperature/Energy equation.

The first step here would be to derive the Temperature/Energy equation that you will be solving.

I would advise you to take a look at the works of I Tiselj, R Monod & G Brillant and Kasagi. they have dealt mainly with turbulent channels channels but I guess the principle is more of less the same.

I reckon that the first steps would be to decide what kind of simulation you will be running and boggle down on defining your governing equation in the context of finite volume methods.

your proposal is doable, but there are many variables to be considered.

In order to gain experience with tinkering in openfoam I would recommend you restrict yourself to simpler geometries with fewer regions. At most 2. Once you get the hang of that, then moving on to more complex geometry would be a lot simpler.


I wish you good luck and don't hesitate to ask for help.

Leonardo,

I appreciateyour response. I have an equation relating the temperature to the heat source (in W/cm^3) already, I guess my main problem is trying to incorporate that into OpenFOAM. I don't know if there is a certain program I need to put it into or what? If so, how? I know a bit about programming, but not enough that I would feel confident just guessing how to go about doing that.

Also, I have been tinkering around with the simpler geometries as well. The solver chtMultiRegionFoam has an example called "multiRegionHeater", but when I look at the blockMeshDict and then the makeCellSets.setSet, the size in the blockMesh is much smaller than the domain in the makeCellSets.setSet. I think this is where I'm going wrong in my meshing with the solids: I don't understand why the vertices were chosen, etc. especially since in the first tutorials, the blockMesh was over the whole domain.

Thanks,

~ Denise

[lhcamilo]

Originally Posted by math.n3rd91

Leonardo,

I appreciateyour response. I have an equation relating the temperature to the heat source (in W/cm^3) already, I guess my main problem is trying to incorporate that into OpenFOAM. I don't know if there is a certain program I need to put it into or what? If so, how? I know a bit about programming, but not enough that I would feel confident just guessing how to go about doing that.

Also, I have been tinkering around with the simpler geometries as well. The solver chtMultiRegionFoam has an example called "multiRegionHeater", but when I look at the blockMeshDict and then the makeCellSets.setSet, the size in the blockMesh is much smaller than the domain in the makeCellSets.setSet. I think this is where I'm going wrong in my meshing with the solids: I don't understand why the vertices were chosen, etc. especially since in the first tutorials, the blockMesh was over the whole domain.

Thanks,

~ Denise

well, I understand where you are coming from, but there is no way to make an omelet without cracking a few eggs. I started with openfoam doing channel flows using the default solver channelFoam. Then I modified modifying the pressure source term in the momentum equation so it would use a constant pressure source rather than a constant mass flow. Consequently I modified the solver so it would calculate turbulent energy budget terms. Once I was comfortable with all that I move on to multi-region solvers.

Comparing the differences between the momentum equations in pisoFoam and channelFoam should give you an idea about how source terms are add. In this case in channelFoam the source term is at the RHS of the momentum equation.

I don't know yet what solver you are going to use, but in case chtMultiRegionFoam there are two files containing the energy equations for the fluid and solid, namely fluid/HEqn.H and solid/HEqn.H . There is another version in the forums called chtIcoMultiRegionFoam that was made by Brazilian guy called Fabio Canesin and that version replaces HEqn.H by TEqn.H.

There are a bunch of slides and documents from the chalmers university on openfoam that you may find helpful.

As for your other question the parameters in makeCellSets.setSet sometimes go over the computational domain defined in the blockMeshDict. OpenFoam does not care if you go over the domain on makeCellSets.setSet, if is going to limit itself to the area inside the domain.

Say I have a domain composed by a 2D square defined by the diagonal {(0,0)(1,1) }, If I want to select just the lower half of the square I can use {(-1,-1)(1,0.5)} just as well as I could have used {(0,0)(1,0.5)}, openfoam is only going to care about the part of your selection that is inside the computational domain.

In any case good luck.

[math.n3rd91]

Thanks! I will give all of those a shot. I appreciate the help.