Hi, I am trying to implemen
 

Hi,

I am trying to implement a solver that calculates a scalar potential over two connected regions. At the boundary the scalar potential must be continuous and there is a specified jump in the gradient.

This is very similar to the problem discussed in the thread

OpenFOAM Message Board: OpenFoam: Running / Solving / CFD: "Heat transfer with solid elements (conduction)"

in which a similar two region Laplacian solver was given (conjugateFoam2.tar.gz). However, in that solver the boundary between the two regions was parallel to the x-axis and the problem is one dimensional, while my problem is 2D and the boundary is more general --- it is a circle. Hence I used the following code (my two regions are called 'wire' and 'space', and I am imposing the value of the potential from region 'space' to region 'wire' and the gradient information with jump goes in the opposite direction.)

label wire_interface = wire_mesh.boundaryMesh().findPatchID("wire_space") ;
label space_interface = space_mesh.boundaryMesh().findPatchID("wire_space" );

// create a reference to the gradient in the wire_space patch of region space
fixedGradientFvPatchScalarField& space_T_patch =
refCast<fixedgradientfvpatchscalarfield>(space_T.b oundaryField()[space_interface ]);
scalarField& space_T_patch_gradient = space_T_patch.gradient();

// calculate vector field of unit normals
surfaceVectorField unitNormalField = wire_T.mesh().Sf()/wire_T.mesh().magSf();

// calculate condition for the grad of space_T on boundary wire_space
// require that grad of space_T is equal to match
volVectorField grad_wire_T = fvc::grad(wire_T);
volVectorField match = -(grad_wire_T - (jumpFactor*R_w.value()*wire_M_w));
scalarField space_T_grad_condition =
(match.boundaryField()[wire_interface] & unitNormalField.boundaryField()[wire_interface]);

// use grad of potential condition (wire -> space)
space_T_patch_gradient = space_T_grad_condition;

// use potential condition (space -> wire)
wire_T.boundaryField()[wire_interface] == space_T.boundaryField()[space_interface];


When I run this solver on a simple domain --- a 1D problem, in which the boundary between the two regions is parallel to the x-axis --- the calculated solution is as expected. However when I apply it to my 2D problem, the solution grows without apparent bound - nearly doubling at every iteration.

Any ideas would be appreciated,
KM

Hi K.! I didn't go through
 

Hi K.!

I didn't go through it in detail (and I don't know where exactly in your code this should be corrected), but is it possible, that the other diffeence between your two grids is that in the first case the grid-spacing on the left and the right of the boundary is the same while in the second case they differ by a factor of 2?

Have you checked the obundary conditons for T for the non-interface boundaries?

Hi Bernhard, Thanks for the
 

Hi Bernhard,

Thanks for the reply.

I was worried about the boundary between the two meshes, i.e., whether the faces of the corresponding cells on the commom boundary of both meshes were aligned, and so I wrote code to:

• Determine the two unit normals fields, using .Sf(), and compare to ensure that the normal for corresponding faces were pointing to opposite directions.
• Determine the face centres, using .Cf(), and compare to ensure that the face centres for corresponding faces were equal.

However, the cell volumes (or just areas since it is a 2D problem) may well be different on either side of the common boundary of the two meshes. In the 1D case I ran a few experiments with different size cells on either side of the boundary (while keeping the boundary faces aligned). Varying the cell volume on either side of the boundary had no effect so I assumed that this was not the source of the problem.
(Bad idea?)

If the relative cell volumes did have an effect how do I add a correction for this when transferring the gradient information between meshes (I think that the transfer of potential information is OK - it is only one line!).

The boundaries conditions for the non-interface are that T is zero.

KM

OK. The mesh stuff was only a
 

OK. The mesh stuff was only a guess. Sorry for asking.

About the BC: it's zero, not zeroGradient? (sorry for asking, but I want to be sure) Does the solution get unbounded in the whole domain or only near your interface? Can you provide us a picture?

are you sure about the code pa
 

are you sure about the code part that compute the normals?

Can you explain better the jump condition? and your BC?
In my experience such a problem can appear if you wrong the sign of the heat flux (scalar gradient).
Are the heat fluxes pointing in the same direction at the boundary (going in or out from the circle)?

Daniele