
[Sponsors] 
June 7, 2011, 15:28 
conjugateHeatFoam error in new commits of 1.6ext

#1  
Senior Member
Ben K
Join Date: Feb 2010
Location: Ottawa, Canada
Posts: 140
Rep Power: 11 
When I try to run the conjugateHeatFoam tutorial in 1.6ext, I get the following error message in log.conjugateHeatFoam:
Quote:
commit 2f89ab8b718f44bda09af0ed421ec6751d7c52fa Date: Fri May 27 20:53:20 2011 +0100 For the version that doesn't give me a problem, git log gives: commit 3b1cb74951fe1630a2b3f39fc8d267622421f54d Date: Fri Nov 26 11:08:01 2010 + 0000 

June 8, 2011, 11:34 

#2 
Senior Member
Ben K
Join Date: Feb 2010
Location: Ottawa, Canada
Posts: 140
Rep Power: 11 
Just for everybody's info, to solve this problem, you need to add
Code:
. . . attached on attachedWalls off 

June 9, 2011, 06:45 

#3 
New Member
babak kamkari
Join Date: Dec 2010
Posts: 26
Rep Power: 7 
p { marginbottom: 0.08in; } Hi Ben
I am also using OF 1.6.ext. There is something wrong with conjugteHeatFoam solver. If you select different transport properties for solid and liquid regions ( for example DTsolid=1e2 & DTcavity=1e5) you will find that the temperature calculated for adjacent node to the coupled boundary in cavity region is not correct. The temperature of adjacent node is almost the same as the value of interface (coupled boundary)!!!!! when I investigate the temperature distribution through the solid and cavity region, I found that the solver use the solid value of DT (DTsolid=1e2) for calculation of temperature for adjacent node to the interface. Or it seems as if the boundary has moved through the cavity region when the solver wants to solve energy equation I also should mention that this does not happen for calculation of velocity and it is OK. I think there must be something wrong in the solver when it attach the patches for solving energy equation!! I also applied your changes (attachedWalls off) but I didn't get any changes. This defect shows its importance when someone wants to calculate the gradient of transport properties on the interface which will faces to wrong results. 

June 9, 2011, 09:02 

#4 
Senior Member
Ben K
Join Date: Feb 2010
Location: Ottawa, Canada
Posts: 140
Rep Power: 11 
kamkari: I don't think there's a problem. I've been using this solver successfully for a while.
If you can make a simple 1D model and send me the results, I'll compare your results against the same model in Comsol. Make the following model using the coupledFvScalarMatrix functionality, then send me the results (concentrations and the gradient of the concentration): Region 1 (x=0 to 0.5): laplacian(Dleft,c) + 0.1 = 0 BC: at x=0 c=100 Dleft=10^5 Region 2 (x=0.5 to 1.0): laplacian(Dright,c) + 0.1 = 0 BC at x=1 c=0 Dright=10^3 

June 9, 2011, 16:31 

#5 
New Member
babak kamkari
Join Date: Dec 2010
Posts: 26
Rep Power: 7 
Hi Ben
Thank you for your kind offer. I solved the diffusion equation with the mentioned boundary conditions (as the first step without source term). Region 1 (x=0 to 0.5): laplacian(Dleft,c) = 0 BC: at x=0 c=100 Dleft=10^5 Region 2 (x=0.5 to 1.0): laplacian(Dright,c) = 0 BC at x=1 c=0 Dright=10^3 As I expected the value of adjacent node to the interface (left region) has been calculated wrong. It follows the same slope as the right region!!! As I said before the solver seems to use the right regionDT value for calculation of the first node of left region!!! I also emailed you the solved case (limitation for upload). Please have a look and let me know if I have made any mistake. I used the conjugateheatFoam solver OF1.6 ext without any manipulation. Thank you in advance 

June 9, 2011, 18:41 

#6 
Senior Member
Ben K
Join Date: Feb 2010
Location: Ottawa, Canada
Posts: 140
Rep Power: 11 
Ok, I see what you mean. I'm getting a similar result as you and the comsol solution is right.
Here's a link to the solver that I'm working with: http://dl.dropbox.com/u/12812/OF/multiRegion2.tgz I've been burned before using silly settings in fvSchemes...I wonder if that's the problem. 

June 9, 2011, 22:57 

#7  
Senior Member
Ben K
Join Date: Feb 2010
Location: Ottawa, Canada
Posts: 140
Rep Power: 11 
kamkari: I've been playing around with the laplacianSchemes and I've found 2 that seem to work better (but not perfectly) for this 1D example (but I have no idea what these schemes are doing).
Try using: Quote:


June 10, 2011, 08:59 

#8 
New Member
babak kamkari
Join Date: Dec 2010
Posts: 26
Rep Power: 7 
Hi Ben
I solved the same problem using Gauss localMax corrected scheme for laplacian. I calculated concentration gradient on the interface (for left side): linear: 4.28 localMax: 159.04 The correct value: 206.23 Although it seems using localMax scheme leads to better results in this case, there is still about 20% deviation from correct value which is a considerable error!!!! In addition the value of interface will change using this scheme. The final goal of any heat or mass transfer analysis is usually calculation of nusselt or Sherwood number on the interface. But this solver cannot provide good results. I really can’t understand where the problem is. Is it from detaching or attaching the mesh or fvSchemes!!?? Regards, 

June 10, 2011, 12:58 

#9 
Senior Member
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,769
Rep Power: 21 
Easy: you are doing something wrong. The solver works, it is validated and the results are correct.
My guess is that you messed up the diffusivity update or boundary conditions. Hrv
__________________
Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk 

June 12, 2011, 00:39 

#10 
New Member
babak kamkari
Join Date: Dec 2010
Posts: 26
Rep Power: 7 
I checked boundary conditions for coupled patches again and again. They are the same as the tutorial case (heated cavity) with different DT for left and right regions:
right { type regionCouple; nFaces 1; startFace 60; shadowRegion solid; shadowPatch left; attached on; } left { type regionCouple; nFaces 1; startFace 59; shadowRegion region0; shadowPatch right; attached on; } And the same for T boundary condition. So this couldn’t be because of boundary condition. I don’t understand what you mean with diffusity update. As I mentioned before I didn’t changed anything in the solver (OF1.6 ext) And I think the conjugateHeatfoam included in the OF1.6 ext is the last version of the solver. 

June 12, 2011, 09:46 

#11 
Senior Member
Ben K
Join Date: Feb 2010
Location: Ottawa, Canada
Posts: 140
Rep Power: 11 
Hrv means the part where you have to do a: diffusivity.correctBoundaryConditions(); I've tested all of this and it's not the problem. I invite anybody who has 15 minutes to please look at my very simple test case that I posted above and let me know what I'm doing wrong. I've even included instructions for how to run/compile the solver:
http://dl.dropbox.com/u/12812/OF/multiRegion2.tgz In my opinion, the problem is in the fvSchemes. The conjugateHeatFoam tutorial is mostly meant as an example of how to use the coupledMatrix functionality. Whenever there's a sharp change in the diffusivity at an interface, harmonic schemes should probably be used. That being said, I haven't had any luck using harmonic schemes in 1.6ext (see below). I did some more comparisons, this time between OF1.5dev and OF1.6ext. These are my findings: 1. I still can't get exactly the right answer at the interface. 2. As can be seen in the attached figures, something has changed in the harmonic scheme between OF1.5dev and OF1.6ext. To get these results, I used the following fvSchemes: Code:
divSchemes { default Gauss linear; } laplacianSchemes { default Gauss harmonic corrected; } interpolationSchemes { default linear; } 4. Using 1000 grid points of course gives a better solution than 100 grid points as long as you use localMin scheme. So far, the only suggestion I can give is to refine the mesh close to the interface and if you're using 1.6ext use the localMin scheme. This shouldn't be necessary, but it's the best I've been able to do so far. One more comment: I admit that I could be doing something wrong, but it seems to me that the gradient is being conserved at the interface, not the flux. When I do a print out of fvc::snGrad(cLeft) and fvc::snGrad(cRight), I get 99.96 and +99.96 at the interface. Clearly if the diffusivity is changing at the interface, the gradients at the interface should not be the same. The condition at the interface should be: cLeft = cRight AND Dleft*grad(cLeft) = Dright*grad(cRight) Last edited by benk; June 12, 2011 at 16:27. Reason: one more thought 

June 13, 2011, 15:44 

#12 
Senior Member
Ben K
Join Date: Feb 2010
Location: Ottawa, Canada
Posts: 140
Rep Power: 11 
For the record, using the harmonic scheme seems to be an issue with earlier builds of 1.6ext, but not newer ones. With a build date of Feb 4 I have issues (this is the build date of the mac binary) but not with the latest build of 1.6ext.


June 16, 2011, 07:09 

#13 
New Member
babak kamkari
Join Date: Dec 2010
Posts: 26
Rep Power: 7 
Hi Ben
As you suggested I started playing with fvScheme to achieve better results than localMax/Localmin on the interface for the previous case with 40 nodes. I faced something strange when using snGrad with localMax scheme and linear scheme. As far as I know snGrad uses the following formula for calculating the gradient of variable normal to the surface: (Concentration of the interface – Concentration of the neighbor cell center) /dx (1) I implemented the following line in the code for snGrad calculation on the both side of interface: ****** Right side *************** label patchi = mesh.boundaryMesh().findPatchID("right"); gradT.boundaryField()[patchi]=T.boundaryField()[patchi].snGrad(); ****** Left side *************** label patchj = solidMesh.boundaryMesh().findPatchID("left"); gradTsolid.boundaryField()[patchj]=T.boundaryField()[patchj].snGrad(); Here is the comparison between the result of snGrad with formula (1): ****** linear scheme*************** Left Side snGrad : 40.8 formula (1): 42.8 Right Side snGrad : 2016 formula (1): 38.9 ****** localMax scheme*************** »Left Side snGrad : 1016 formula (1): 1590 Right Side snGrad : 1970 formula (1): 443 I expected that the snGard and formula (1) results to the same answer. To check the implementation of snGrad I added the snGrad as above to the boussinesqBuoyantFoam solver (solving natural convection in a single domain) and the results of snGard and formula (1) were the same. Would you please check it and let me know if I am doing something wrong. Tank you in advance Last edited by kamkari; June 17, 2011 at 02:30. 

June 20, 2011, 10:58 

#14 
New Member
babak kamkari
Join Date: Dec 2010
Posts: 26
Rep Power: 7 
p { marginbottom: 0.08in; } I think that I made a mistake for calculating the sngrad on the right region. The correct code for right region is:
label patchj = solidMesh.boundaryMesh().findPatchID("left"); gradTsolid.boundaryField()[patchj]=Tsolid.boundaryField()[patchj].snGrad(); now the results are as the following: Left Side snGrad : 40.8 formula (1): 42.8 Right Side snGrad : 40.8 formula (1): 38.9 ****** localMax scheme*************** »Left Side snGrad : 1016 formula (1): 1590 Right Side snGrad : 1016 formula (1): 443 There is still difference between formula (1) and snGrad. I found something interesting. SnGrad unexpectedly uses two point on the right and left side the interface for calculation of the surface normal gradient!!!!! This couldn't be right, because it should use the value of interface and its adjacent node value. 

June 20, 2011, 16:46 

#15 
New Member
babak kamkari
Join Date: Dec 2010
Posts: 26
Rep Power: 7 
Hi Ben
I think that the snGrads you found for right and left regions in your last post are not correct, since the sngrad surprisingly is using two nodes located in right and left of the interface (you can check it by hand calculation)!!! Obviously if you use snGrad for linearMax /linearMin it shouldn’t leads to the same value for right and left regions but snGrad gives the same values. I tried to solve it by detaching the patches: # include "detachPatches.H" label patchi = mesh.boundaryMesh().findPatchID("right"); gradT.boundaryField()[patchi]=T.boundaryField()[patchi].snGrad(); label patchj = solidMesh.boundaryMesh().findPatchID("left"); But I get the following error: > FOAM FATAL ERROR: Requesting patchToPatchInterpolation in detached state From function const patchToPatchInterpolation& regionCouplePolyPatch:atchToPatch() const In file meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.C at line 118. FOAM aborting Regards 

June 20, 2011, 20:30 

#16 
Senior Member
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,769
Rep Power: 21 
Hahaha, you really made me laugh: so the gradients are correct, just like the doctor says.
As for your formula, you are wrong:  in the detached state, the gradient is face centre minus cell centre divided by the distance, just like on any normal patch  in the attached state, the snGrad is right CELL value minus the left CELL value divided by the distance, just like on any coupled patch. So, as I said all along, the code is exactly right, conservative and consistent. Any other behaviour would be inconsistent. Hrv
__________________
Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk 

June 21, 2011, 15:17 

#17 
Senior Member
A_R
Join Date: Jun 2009
Posts: 118
Rep Power: 8 
Dear kamkari
You mean that snGrad() uses central difference scheme(but the openfoam site says different definition). So, is it really normal to surface? The second problem is there are two regions having different conductivity. Can the snGrad() understand it? I think the source of the problem is the solver. It shouldn`t really differ between two region. It is a really childish bug. 

June 21, 2011, 16:26 

#18 
New Member
babak kamkari
Join Date: Dec 2010
Posts: 26
Rep Power: 7 
Dear Prof. Jasak
First I have to thank for your time and participation in this discussion. You still insist that the code is exactly right but I think it worth putting a little time for modification of code for the next extended version (1.7ext). There are some notes: 1 Although your solver (conjugateHeatFoam) is easier for handling and may be faster than chtmultiregionFoam but without the capability of right prediction of gradients at the interface it can’t be much useful for heat or mass transfer problems. 2 I have used snGrad() in attached state and resulted to the same value for both regions. It means noting for a heat transfer problem with different DT for regions. 3 I am sure that what Doctor Ben has calculated for Right & Left gradients was in attached state which obviously leaded to the same problem and it is not what we are looking for. 4 If someone could calculate the snGrad() on the interface in detached state and get reasonable values for interface gradients, I mean that DTleft/DTright=snGrad()right/snGrad()left, you could be confidently sure that the code is exactly right. 5 I have calculated the gradient on the interface for both regions for the simple 1D case suggested by Ben like this with a calculator: Gradient of the variable on the interface for Right/Left region = (face center value – cell center value)/dx The results were far from the result of other software of analytical solution !!! 6 Whatever happens in the solver is related to the interpolation on the interface because as the Ben said by selecting localMax or Loacl Min for laplacian scheme the results get better but still not acceptable (60%100% deviation from the right value of interface gradient). The attached pics of previous posts show this clearly. Regards 

June 22, 2011, 00:25 

#19 
New Member
Hadi
Join Date: Nov 2009
Location: Iran
Posts: 2
Rep Power: 0 
Dear Babak,
I have recently installed OpenFOAM 1.6ext and faced with the same problem as you mentioned. I think that your evaluation about the code is true. More investigation about this bug is required by the developers. Good luck 

June 22, 2011, 04:35 

#20 
Senior Member
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,769
Rep Power: 21 
This is starting to annoy me  make sure you don't show up for an exam at my class because you're going to fail.
1 in the attached state, the gradient across the Coupledinterface equals right minus left CELL value and includes the nonorthogonal correction and all other goodies 2 in the detached state the face minus cell gradient applies. I have state this in the snGrad message above 3 the flux is calculated AT THE INTERFACE, and is a product of the snGrad and the diffusivity at the face. The diffusivity at the face is calculated to be consistent on both sides 4 how to calculate diffusivity: in case of jump in diffusivity in two materials, the FACE diffusivity is calculated based on serial resistances of each material. This implies harmonic averaging at the interface, achieved by the regionCoupling b.c. ON THE DIFFUSIVITY VARIABLE. if your FACE diffusivity is. Not identical on both sides, you are doing something wrong  probably not correcting the diffusivity b.c. correctly. 5 in case of ADDITIONAL RESISTANCES AT THE INTERFACE, such as wall functions from turbulence model, insulation or similar, further serial resistances are added AT THE INTERFACE ITSELF. Note that cell values of diffusivity are still material properties of each material, including any jump states that may occur. This is handled specially, by manipulation of surface diffusivity in the laplacian. 6 in case of radiation or similar surface effects, further treatment is needed  but I told you enough already. Is that now clear? The code is correct, it has been validated and it works properly. To conclude with a personal note: I am a University professor and a professional CFD developer who wrote a million lines of quality CFD code that you can look at and check yourself. I DO KNOW how to implement a region coupled boundary condition (and you don't). So the next time you get into this situation ask yourself are you really certain enough in what you are saying to tell Hrv he does not know how todo his job. If you are, and you happen to be wrong, I will come down upon you like a ton of bricks! So let's hear it now: go to your diffusivity fields on both sides of the interface, look at the face values and tell me if they are the same. If not, why didn't you update them correctly with a correct boundary condition? Hrv
__________________
Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Test directory missing in OpenFOAM 1.6 ext  andrewryan  OpenFOAM  2  March 20, 2011 16:41 