CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Setting normal component of a vector to zero

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By ARTem

Reply
 
LinkBack Thread Tools Display Modes
Old   July 19, 2013, 15:48
Default Setting normal component of a vector to zero
  #1
Member
 
Chris
Join Date: Aug 2012
Location: Calgary, Alberta, Canada
Posts: 73
Rep Power: 4
ChrisA is on a distinguished road
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:

Code:
               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.
ChrisA is offline   Reply With Quote

Old   July 22, 2013, 03:23
Default
  #2
Member
 
Artem Shaklein
Join Date: Feb 2010
Location: Russia, Izhevsk
Posts: 43
Rep Power: 7
ARTem is on a distinguished road
Why don't you consider setting up diffusive flux of your variable T to zero by using zeroGradient (dT/dn = 0) boundary condition?
ARTem is offline   Reply With Quote

Old   July 22, 2013, 03:29
Default
  #3
Member
 
Chris
Join Date: Aug 2012
Location: Calgary, Alberta, Canada
Posts: 73
Rep Power: 4
ChrisA is on a distinguished road
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).
ChrisA is offline   Reply With Quote

Old   July 22, 2013, 05:36
Default
  #4
Member
 
Artem Shaklein
Join Date: Feb 2010
Location: Russia, Izhevsk
Posts: 43
Rep Power: 7
ARTem is on a distinguished road
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

Code:
    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
Code:
        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 likes this.
ARTem is offline   Reply With Quote

Old   July 22, 2013, 09:32
Question one additional question...
  #5
New Member
 
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 29
Rep Power: 7
andyru is on a distinguished road
@ARTem:

Many thanks for the explanation.

One question. After:
Code:
...
forAll(mesh.boundary()[patchi], facei)
        {
            igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];
        }
}
which creates the boundary diffusive flux,
the next function call is:
Code:
gGrad.correctBoundaryConditions();
with:

Code:
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
andyru is offline   Reply With Quote

Old   July 22, 2013, 12:21
Default
  #6
Member
 
Chris
Join Date: Aug 2012
Location: Calgary, Alberta, Canada
Posts: 73
Rep Power: 4
ChrisA is on a distinguished road
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.

Last edited by ChrisA; July 22, 2013 at 15:15.
ChrisA is offline   Reply With Quote

Old   July 23, 2013, 03:03
Default
  #7
Member
 
Artem Shaklein
Join Date: Feb 2010
Location: Russia, Izhevsk
Posts: 43
Rep Power: 7
ARTem is on a distinguished road
Hello, andyru.
Don't forget about functions with same name but taking different arguments.
Code:
gGrad.correctBoundaryConditions();
takes no arguments
and
Code:
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.

Quote:
Originally Posted by ChrisA View Post
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.
ARTem is offline   Reply With Quote

Old   July 23, 2013, 12:36
Default
  #8
Member
 
Chris
Join Date: Aug 2012
Location: Calgary, Alberta, Canada
Posts: 73
Rep Power: 4
ChrisA is on a distinguished road
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.
ChrisA is offline   Reply With Quote

Old   July 24, 2013, 03:01
Default
  #9
Member
 
Artem Shaklein
Join Date: Feb 2010
Location: Russia, Izhevsk
Posts: 43
Rep Power: 7
ARTem is on a distinguished road
Quote:
Originally Posted by ChrisA View Post
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
Code:
surfaceScalarField DTf("DTf", fvc::interpolate(DT));
DTf.write();
Check out different interpolation schemes. I'm sure, you will get "1" at all internal faces.
ARTem is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Setting Opening Boundary to Move Normal to Itself NCle CFX 10 April 20, 2012 07:54
Cells with t below lower limit Purushothama CD-adapco 2 May 31, 2010 21:58
my stl surface is seen as just a line rcastilla OpenFOAM Meshing & Mesh Conversion 2 January 6, 2010 02:30
NACA0012 geometry/design software needed Franny Main CFD Forum 13 July 7, 2007 15:57
Warning 097- AB CD-adapco 6 November 15, 2004 05:41


All times are GMT -4. The time now is 01:27.