[hjasak]

Hello,

Fields to couple:

Clearly, you must couple the solution variable, say temperature or energy. You also need to couple diffusivity (look at your laplacian): this is because both sidesw need to "see" identical diffusivity on the interface. An "early" error I did was to use a dimensionedScalar as diffusivity - the code runs, but the result is wrong because the left side sees one value and the right a different one.

There is no need to couple rhoCp - it only appears as a cell-based term (in ddt) and not in flxues.

Forget about the manual: the time-step size depends only on your convection term. If you have diffusion-on-diffusion, you can jump straight into steady-state.

Enjoy,

Hrv

[mike_jaworski]

Hrv,
Thanks for the information on the time-step material. This answers my question.

About the diffusivity coupling I have an additional question:
If one allows for spatial and time variation of rho*Cp quantities, the PDE is given as:

d/dt ( rho * Cp * T ) + d/dx_i (u_i rho * Cp * T) = d/dx_i ( k * d/dx_i ( T ) ) + some other terms

If one assumes constant rho and Cp, then the problem reduces, I believe, to the one you wrote where DT = k/(rho*Cp) and DT can vary with position by varying k only. If, however, rho and Cp are spatially varying, then they have to remain in their respective terms above.

My question is that with the spatial derivative in the flux term, does this make it necessary to couple rhoCp in addition to k? I'm not sure because the rhoCp flux term across the solid-liquid interface should automatically be zero with a u=0 boundary. But I think the question is relevant since this isn't the most general coupling one can imagine.

What do you think?

Regards,
Mike

[hjasak]

Easy. For the case of conjugate heat transfer, there is no flow through the solid-fluid interface: it is a wall, right? Thus, u = 0, which sort of kills your problem.

If you have a more complex interface, you will get a flux from somewhere, meaning a pressure equation. In that case, it will be the pressure that will be coupled and you will get the consistent flux on both sides naturally.

Clear?

I actually wrote a coupled fluid-fluid solver, but that's another story...

Enjoy,

Hrv

[mike_jaworski]

Hrv,
Thanks again for enlightening me. Being a beginner to a lot of this (including coding), OpenFOAM is a bit daunting and intricate to take in all at once.

-Mike

[mike_jaworski]

Hrv,
I think I've found a bug of sorts in the solver. I noticed this in my own edited/altered solver based on conjugateHeatFoam, but checking things more closely with the tutorial file, I see the same thing.

Basically, it appears that the velocity boundary condition of (0, 0, 0) velocity on the region coupling wall is not being applied correctly.

conjugateCavity velocity solution:
https://netfiles.uiuc.edu/mjaworsk/s...teCavity_U.png

and if you take the exact same mesh and U boundary conditions and run it through icoFoam, you get:
https://netfiles.uiuc.edu/mjaworsk/s...gateTest_U.png

Looking at the two files there appears to be a non-zero velocity on the right wall (which is the region coupling BC in conjugateCavity) whereas the appropriate solution with icoFoam shows zero velocity there.

This same behavior showed up in a larger run I had performed and resulted in anomalous heating of a material there. I thought it was a visual artifact from the FoamToVTK problems, at first, but dxFoam is showing the same thing. My solver is only slightly edited so I went back and checked the original and it's coming up there too.

The test case I ran which caught my attention is shown in the following two images:
https://netfiles.uiuc.edu/mjaworsk/s...hium-5mm_U.png
https://netfiles.uiuc.edu/mjaworsk/s...uid-test_U.png
The far right and bottom boundaries should be zero velocity as shown in the "test" file. those two boundaries are also temperature regionCoupling (as well as k/DT) and they appear to have non-zero velocities. Additionally there's anomalous heating there, but I'd lay that up to this odd velocity.

Maybe you already saw this, but it appears to be some sort of bug.

Regards,
Mike

[hjasak]

Hello,

Yes, I know what is going on. In short, the regionCouple patch exists in two forms: coupled and uncoupled. Each of those comes with its own set of geometric ocefficients, which means you have to tell the couple and decouple itself. Thus, there is a function in the regionCouple, which says:

//- Attach regions
void attach();

//- Attach regions
void detach();

Once you attach or detach the couple, you need to force the re-calculation of weights and interpolation coefficients.

The pictures you are showing say that you are trying to solve the momentum equation with boundary conditions on the conjugate boundary, but you didn't detach the patch.

Hrv

[mike_jaworski]

Hi Hrv,
Thanks for the reply. So it look slike a detach() command needs to be given before the momentum equation is solved and then attach() needs to given before the energy equations are solved?

-Mike

[mike_jaworski]

Hrv,
Sorry again: is there an example on the usage of attach and detach... or how should the conjugateHeatFoam solver be altered so that it correctly solves the momentum problem? I'm still learning the ins and outs of C++ and this quickly left my zone of knowledge.

Do I simply make it:

detach();
# include solveFluid.H
attach();
# include solveEnergy.H

??
or does it need something like:
solid.mesh.detach();

??
Like I said, I'm still learning the C++ language so that may all be incredibly stupid, so I appreciate your help.
Also, I'm going to find a suitable benchmark case for the solver, something like Blasius flat plate flow or stagnation flow or the like to compare this against, so this work won't go into a black hole.

Regards,
Mike J.

[panara]

Hi All,

is the integrating solver working also in parallel?
If yes how should be the grid decomposed?

Daniele

[mike_jaworski]

This is a hack, but I'm trying and getting stuck.

I'm trying to decouple and recouple the meshes before solving the fluid equations. This will hopefully correct the present errors in this solver.

I've added to solveFluid.H the following
*********************
#include "coupledFvPatch.H"

label interface = mesh.boundaryMesh().findPatchID("plate");

makeWeights( U.boundaryField()[interface] ) = 1;

// Regular fluid equations go in here

makeWeights( U.boundaryField()[interface] ) = 0.5;

*********************

This however generates the error that 'makeWeights' was not declared in this scope.

Any thoughts on how to fix that? Is this even the right approach? Also, the patchID thing is something I started out doing to test with a specific case. There must be some elegant way of generalizing this...

By the way, regionCoupleFvPatch.H does not contain any "attach" or "detach" functions. where do those live?

[vtk_fan]

Hi Mike,

Were you able to resolve the discontinuity issue in your example studies?

[mike_jaworski]

Hi Mark,
Yes. I believe the first issue had to do with the sample utility and not the solver itself, though I must admit that after determining the temperature and slopes were the correct values, I didn't continue testing it.

As for the velocity problem, this is why I'm trying to decouple the regions before solving the fluid equations.

I have made a benchmark case study using the Blasius flat-plate flow solution as the standard. I modified icoFoam (icoHeatFoam) to perform thermal transport and determined the appropriate solution geometry to get the solution close to the Blasius flow. The constant temperature flat-plate problem using the modified icoHeatFoam is shown here:
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Blasius-flatPlate/gnuTheta.pn g
This is for water (Pr = 7) liquid, theta is the dimensionless temperature difference T-Te/(Tw-Te) and eta is the dimensionless similarity variable.

If you do the same problem, pretty much, with conjugateHeatFoam, however, you achieve this:
https://netfiles.uiuc.edu/mjaworsk/shared/OpenFOAM/Blasius-flatPlate/gnuTheta_co njugate.png

Solution accuracy suffers quite a bit because of the velocity issue at the boundary. I'd like to determine how to fix the conjugate solver and then write a tutorial on the whole thing and submit the case files, but I'm stuck trying to delve into the vastness that is OF. If you can help, that would be appreciated.

Regards,
Mike J.

[hjasak]

Heya,

Sorry, this is really my fault. I was going to write a little piece of code for you but it's just too busy (big things are afoot!)

This is what we will do: couple the mesh before the energy equation and decouple it afterwards. I will write some pseudo-code for you - please test it and let me know.

So, before

# include "solveFluid.H"

do something like:

const polyPatchList& patches = mesh.boundaryMesh();

forAll (patches, patchI)
B
if (isType<regioncouplepolypatch>())
B
regionCouplePolyPatch& rcp = refCast<regioncouplepolypatch>(patches[patchI]);

// Detach it here
rcp.detach();
B
B

//Force recalculation of weights
mesh.surfaceInterpolation::movePoints();


Formatting is eating my brackets so I've used B-s. Also,m beware of capitalisation, that got eaten as well...

Before the energy equation, do rcp.attach() and all should be well.

Looking forward to hearing about your results,

Hrv

P.S. I will pack this up properly into a class or a function, with named patches to be coupled and decoupled etc...

[hjasak]

Forgot: you will have to do the same on the solid mesh before you do mesh.surfaceInterpolation::movePoints();

In any case, tell me if you're stuck and I'll help out.

Hrv

[vtk_fan]

Hi Mike,

I am delving into the OF code organization now, and the internals of the discretization operators. However, it will be at least 2 weeks before I can assist in debugging this problem. Hope Hrvoje's suggestion above settles the issue.

Would you be interested in working on the DOM model? If so, maybe we can collaborate on that.

[mike_jaworski]

Hrv,
Thank you for the snippets. I will be able to try them out this weekend and report back next week on the results.

-Mike J.

[mike_jaworski]

Good news and bad news, I suppose.

First the good: I think I understand more of what is going on in OF here

Bad news: There's a road block error.

code I've added:

const polyPatchList& patches = mesh.boundaryMesh();
const polyPatchList& solidPatches = solidMesh.boundaryMesh();

...
(enter runTime loop)
...

forAll (patches, patchI)
{
if (isType<regioncouplepolypatch>(patches[patchI]))
{
regionCouplePolyPatch& rcp =
refCast<regioncouplepolypatch>(patches[PatchI]);
rcp.detach();
}
}

Those same lines (almost) are repeated for solidPatches

and then mesh.surfaceInterpolation::movePoints();

is put in.

The same loops are run after the fluid equation is solved except they have rcp.attach();

Now, here's the wall I run into:

my_conjugateHeatFoam.C:94: instantiated from here
/home/mjaworsk/OpenFOAM/OpenFOAM-1.4.1-dev/src/OpenFOAM/lnInclude/typeInfo.H:99: error: cannot dynamic_cast '(const Foam::polyPatch&)((const Foam::polyPatch*)r)'(of type 'const class Foam::polyPatch&') to type 'class Foam::regionCouplePolyPatch&'(conversion casts away constness)
make: *** [Make/linuxGccDP0pt/my_conjugateHeatFoam.o] Error 1

There's a similar complaint centering on typeInfo.H:108

Ok, so as I read it, the refCast function doesn't know what to do with the regionCouplePolyPatch class yet? The older conjugate solver had this same error and solved it by editing one of the solver specific files. I was browsing around the typeInfo.H file, however, and I don't see any information that would change how this works? Admittedly, however, this rabbit hole is really deep and I've not looked at this material before.

I wanted to check things by trying a second method for a specific case I created. I know the patch names that are coupled in that case so I tried an alternative route:

*******
label fluidPatch = mesh.boundaryMesh().findPatchID("plate");
label solidPatch = solidMesh.boundaryMesh().findPatchID("plateTop");

regionCouplePolyPatch& rcpFluid = refCast<regioncouplepolypatch>(mesh.boundaryMesh()[fluidPatch]);
regionCouplePolyPatch& rcpSolid = refCast<regioncouplepolypatch>(solidMesh.boundaryM esh()[solidPatch]);

...
//And then inside the runTime loop execute
rcpFluid.detach();
rcpSolid.detach();
mesh.surfaceInterpolation::movePoints();
...
**************************

This, however, results in the same error about not being able to dynamic_cast.

Any ideas?

-Mike J.

[mike_jaworski]

Mark,
Radiation with participating media is an interesting problem, however it's a little outside my direct research focus at the moment. Perhaps if there are areas where it overlaps with the conjugate solver I could provide input such as I'm able(I am a newbie at both OF and C++ ).

In a year's time, I may be more interested in adding volumetric radiation sources to a solver, so there may be something down the way in this, right now, however, I have a pretty full plate (I'm not even a computationalist, I'm an experimentalist!)

[hjasak]

OK, not a problem. I have changed the access function in SVN. In the meantime, you can also try:

const polyPatchList& patches = mesh.boundaryMesh();

forAll (patches, patchI)
B
if (isType<regioncouplepolypatch>())
B
const regionCouplePolyPatch& rcp = refCast<regioncouplepolypatch>(patches[patchI]);

// Detach it here
const_cast<regioncouplepolypatch>(rcp).detach();
B
B

//Force recalculation of weights
mesh.surfaceInterpolation::movePoints();


Tell me what happened,

Hrv

[mike_jaworski]

Hrv,
Thanks for the added suggestions. I added the code as you suggested. There remain compiler errors.. First, the dynamic_cast error remains in there as before(conversion casts away constness). An added error message appears, as follows:

my_conjugateHeatFoam.C:97: error: invalid use of const_cast with type 'Foam::regionCouplePolyPatch', which is not a pointer, reference, nor a pointer-to-data-member type

Do I understand this right that because regioncouplepolypatch is so new, the lists of items like pointers and valid polypatches has to be updated as well, to let them know that it's an available/valid patch type? (I'm still trying to understand the machinery behind the scenes here.)

Since you have updated the SVN, I will update and see if I can find the appropriate subdirectories to recompile (saves a lot of time) and try again.

Regards,
Mike J.