Adding heat source to chtMultiRegionFoam
Dear FOAMers,
in order to simulate a conduction + convection problem with openFoam, I added a heat source to the standard chtMultiRegionFoam release. This was straightforward: in solveSolid.H I added the heat source H [W m-3] contribution to the general conduction equation: Code:
{ After some more investigations, I noticed that the problem may be due to the coupling between interfaces. Indeed, the solidWallMixedTemperatureCoupled formulation calculates Twall as: Code:
scalarField Twall Code:
scalarField Twall I am thinking that I need a more general review of solidWallMixedTemperatureCoupled.C, but I am asking for suggestions and hints before going further:
Regards, mad |
coupling conditions
Following a presentation of the coupling made by HRV, I ended up with the following.
In the case of coupling between two regions A and B, both of them with a heat source H, the coupling condition should be something like:
I was thinking to implement this as: Code:
tmp<scalarField> myHDelta = 2*H()/patch().deltaCoeffs(); Is there anyone that can comment on that? Is this the right approach or am I on the wrong path? thank you mad. |
comments? hints? ideas?
Ok, the ideas presented above were unsuccessful, probably due to my poor C++ and / or something on the math formulation of the problem that I am missing.
Someone else agrees with me that the solidWallMixedTemperatureCoupled is not general enough to allow a heat source inside one or more regions. Unfortunately, it seems like that everybody facing this problem has given up and turned to a fixedValue approach. That is what I probably will do as well, at least for the moment. On the other hand, I am really interested on a new coupling formulation, and I'd like to spend some of my time on it. However, I need some general ideas on how to approach the problem, since I am still not completely into it. The first step, for example, can be some comments on the ideas described above... :rolleyes: regards, mad |
Hi,
first I would suggest you to check if the problem is really the boundary condition, did you try to compute the continuity error? For "engineering" porpoise the coupling mixed condition gives good result, but maybe the source term is giving troubles, did you try to solve implicitly, using SuSp()? If the power source is generated also along the boundary, which maybe is unphysical, you should add it in the balancing eq. otherwise the original should do the trick. Andrea |
1 Attachment(s)
Hi Andrea and thank you for replying!
Quote:
Quote:
Quote:
cheers, mad |
Hi,
about source term read this: http://www.cfd-online.com/Forums/ope...r-details.html If you don't want source at boundary create a vol Field with 0 fixed value boundary condition and add it to the equation to be solved, maybe this is worth a try. To compute errors you should write your own little code... not too difficult. Quote:
Hope it helps Andrea |
Hello,
Quote:
Code:
internalField uniform 500; I will try the SuSp() as soon as I understand where I should change the code.. cheers, mad |
Hi,
ok then, I'm curious, how much this heat source is large? Have you tried to lower H's value and see if there's any hope to get convergence? Quote:
Code:
SuSp(H) Code:
H.SuSp() |
Hi Andrea,
Quote:
Some thoughts about that... Let us for a moment drop the heat source contribution and consider the (working) case if I impose a fixedValue temperature. From what I can understand, the solidWallMixedTemperatureValue copies the temperature value from the neighbour region, thus from the "cold" region, on my cell in the "warm" region. On the "cold" region, a zeroGradient BC is applied. In different (and easier) words the temperature of the "cold" region is diffused on the "warm" region and viceversa. This works fine if I do not have any heat source. When I add H, the contribution of the "cold" region on the "warm" region is not enough, since the temperature gradient obtained on the boundary cell by coping "cold" cell value onto the "warm" cell value is lower than the heat source contribution. In other words, the temperature diffusion in the "cold" region works fine, but not in the "warm" region. Therefore, the "warm" region is somewhat isolated from the outside, and its temperature increases until the simulation ends. This is, at least, the explanation I am giving to the weird results I am getting. Quote:
Code:
tmp<fvScalarMatrix> TEqn In the meanwhile, another question come into my mind. chtMultiRegionFoam is a compressible solver, while I want to consider my simulation as incompressible. May I get problems if I use it in my case? I would say no, but maybe there is something that I am missing. Thanks for your time, regards, mad |
Hi again,
maybe you know about that but there was a sort of "bug" in the coupling patch in previous releases of OF, in 1.5 was fixed. Quote:
I still would like to see the coupling error between the two patches, can you try something as Code:
Info << "error: "<<Ts[i].boundaryField()[patchNumber]-Tf[i].boundaryField()[patchNumber]<< endl; you're welcome, Andrea |
Hi Andrea,
Quote:
Quote:
Quote:
Quote:
Quote:
Still working on the heat source definition. This is the error I am getting when using fvm::Su(H, T): Quote:
regards, mad |
Hi,
to check the coupling properly you should also check the heat flux, and for this you have to add two code lines :eek:, don't be too scared... If the heat is not flowing to the next regions of course the temperature will rise whatever the source implementation, I think you had a feeling about this too. I wait for some new updates, good luck! Andrea |
Quote:
Quote:
In any case, where should I add the lines you suggested? In solidWallTemperatureCoupled.H or in chtMultiRegionFoam.H? Is Ts the solid temperature at the boundary and Tf the fluid temperature? If so, what does change if I add the lines or check the T value at the coupling region? mad |
Ok then,
from what I see it seems fine except for a little heat flux error, maybe not enough for a totally wrong simulation, is this error so little in the first time steps? The problem is that you get unrealistic too high values inside the source domain, right? I noticed that the coupling gives slower temperature rising than it should, maybe you "accumulate" too much heat in the beginning (I'm just thinking aloud now) Andrea |
Hi andrea, and thank you to keep replying.
Quote:
Quote:
Quote:
1. the temperature is too high inside the source region. As a consequence, the temperature is too high inside the "cold" region as well. 2. The time vs temperature has not the theoretical tendency: in my simulation, I have a linear time vs temperature! Quote:
Quote:
1. from 0 to 25000 s with H = 5 W/m^3 2. from 25000 to 50000 s with H = 50 W/m^3 3. from 50000 to 100000 s with H = 500 W/m^3 no success. The only things that changes is the time vs temperature line angular coefficient. The theoretical curve rising is not matched. mad |
Hi,
so the coupling seems to work, even if with the "small" enough error in heat flow that unlikely it is the problem. So what about the others boundary conditions? It seems that you possibly have a zeroGradient somewhere in the colder domain. Let me know. Bye Andrea |
1 Attachment(s)
Hello Andrea,
find attached a simple figure of my domain. I have zeroGradient BC on K, cp and H (heat source) both at the interface between cube1 and cube2 and on the other vertical patch. As you can see, the grid is rough: should this be a problem? I guess that the solution may be wrong to the poor grid quality, but not SO wrong... In the meanwhile, I found a working case. If I add the heat source to the multiRegionHeater tutorial, everything works as it should: time vs temperature is not linear but it has the well known curved tendency. Note that the grid is really really rough on the tutorial. I do not now if I should be happy with that (OF does not betray us) or not happy (there is something wrong with my model). In any case... In order to find what is wrong with it, I cross checked what changes between the two models (tutorial and my own model):
cheers, mad |
Well,
I agree with your analysis, and then the obviously thing to do is drop the cyclic patch type... then we can re-think on "why you cyclic?" fingers crossed Andrea |
So...
cheers, mad |
ok,
don't loose your temper, the thing is weird but maybe we're still missing something. Can you post the case, along with the solution? I will take a look at it cheers Andrea |
This "WE" implies that there is someone else loosing sleeping hours for this as me??? ;)
Do you also need the solver? I will email them as soon as possible. thanks mad |
Hi mad,
I've looked at it. The problem maybe the symmetry/zeroGradient condition on the left side of cube1, that makes T unbounded, while a fixedValue works. My advise will be to try a cube2-cube1-cube 2 configuration. I opened this case and the tutorial, the cube case looks 10 times greater than it should be, but maybe is a paraview ocus-pocus. let me know Andrea |
Hello Andrea,
Quote:
Quote:
Quote:
Thanks for your support, at least now I have a possible answer to the problem! :rolleyes: ciao mad |
Code:
cube1 is a cube of 1 meter long. It seems to me that it is like that, no? Or did you compare it with the size of the tutorial? I have made them independently, so maybe I should scale my case to match the tutorial dimensions... Code:
Thanks for your support, at least now I have a possible answer to the problem! ok |
1 Attachment(s)
Hello,
Quote:
but... If I use a condition like the one shown in the attached figure... well, everything works! It seems like the solver needs a fixedValue somewhere in the domain to limit the heat source contribution... Indeed... coming back to the modified tutorial example. If you modify the minY BC from fixedValue 500K (or 300K) to zeroGradient BC, (and set U 0 0 0 in the fluid regions) well... The time vs temperature is ok in the first 100 simulated seconds, but it becomes linear after on! Thus it is not the "symmetry" or the "cyclic" or some special condition that does not work, but it is a "heat source limitation" problems. I am going to think to it. Does it has any physical meaning connected with the BC we are fixing or is it a computational problem? What is your opinion on that? mad |
Hi,
I agree, but maybe is more subtle than this. Of course source should be limited, but putting a fixedValue on the left side of cube1 and putting zeroGradient on the right of cube2 just works fine... hence, in my opinion, coupling condition misbehaves... Andrea |
1 Attachment(s)
Hi Andrea,
Quote:
One step more... I created a new case, like the cube2-cube1-cube2 case of yesterday and added a very large solid surrounding domain, where I fixed the temperature. Something similar is shown in the attached figure. In this case, we have NO BC on cube1, and its solution depends only on what is happening there. In addition, there is no involving convection, that may represent an additional problem. Now, we have two hypothesis on why the modified version of chtMultiRegionFoam fails:
After running as usual, I have result as usual: linear time vs temperature, hence 2 is true and, as you suggested yesterday, Quote:
Therefore, we are at the same point of the 20th of July, when I started this thread. How to improve the coupling condition? Or, is there a way to limit the heat source without fixing the temperature on one side of it? regards, mad |
Hi mad,
I asked the OF support team about the heat source issue. They told me I should take a look at the porousExplicitSourceReactingParcelFoam solver in OF-1.7.x There is a energy source term in createExplicitSources.H: Code:
Info<< "Creating energy source\n" << endl; Regards, Toni |
Ehi!
Quote:
Of course, I will have a look to it during the next week or so. And please report everything. The three of us can do a good job, I hope... :cool: mad |
Quote:
I looked at the file you pointed out. It seems that the createExplicitSources.H, which recalls timeActivateExplicitSource.H, is used only to define a heat source that can be activated / deactivated whenever you need it. In my case, the heat source is on and does not change for all the time, so this approach seems useless. Indeed the explicit source can be defined as I posted above. In any case, Andrea and I found out that is the coupling condition that does not work properly when defining a heat source and not fix at least one of the boundary values. If the temperature is fixed in one of the boundaries, the heat source defined as above gives reasonable results. The problem sounds deeper than it looks like. Hope I can still get some sort of support from the community, but at the moment none has suggested improvement for the coupling yet. Regards, mad |
Hi,
I will take a look at this also, can you send me the solver you have for now ??? You can send it to: fc [at] canesin [dot] com My problem is a little like yours, I have to add an heat source term that reads a variable cp and a variable value that both depends on another equation I solve in the domain.. (it is a triple coupled problem, magnetic+heat+fluid). I will be back in here for the next developments of my solver... hope to work together. |
Just send it.
Quote:
Cheers, maddalena |
Quote:
mad |
I have run a preprocessor in it and at the moment I'm trying to understand better the solver.
I agree it seams to do the interface like Patankar book. |
Quote:
tmp<fvScalarMatrix> TEqn ( fvm::ddt(rho*cp, T) - fvm::laplacian(K, T) - S ); Then added a field to it: PrtList<volScalarField> Ss(solidRegions.size()); Ss.set ( i, new volScalarField, IOobject ... ) ... In fact have copied temperature T field and just changed to "S". Then I have changed the flow solver to a incompressible one and removed the PIMPLE loop to a PISO one... now I'm trying to validate with Blasius and planewall.. but I'm having some problems with mesh =/ ... |
Hi, Canesin,
I am interested in using incompressible solver for fluid in stead of the compressible solver in the current chtMultiRegionFoam solver. I am wondering if you can send me your solver so that I can learn from you? What is your reason for using PISO instead of PIMPLE? Pei phsieh2005@yahoo.com Quote:
|
Quote:
I`m really sorry for not being able to send my solver.. I have an confidentiality agreement with the sponsor of my project .. But as soon as I have submited an paper to some journal I can publish some code here.. I have used an PISO loop because my problem is not steady-state, it has no steady-state configuration, it can be only determined an periodic developed region... In that case I use an PISO loop since PIMPLE is only to permit the use of large time-step, with don't make sence for me.. removing the PIMPLE to a PISO loop I can remove the outer loop (the "nOuterCorr" for loop) and improve the speedy of the solver. |
Thanks Canesin!
Pei |
Hi Maddalena and Andrea,
Could you please show me how to compute the wall heat flux at every iteration when using the chtMultiRegionSimpleFoam (OF 1.7.0)? I am a complete beginner with OpenFoam, and absolutely weak with C++. Thank you very much. Best Regards, Stefano |
Quote:
I commented the if selection while leaving uncommented the inner test in solidWallHeatFluxTemperature. Something like: Code:
// if (debug) mad |
All times are GMT -4. The time now is 10:25. |