[ChrisA]

In my efforts to implement a thermal diffusion model I've run into another snag. Essentially what I think is happening is I have species diffusing through the wall. I suspect this is because I was severely underpredicting the wall concentration when I was giving no special treatment to the diffusion coefficient at the wall.

Initially I tried just setting the diffusion coefficient at the wall boundary to zero, but that seems to lead to an over prediction of the wall concentration. I suspect it isn't working because that doesn't allow for species to diffuse away from the wall. So what I need to do is check the direction of the species transport and if it is into the wall set the wall normal component to zero, otherwise let it happen... here's what I have so far:

               volVectorField JTherm = -DT*fvc::grad(T);
 	      	 	forAll(JTherm.boundaryField(),patchi)
 	      	 	{
 	      	 	if (isA<wallFvPatch>(mesh.boundary()[patchi]))
 	      	 	{
 	      	 		const vectorField& surfaceNormal = mesh.boundary()[patchi].nf();
 	      	 		forAll(JTherm.boundaryField()[patchi],celli)
					{
						scalar dotprod = surfaceNormal[celli] & JTherm[celli];
						if(dotprod < 0)
						{
//set normal component to zero

Where I run into an issue is I'm having problems setting the normal component of this vector to zero, is there an easy way to do that? (I also might have my signs mixed up as far as the vector of interest goes, I'll check that once I get the core code down).

Edit: well that was a bit silly, two cross products should do exactly this.

[ARTem]

Why don't you consider setting up diffusive flux of your variable T to zero by using zeroGradient (dT/dn = 0) boundary condition?

[ChrisA]

You're right that would work for cases where T (temperature) has a zero gradient at the wall, however this is not true for the case I'm trying to run (fixed wall temperature based on experiments).

[ARTem]

Hello, ChrisA.

Ok, I get your point)

When you computes JTherm = -DT*fvc::grad(T); non-zero diffusive flux is already in cell adjacent to boundary. And JTherm boundary will not be used in further computations (strictly speaking, it will be used, but will make changes only for variable fields with calculated bc types).

Why you are trying to take gradient of vector? Gradient means changing of your variable across each of three dimensions. So, you want to get 3x3 tensor matrix using gradient operation (∂T1/∂x, ∂T1/dy, ∂T1/∂z, ∂T2/∂x....)?

Anyway, if you want to compute grad(T) without wall diffusive flux affection, just dig dipper: fvc::grad(volField) with Gauss scheme is computed in gaussGrad.C (you can find file by "locate gaussGrad.C"). To be more precise, function to search for is gradf.

This part of function operates with boundaries

    forAll(mesh.boundary(), patchi)
    {
        const labelUList& pFaceCells =
            mesh.boundary()[patchi].faceCells();

        const vectorField& pSf = mesh.Sf().boundaryField()[patchi];

        const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];

        forAll(mesh.boundary()[patchi], facei)
        {
            igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];
        }
    }

And this part creates boundary diffusive flux for your variable

        forAll(mesh.boundary()[patchi], facei)
        {
            igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];
        }

If you want to get rid of your problem, you can create another function, e.g. gradfZeroBC, copy in ot all from gradf except boundary treatment, build way up to fvc::gradNoBC in fvcDiv.C and compile library. This will do the job.

[andyru]

@ARTem:

Many thanks for the explanation.

One question. After:

...
forAll(mesh.boundary()[patchi], facei)
        {
            igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];
        }
}

which creates the boundary diffusive flux,
the next function call is:

gGrad.correctBoundaryConditions();

with:

forAll(vsf.boundaryField(), patchi)
    {
        if (!vsf.boundaryField()[patchi].coupled())
        {
            vectorField n =
                vsf.mesh().Sf().boundaryField()[patchi]
               /vsf.mesh().magSf().boundaryField()[patchi];

            gGrad.boundaryField()[patchi] += n *
            (
                vsf.boundaryField()[patchi].snGrad()
              - (n & gGrad.boundaryField()[patchi])
            );
        }
     }

Why do we loop again over the boundaryFields?
For what does this part stand for?

Best regards,

Andy

[ChrisA]

Before I start I'd just like to note that I don't 100% understand how OF handles boundaries so please forgive any ignorance included below.

This is good discussion, probably useful for those who are having issues setting species diffusion to zero at boundaries (mass flow inlets for example) among other things but I think you are missing what I'm trying to do, which is mostly my fault.

I should have been more clear with what I'm trying to implement. I'm trying to use a species thermal diffusion model, which is a species tendency to travel along temperature gradients. This is why you see a grad(T) in my code. T is temperature and its "transport", for lack of a better word, of course gets handled in the energy equation. With the model I end up with a mass flux, JTherm, based on the temperature gradient. The next step in the code is to compute the transport of species Yi by taking the divergence of JTherm. While manipulating the way the gradient is taken at the wall will fix the species transport it will impact the energy transport in a non-physical way (I think?).

Now if I understand the first part of what you said right (the bit about diffusive flux being in the cell adjacent to the boundary). I don't think this is the case because when I change DT at the boundary (example set it to zero or the internal field value) the solution changes. If what you're saying is true wouldn't the solution be independent of DT's boundary value?

I'm also not sure where I'm taking the gradient of a vector? Perhaps inherent in the confusion as to what T is in this case? If that is the case I apologize for not being clearer from the start.

[ARTem]

Hello, andyru.
Don't forget about functions with same name but taking different arguments.

gGrad.correctBoundaryConditions();

takes no arguments
and

template<class Type>
void Foam::fv::gaussGrad<Type>::correctBoundaryConditions
(
    const GeometricField<Type, fvPatchField, volMesh>& vsf,
    GeometricField
    <
        typename outerProduct<vector, Type>::type, fvPatchField, volMesh
    >& gGrad
)

takes 1 argument.

2ChrisA. Sorry for grad of a vector stuff, my mess)
So, you have an equation
∂(ρ Yi)/∂t + div(ρ Yi U) - div( (µ/Sc + µ_t/Sc_t) grad(Yi) ) = div ( DT grad(T) ) ?
Strange, but DT = 0 at walls should do the job.

Originally Posted by ChrisA

Initially I tried just setting the diffusion coefficient at the wall boundary to zero, but that seems to lead to an over prediction of the wall concentration. I suspect it isn't working because that doesn't allow for species to diffuse away from the wall.

But you've set DT = 0 just at walls? Because when div(DT grad(T)) is being computed, DT is interpolated to internal faces and boundary values have no part in it.
You can try to do next (with DT=0 at walls). Consider just 1d case with 2 walls. Firstly, set, e.g., T_left = 400, T_right = 100. When steady-state is achieved, switch values between T_left and T_right and continue a simulation. Compare intermediate and final solutions. Is there a flux from cell near to wall in domain with DT = 0 at walls?

About gradT and energy equation. If you create new set of functions and will call fvc::gradZeroBC(T), with original functions untouched, all will be fine.

[ChrisA]

ARTem,

That is the general idea behind my species transport equation yes. It's known as Soret diffusion and is often neglected due to its relatively high computation cost but it plays a non negligible role in the diffusion of notably different sized molecules when a temperature gradient is present (that's an oversimplification of when its non-negligible but that's the case for what I'm working on).

I see what you mean now with the grad operator.

DT=0 at the walls ended up severely over predicting the wall concentration. When the DT value is interpolated to the internal faces does it use the wall value at all in that interpolation? I ask because I think the over prediction is because at some point things want to start diffusing away from the wall. The equation is in such a form that the sign will flip if the concentration of the diffusing species gets high enough, so what was once diffusing to the wall will diffuse away till steady state and whatnot.

I should also note that I've been suspecting the boundary as the source of my issue as the concentrations away from the wall are lining up quite well with simulations in previous work on the same problem.

[ARTem]

Originally Posted by ChrisA

ARTem,
DT=0 at the walls ended up severely over predicting the wall concentration. When the DT value is interpolated to the internal faces does it use the wall value at all in that interpolation?

It's very strange for me. You can check this out for youself:
create 1d mesh with 2 nodes (or 3d with 8 nodes - 2 in each directions - to be sure by 100%). Create volField DT with, e.g. value "1", and in one boundary (or in all boundaries) set DT to be "0". Then call for

surfaceScalarField DTf("DTf", fvc::interpolate(DT));
DTf.write();

Check out different interpolation schemes. I'm sure, you will get "1" at all internal faces.