[mike_jaworski]

Update with SVN revision 492:

Errors above are unchanged.

Regards,
Mike J.

[hjasak]

I was hoping you had some basic C++ knowledge, which would allow you to fix this yourself, based on the snippet above... In any case, I have now done it myself and checked it in.

Hrv

[mike_jaworski]

Hrv,
"Basic" is a relative term. I'm learning both C++ and OF at the same time. Your help and patience is greatly appreciated.

Regards,
Mike J.

[mike_jaworski]

Hrv,
The updated solver compiles, and I was eager to look at what the missing steps were and run it through the benchmark cases I have, however it produces a segmentation fault.

Now I'm feeling a bit awkward because I keep having difficulties and am not much help. However, I have isolated the problem by elimination. You can comment out everything and add parts in and see where the the fault occurs. I did this and it's the rcp.detach(); function call.

I'm going to try and understand the functions in regioncouplepolypatch and compile the debug version tomorrow (need to sleep tonight). Maybe a solution will present itself though on a faster time-scale than the (steep) learning curve I'm on.

Regards,
Mike J.

[hjasak]

Heya,

I have first tidied up the signalling on coupled/decoupled state of region-coupled patches and fields - this gets rid of segmentation fault (please accept my apologies).

I have then had a long careful look at this issue of non-zero vectors on a coupled wall. The conclusion is that cell-centre and face-centre fields are all fine and what you are looking at is basically due to point-based post-processing. let me explain:

- paraview and other similar tools want point-based fields. In order to do the best possible job, OpenFOAM does the intepolation

- for internal points, interpolation uses cell values only (inverse distance weighting), which would give strange-looking values on the boundary

- to overcome this, after the cell-to-point interpolation is done, I visit all boundary points and do the interpolation FROM PATCH FACES to points, which over-rides the first attempt. This is how you get a zero on the wall for the velocity

- you should not be doing this correction if the patch is coupled: in that case, interpolated cell values give you the right number.

- what about a region-coupled patch? Well, it can be both coupled (use interpolated cell values), or decoupled (over-ride with face values). Therefore, I cannot tell which to do: for temperature, I want it coupled and for velocity decoupled (in fact, it's no that bad because the face values for temperature are identical on both patches).

- how does the post-processor know if the patch is coupled or not? Well, it will read constant/polyMesh/boundary file and look for:

right
{
type regionCouple;
nFaces 10;
startFace 200;

shadowRegion solid;
shadowPatch left;
attached on;

}

This one is attached, I do interpolation without over-ride and you get non-zero vectors on the wall:



If I tell the mesh it is detached:

right
{
type regionCouple;
nFaces 10;
startFace 200;

shadowRegion solid;
shadowPatch left;
attached off;

}
I will do an over-ride with wall values and get a zero velocity on the wall:



Note that I am showing unscaled vectors and that an x-facing blue arrow basically means (0 0 0). The code is now attaching and detaching in each time-step, so what you have in the dictionary is just a "static" state.

So much from me: all you have been chasing is a bit of a strange feature in interpolation from the correct CFD solution. Quite interesting though - I've enjoyed my dif through the code

In summary: there's nothing wrong: this whole attach-detach business isn't really necessary in terms of discretisation/numerics accuracy, but at least we now have nice piccies!

Enjoy,

Hrv

[mike_jaworski]

Hrv,
First, this is really good to know about what's going on in post-processing and the nice pictures are not to be under-estimated, in my opinion. I had shown a picture with this sort of odd behavior in it at my preliminary examination. I had to explain that, at the time, I thought it was post-processing but would benchmark the code to be sure. So, here I am, but it avoids all sorts of questions in presentations, etc.

Second, I really do appreciate all this time you're putting in, just as I hope that the benchmarking I'm doing is bettering the code in the way I can contribute right now.

So, now we come to it: that we're not quite out of the woods yet. Here are the results of more benchmarking I have done. They show that the conjugateHeatFoam solver is not as accurate as a modified icoFoam but that this is possibly due to calculations involving the DT field.

I am using the Blasius analytical solution for boundary layer growth as a benchmark for all these cases.

I model this with a rectangular region with a short entrance length. The entrance has an inlet velocity of 0.2m/s, the top wall moves at that same velocity and the entrance length has zeroGradient applied. The flat plate is a no-slip wall with 0 velocity. The outlet is fixedValue pressure at 0. The fluid region is kept exactly the same for all these comparisons. In the case of the conjugate solver, there is an additional solid mesh at the plate boundary which is 1mm of material. The fvSchemes and fvSolution files were matched for these runs.

Here are the case files for inspection:
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Blasius-flatPlate/icoHeatFoam -flatPlate2.tar.gz
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Blasius-flatPlate/conj_flatPl ate.tar.gz

First, the good results come from a modified icoFoam which tracks temperature which I called icoHeatFoam for these purposes.
Here are two graphs showing the comparison to velocity and dimensionless temperature (T(eta) - Te/ (Tw-Te)) where Te is the free-stream temperature and Tw is the wall temperature.
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Blasius-flatPlate/compareVel_ icoHeatFoam.png
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Blasius-flatPlate/compareThet a_icoHeatFoam.png

These are really good results, in my opinion. There is a slight flow acceleration due to the growing boundary layer and continuity, however, the temperature solution is nearly exact. I have run this same case with a more refined mesh and the theta solution does become more exact with sufficient cells.

In order to approximate the constant temperature boundary used in the Blasius solution, I modeled the solid as a very thin(1mm) block of material with the material properties of copper. There are constant temperature boundaries on all sides of the block except, of course, the regionCoupling BC. Before using the sample utility, I altered the boundary file to disconnect the regions so it would output the correct values.
Here are the comparisons of this solution with the analytical one:
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Blasius-flatPlate/compareVel_ conjugateHeatFoam_Cu.png
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Blasius-flatPlate/compareThet a_conjugateHeatFoam_Cu.png

The velocity profile is not quite as exact as icoHeatFoam gives, but it is not way off the mark. It is much closer to the analytical than before the attach/detach business, even with the "static" state item we have found in the post-processing (so, good: progress!). The theta solution, is far off the mark here, however. I will note that the temperature in the block does not vary more than a couple degrees C (out of 75C total temperature range) and playing with Tw in the plotting does not sufficiently alter the graphs for this to account.

However, in order to remove that doubt, I altered the material properties of the solid to better approximate a constant temperature wall by making it 100x as thermally conductive as copper. This should have absolutely no effect on the fluid solution since I am not altering those values and the analytical solution only assumes constant temperature wall conditions. More conductive material should better approximate that condition.

Here is the results for velocity and temperature:
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Blasius-flatPlate/compareVel_ conjugateHeatFoam_100Cu.png
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Blasius-flatPlate/compareThet a_conjugateHeatFoam_100Cu.png

As expected the velocity solution is not altered at all by the change in thermal material properties. The temperature solution, however, is radically altered! Again, the change in thermal diffusivity of the solid should not have altered the temperature solution in the fluid since the constant wall temperature BC was already approximated by copper or a more conductive (but still finite) material.

It seems to me that there must be some calculation being made dependent on the DT fields that is affecting the solution in a way that it ought not to. I have not looked at this machinery yet, but the behavior of the solver with these adjustments to material properties point me in that direction.

In conclusion, I am heartened not just by the progress in my own understanding of OpenFOAM, but that the velocity issues are mostly fixed and that is a step forward. It is also great to see OpenFOAM achieve a benchmark as closely as it does with the icoHeatFoam solver. The temperature solution for the conjugate solver, however, seems to be tied to the DT fields in such a way that unphysical results are being produced.

Because of some of the modeling issues with making a sudden transition to a flat plate, I intend to create a stagnation flow benchmark, again based on the Blasius analytical flow solution (I follow Frank White's "Viscous Fluid Flow" 3rd Ed. 2003). I think this will complement things nicely and remove some of my last doubts due to the flow acceleration issue, etc.

Regards,
Mike J.

[jwp]

Hi!
I have a problem with the conjugateHeatFoam-solver (or with me?)
My case is a plate of copper as the solid and water as the fluid. There is a internal spring of heat in the solid. The fluid flow through the solid like a serpentine and normally it's fine. But the heat transfer through one boundary isn't correct.
The interesting is, that all boundaries have the same conditions (and the same name) and the mesh is okay! (I have created a new mesh and another boundary is now not okay - and the former is okay)

Any ideas? I could show a picture but I don't now how to upload it.

Jens

[mike_jaworski]

Jens,
Hrv discovered that some of the post-processing steps can result in visual artifacts. First, make sure you have the latest subversion revision (512 came out several days ago). Also, edit your boundary file so that the regionCouple patches are not attached:

shadowRegion a_Region;
shadowPatch a_Patch;
attached off;

This will avoid the interpolation problems that can make the data look funny in a picture. The latest release of the conjugateHeatFoam solver won't complain if you leave the boundary as "off" normally. You'll have to do this for all the patches that are of type regionCouple.

There may be other quantitative problems with the simulation and this should at least let you know if that is occurring.

Regards,
Mike J.

[jwp]

Thanks Mike!

Hmm, the problems are in the heat transfer, so I thought, I had to use attached on.
But in fact, I didn't use the latest subversion.

[jwp]

I've solved the problem by defining all coupled patches separately. Before that I had only one patch for the solid interface and one for the fluid. I don't know, whether it didn't work, but now it works fine!

Jens

[mike_jaworski]

I have constructed another benchmark case for this solver and it shows the same odd behavior with respect to changing solutions when it should not be changing.

Introduction:
The case is a planar stagnation flow for which a similarity approach can be taken to describe the boundary layer that forms at the plate. The typical approach is to use the inviscid flow solution to yield the free stream velocity beyond the boundary layer, and then normalize the boundary layer velocity to the free stream. The analytical solution for the temperature profile assumes a constant temperature wall, as with the Blasius solution shown above.

The results of the case study will be broken down similarly to the flat-plate flow solution given earlier: first a demonstration of the correct solution behavior with a non-coupling code and second, the results with conjugateHeatFoam. The case files are contained here:
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Similarity-Stagnation/icoHeat Foam-planarStagnation.tar.gz
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Similarity-Stagnation/conjuga te-planarStagnation.tar.gz

Results:
First, the simply modified icoFoam (icoHeatFoam) provides good velocity and temperature agreement to the analytical solution. The overall solution velocity field is shown in the following figure:
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Similarity-Stagnation/icoHeat Foam-planarSolution.png
and quantitative comparison is shown in these two figures:
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Similarity-Stagnation/compare Vel_icoHeatFoam_planar.png
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Similarity-Stagnation/compare Theta_icoHeatFoam_planar.png

The conjugateHeatFoam solution provides similar results for the material properties of water on copper (with a single fixedValue temperature boundary condition on the bottom of the block. The overall temperature solution is shown in this figure:
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Similarity-Stagnation/conjuga teHeatFoam-planarSolution.png
The quantitative comparison with the analytical solution is shown in the following two graphs:
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Similarity-Stagnation/compare Vel_conjugate_planar.png
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Similarity-Stagnation/compare Theta_conjugate_planar.png
which are subtly different from those created by the icoFoam based solver. They are, however, fairly accurate.

In the previous test, a certain non-physical behavior was observed in the solution. Namely, that the temperature solution was affected by altering the material properties of the solid in the direction of more conductive. This is a non-physical result because the analytical solution assumes a constant wall temperature. A MORE conductive medium (high diffusivity) is more akin to a constant temperature situation as it can maintain smaller and smaller temperature gradients across it for the same power flux. However, as was seen in the case of the flat-plate flow, the solution is again modified by making the plate material more conductive as shown here:
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Similarity-Stagnation/compare Theta_conjugate_planar_Cux100.png

Discussion:
Again, OpenFOAM shows it can provide accurate solutions in the benchmark cases for a simply-modified icoFoam solver. This is clearly shown in the icoHeatFoam solver runs. However, the conjugateHeatFoam solver appears to suffer from some bugs which prevents it from providing reliable results. The coupling between the solid material and the liquid is affecting the solution in a non-physical way and because of this, the apparent similarity to the temperature solution seems like happy-chance. The flat plate flow case showed more definitive differences with the copper plate material, however both solutions exhibit the non-physical dependence on the plate material property.

Conclusions:
I am going to look at the conjugateHeat solver that Danielle wrote a while back which iteratively solves the temperature fields in the two regions. I'm not at a point yet where I can debug the equation coupling (coupledFvScalarMatrix) portions of OpenFOAM, though I do feel comfortable at the higher-level of things with the iterative method from what I've seen so far.

An alternative is for someone to produce a range of applicability for the solver which can be tested and the accuracy determined against analytics or other benchmarks.

Regards,
Mike J.

[hjasak]

You do not need to look at the fvScalarMatrix, all is well there. I do however have a question and an idea.

The only possible source of the difference is the value of diffusivity in the energy equation right on the fluid-solid interface. I've been thinking about this for a while and it is not clear if the answer is the same for all fields and all solutions.

Basically, I currently use linear interpolation to evaluate the value of DT on the interface, which from your results is clearly wrong. I suspect harmonic interpolation should be the right way, but that will not work for vectors and tensors. This is very much in line with your results, where the neat transfer is massively over-estimated for a large jump in diffusivity/conductivity.

From your tests, can you tell me what is the value of diffusivity on the wall face:
- on the wall?
- separately on the fluid and solid side?

This is an easy fix, but I need to know what's the right thing to do...

Please let me know,

Hrv

[jwp]

Hi!

Do anybody know, whether it's okay, when I use a steady state flow by commenting out the line
# include solveFluid.H
in conjugateHeatFoam.C ?
And the other thing is: I want to use an internal source by adding to the solid equation, like:


TEqns.set
(
1,
new fvScalarMatrix
(
fvm::ddt(Tsolid) - fvm::laplacian(DTsolid, Tsolid) - ( I /(Ksolid/DTsolid) )
)
);

with I in W/m³

It is okay?

Jens

[stephan]

Hi Hrvoje,

while calculating with your solver (still nice!) almost the same problem as michael experienced came up. using 1.0*10^-4 on one side and 1.0*10^-7 on the other side of a coupled patch gives huge differences between geometric and linear interpolation.
i am not sure what to do with tensors and vectors but i think for most problems geometric interpolation for scalar should be fine.
could you please give me a hint where i have to look in the sources in order to change the interpolationscheme to geometric.
thanx in advance
best regards
stephan

[jwp]

Hrv, I think I have similar problems like you wrote.
I have a jump in the diffusivities of a factor about 1000 (10^-7 and 10^-4) and DT-values on the interface of about 10^-6 till 10^-5
So I think harmonic interpolation would be a better solution, isn't it?
Is there a way to solve this problem? I use isotropic materials, so I have scalar values.

Jens

[hjasak]

Well, here it is. Have a look at

fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C

or in SVN:

http://openfoam-extend.svn.sourceforge.net/viewvc/openfoam-extend/trunk/Core/Ope nFOAM-1.4.1-dev/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/re gionCoupleFvPatchField.C?view=markup

template<class>
void regionCoupleFvPatchField<type>::evaluate()
{
Field<type>::operator=
(
this->patch().weights()*this->patchInternalField()
+ (1.0 - this->patch().weights())*this->patchNeighbourField()
);
}

I'm on line 168: this is linear interpolation.

This is where you can change the interpolation. My problem is that harmonic interpolation is really properly defined only for scalars, and this code is templated. The suggestion I've got is to take a magnitude of the variable, do harmonic average on the magnitude, then conver the harmonic average into weights

phi_f = f_x phi_P + (1 - f_x) phi_N;

where f_x is calculated as

f_x = (phi_harmonic - norm(phi_P)/(norm(phi_N) - norm(phi_P));

and

phi_harmonic = 1/(1/norm(phi_P) + 1/norm{phi_N))

stabilised and weighted if you feel like it.

This also needs escapes when the denominator is zero etc.

So, here's the code, here's the idea. I have no time to implement and test this right now - does any of you (especially the guys who already have validation cases) feel like doing it?

Please let me know + "please contribute"

Hrv

[stephan]

Hi Hrvoje,

actually Jens needs this pretty soon and i am quite busy till end of next week - so i guess he will do it. (we working at the same university) I cannot reach him rightnow but in case he wont i will do it (or let is say this: we will both try )
best regards
stephan

[stephan]

Hi Hrvoje,

one thing was lost in the last posting:
thanks a lot for fast respond and your ideas!
br
stephan

[hjasak]

Well, I want it as well, but life is hectic. Whoever moves first, please keep me posted.

Hrv

[mike_jaworski]

Hrv (and others),
I started trying to pull the diffusivity out of the solutions I had. I just downloaded the updated svn and I'm compiling now, so I can test the updated libraries/etc. It'll be later today/tonight when I can recheck it.

As for checking on the behavior of the diffusivity on either side of the boundary, it *appears* that this is, indeed, being evaluated incorrectly. I'm a little less than certain, however, due to the issues of how the sample utility might be interpolating the fields... though I suppose this is all part of the overall problem. Still, like I said, I'm a little less than certain, though I did make sure the regionCouple patches were all "off" before running the utility and this shouldn't make any difference at all in the case of the temperature field. The solid side looks like this:
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Similarity-Stagnation/solidEx amineTemp.png

The solution *should* be a simple linear line since it's hit steady state and there *shouldn't* be a discontinuity in the slope from material properties. Interestingly, there are 4 cells across the bottom 1 mm of material which indicates that the diffusivity is being evaluated incorrectly over half of the cell (from center to interface).

The diffusivity is evaluated as 6 times *less* than the rest of the domain. There's a similar discontinuity on the fluid side. and it's being evaluated as *more* thermally diffuse.

Again, I'm cautious about this analysis since I'm still fuzzy on some of the interpolating aspects of the sample utility and others.

As for extending this to vectors and tensor quantities, I have a thought on how to test this. I think a solver can be worked up to utilize anisotropic conductivity. I can do this and create the test case. I can create the top-level solver, the interpolation methods of OF is something I'm still learning, truth be told.

Ah but looking at the your post, it doesn't look like it's been updated yet. Still, good to have the freshest svn I suppose!

-Mike