# Calculation Tangential & Normal gradients of a Scalar near a surface.

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

 September 9, 2016, 00:50 Calculation Tangential & Normal gradients of a Scalar near a surface. #1 New Member     Tony Join Date: May 2016 Posts: 22 Rep Power: 10 Hi Everyone (First post ), I was wondering if I could get some advice on a how to approach a code snippet. (I am only just starting to program with OpenFOAM, please forgive any major oversights, please feel free to point them out though.) For my problem, I am interested in the tangential and normal components of Electric field to the surface of a patch. Electric field strength can be found the through the following calculation: I am currently having the following difficulties:I believe I would like to utilize boundaryField() as I only need the values for the cells in contact with the patch. However, this seems to be unavailable for volScalarFields which I need perform the calculations. (examples like wallGradU.C and the linked posts based the calculations off of volVectorFields). Due to this, I am currently unable to utilize both boundaryField and snGrad to target the wall and compute the tangential gradient. Considering the post:http://www.cfd-online.com/Forums/ope...tml#post186010, is snGrad appropriate? Using an example comparable to the post, assuming the wall has a no-slip condition applied, the change in the velocity across the surface should be zero resulting in no gradient. In my case, if I were to apply this function to a conductor (equipotential surface) I would expect . Finding across boundary cells would be nice though to observe the rate of change in Etan. One step at a time though. Links for further investigation of these boundary cells would be greatly appreciated. For the normal component in my snippet, would I need to interpolate the E to get the value at the surface or is this already performed via boundaryField. (http://www.cfd-online.com/Forums/ope...hi-patchi.html) Since I am calculating values at the surface, am I considering the correct type? (Should I be considering surfaceVectorField?) Considering the previous post and wallGradU.C, it appears I should keep it as a volVectorField. My code snippet: (Note: this snippet does not build unless Ue is a volVectorField.) Code: //Ue is the volScalarField potential determined in the solver. volVectorField E ( IOobject ( "E", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), -fvc::grad(Ue) ); volVectorField Enorm ( IOobject ( "Enorm", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh, dimensionedVector ( "Enorm", E.dimensions(), vector::zero ) ); volVectorField Etan ( IOobject ( "Etan", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh, dimensionedVector ( "Etan", E.dimensions(), vector::zero ) ); //Create references for Calculation forAll(Etan.boundaryField(), patchi) { Etan.boundaryField()[patchi] = Ue.boundaryField()[patchi].snGrad(); } forAll(Enorm.boundaryField(), patchi) { Enorm.boundaryField()[patchi] = E.boundaryField()[patchi]-Etan.boundaryField()[patchi]; } Note: For this calculation, I have been using the following application for a basis: Code: $FOAM_APP/utilities/postProcessing/wall/wallGradU/wallGradU.C I have been looking into the following threads for guidance: http://www.cfd-online.com/Forums/ope...tml#post186010 http://www.cfd-online.com/Forums/ope...hi-patchi.html  September 11, 2016, 12:24 #2 Senior Member Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15 Hey Tony, Since you haven't received any response yet, I wanted to introduce you this post. http://www.cfd-online.com/Forums/ope...tml#post254946 Just bear in mind that the opposite cell center should help you to create the normal component. Tangential component will be derived from the subtraction of and However, one thing does not match in your efforts. a volScalarField does not have a tangential or normal component.  September 11, 2016, 16:41 #3 New Member Tony Join Date: May 2016 Posts: 22 Rep Power: 10 Hi Babak, I have looked into your referenced post. Unfortunately, to perform that approach it looks as if the normals are determined with respect to the location of cell face, center, and opposing face. I should have mentioned that I will be working with a mesh that is not necessarily regular (I will have some oblique cells at some patches). Attached is my current snippet: The main changes are that I changed the type of Etan and Enorm to surfaceVectorFields and I am attempting Code:  forAll(Enorm.boundaryField(), patchi) { Enorm.boundaryField()[patchi] = E.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi] /mesh.magSf().boundaryField()[patchi]; } // forAll(Etan.boundaryField(), patchi) // { // Etan.boundaryField()[patchi] = // E.boundaryField()[patchi]-Enorm.boundaryField()[patchi]; // } In this method I am trying to get the E at the boundary, then perform an inner product with constructed normal vector of the same surface. (This idea is from$FOAM_APPS//utilities/postProcessing/wall/wallShearStress/wallShearStress.C) Unfortunately, I still have having issues building this. Generally my errors consist of "is not a member of" or lack of initialization. There first error for this compile is : Code: ElectricComponents.H:55:47: error: no match for ‘operator=’ (operand types are ‘Foam::fvsPatchField >’ and ‘Foam::tmp >’) Enorm.boundaryField()[patchi] = Note: I have also edited Etan and Enorm to surfaceVectorFields : Code:  surfaceVectorField Enorm ( IOobject ( "Enorm", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh, dimensionedVector ( "Enorm", E.dimensions(), vector::zero ) ); surfaceVectorField Etan ( IOobject ( "Etan", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh, dimensionedVector ( "Etan", E.dimensions(), vector::zero ) );

 September 11, 2016, 17:37 #4 Senior Member   Bobby Join Date: Oct 2012 Location: Michigan Posts: 454 Rep Power: 15 Hey Tony what I know is that: within a patch, the face normals are given by: Code:  patch.nf() and for looping over different patches (or specific patch ID) the following snippet gives face normals: Code:  mesh.Sf().boundaryField()[patchID] / mesh.magSf().boundaryField()[patchID] Now, in your case: Code: volVectorField E= -fvc::grad(Ue) forAll (E.boundaryField() , patchi) { const vectorField& surfaceNormal = mesh.boundary[patchi].nf() forAll(E.boundaryField() [patchi], celli) { scalar Enorm= surfaceNormal [celli] & E[celli] } } Try the above snippet without making it super-complicated

 September 11, 2016, 18:20 #5 New Member     Tony Join Date: May 2016 Posts: 22 Rep Power: 10 Hi Babak, I didn't mean for it to become complicated. I am just trying to make the documentation of the problem easier so it can be referenced easier later. On a more complicating note, I realized that I shouldn't be doing an inner product as my objective is a vector. I would be looking to do element-wise multiplication , where is the normal vector. That is unless the is a simpler way to get the normal component of with respect to a boundary. Concerning your suggestion, I inserted your snippet: Code: forAll (E.boundaryField() , patchi) { const vectorField& surfaceNormal = mesh.boundary[patchi].nf(); forAll(E.boundaryField() [patchi], celli) { scalar Enorm= surfaceNormal[celli] & E[celli]; } } and I received the following error: Code: ElectricComponents.H: In function ‘int main(int, char**)’: ElectricComponents.H:55:59: error: invalid types ‘[Foam::label {aka int}]’ for array subscript const vectorField& surfaceNormal = mesh.boundary[patchi].nf(); ^ ElectricComponents.H:58:13: warning: unused variable ‘Enorm’ [-Wunused-variable] scalar Enorm= surfaceNormal[celli] & E[celli]; ^ I tried changing the type from vectorField to volVectorField and surfaceVectorField and produced the same error.

September 11, 2016, 23:09
#6
New Member

Tony
Join Date: May 2016
Posts: 22
Rep Power: 10
I think I am getting closer.

I am working with the following code snippet:
Code:
	const surfaceVectorField& faceNormal = mesh.Sf()/mesh.magSf();
forAll (E.boundaryField() , patchi)
{
Enorm.boundaryField() [patchi] = faceNormal.boundaryField() [patchi];
}
This idea came from:
Code:
\$FOAM_SRC/fvMotionSolver/motionDiffusivity/directional/directionalDiffusivity.C
The code successfully builds with however I get an error when running a case.
I chose to verify with a simple case. I've attached a .png of the simple grid (a rhombohedron)

I changed the E to be a constant vector <1,0,0> and all other settings for the case are zero, zeroGradient, or empty.
Running the case, I get the following error:
Code:
--> FOAM FATAL ERROR:
different patches for fvsPatchField<Type>s

From function PatchField<Type>::check(const fvsPatchField<Type>&)
in file /home/anthony/OpenFOAM/OpenFOAM-2.1.1/src/finiteVolume/lnInclude/fvsPatchField.C at line 155.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  Foam::fvsPatchField<Foam::Vector<double> >::operator=(Foam::fvsPatchField<Foam::Vector<double> > const&) at ??:?
#3
at ??:?
#4  __libc_start_main in "/lib64/libc.so.6"
#5
at /home/abuild/rpmbuild/BUILD/glibc-2.19/csu/../sysdeps/x86_64/start.S:125
Anyone have any ideas?
Attached Images
 SimpleMesh.PNG (24.5 KB, 18 views)

 September 12, 2016, 09:27 #7 New Member     Tony Join Date: May 2016 Posts: 22 Rep Power: 10 This solution almost works. The issue is that it interpolates the values of E to the wall. In the case of a wall at equipotential, this fails as Etan is non-zero. (Etan should be zero at an equipotential surface as all values along the surface are the same. ) Code:  volVectorField E ( IOobject ( "E", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), -fvc::grad(Ue) ); surfaceVectorField Enorm ( IOobject ( "Enorm", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), cmptMultiply(fvc::interpolate(E), mesh.Sf()/mesh.magSf()) ); surfaceVectorField Etan ( IOobject ( "Etan", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), fvc::interpolate(E)-Enorm ); forAll (Enorm.internalField() , patchi) { Enorm.internalField() [patchi] = vector::zero; } forAll(Etan.internalField(), patchi) { Etan.internalField()[patchi] = vector::zero; } To clarify, I am trying to obtain the scalar values at the surface, then perform a gradient with respect to the surface norm. Note: Learned about cmptMultiply for element-wise vector multiplication, Nice. Note: I found a thread looking exactly for Etan. Unfortunately, the thread is unresolved since 2013. http://www.cfd-online.com/Forums/ope...d-surface.html Last edited by AnthonyP; September 12, 2016 at 13:50. Reason: Additional thread found.

 October 2, 2016, 14:40 #8 New Member     Tony Join Date: May 2016 Posts: 22 Rep Power: 10 Hi Everyone, I just wanted to share my current snippet. They still interpolate the volVector to the wall, but I am still trying to find a way to perform the gradient directly at the wall (ex. gradient performed on surfaceScalar field or pointField) Code: //Ue is the volScalarField potential determined in the solver. volVectorField E ( IOobject ( "E", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), -fvc::grad(Ue) ); surfaceVectorField Enorm ( IOobject ( "Enorm", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), (fvc::interpolate(E) & mesh.Sf()/mesh.magSf()) * mesh.Sf()/mesh.magSf() ); surfaceVectorField Etan ( IOobject ( "Etan", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), fvc::interpolate(E)-Enorm ); surfaceVectorField EnormBound ( IOobject ( "EnormBound", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), Enorm ); surfaceVectorField EtanBound ( IOobject ( "EtanBound", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), Etan ); //These zero out the internal field forAll (EnormBound.internalField() , patchi) { EnormBound.internalField() [patchi] = vector::zero; } forAll(EtanBound.internalField(), patchi) { EtanBound.internalField()[patchi] = vector::zero; } babakflame likes this.

October 6, 2016, 12:57
#9
New Member

Tony
Join Date: May 2016
Posts: 22
Rep Power: 10
Unfortunately, this code snippet is still very sensitive to the cell orthogonality.
for example if and

I have attached the error of E magnitude (which is also Etan in this case) for a geometry composed of orthogonal cells and one with non-orthogonal cells. Both have uniform spacing (simple rectangle and rhombus). The left side of each shape is positioned at x=0 or very close (1e-15). Note that the error is very high for the first cells to which are near the discontinuity; but in the rectangle case, it quickly drops to 0% error; while rhombus maintains a constant 12% error.

I am still trying to figure openFOAM out and I am open to any hints to remedy this orthogonality issue. If I find a solution in the meantime, I will post it.

*Note: if I am posting this in the wrong sub-forum please let me know.
Attached Images
 ErrorRect_E1overX.PNG (16.2 KB, 30 views) ErrorRhom_E1overX.PNG (20.0 KB, 20 views)

 January 17, 2017, 02:31 #10 New Member     Tony Join Date: May 2016 Posts: 22 Rep Power: 10 Hi Everyone, I wanted to post an update with some fresh snippets. I believe I am really close to a solution, but it is still evading me. (not matching with validation). Values in my case are showing to be too large. I am using laplacianFoam as a base solver and running to steady state. The geometry is a wedge in prolate spheroidal coordinates. Top scalar is 1 at xi=0.9915 and bottom face is xi=0.98706 and T(0.98706) from the paper. validation: validation.PNG source: http://ideaexchange.uakron.edu/cgi/viewcontent.cgi?article=1084&context=polymer_ideas (eqs 9 and ~13) case: E.PNGetan.PNGen.PNG I am using the following code snippet to perform my tangent/normal calculation. Code: volVectorField refVector ( IOobject ( "refVector", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), mesh //mesh.C() //temporary space for unit space vectors ); volScalarField distToWall ( IOobject ( "distToWall", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), refVector.component(vector::X)//dummy list to overwrite ); volVectorField E ( IOobject ( "E", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), -fvc::grad(T) //stole from laplacian solver ); Find the position with respect to surface (I want the cell to know the closest patch face and its normal) Code: const word electrode = "electrode"; const word body = "nozzleBody";//counterPotential const polyBoundaryMesh& patches = mesh.boundaryMesh(); const volVectorField& cells = mesh.C(); // get the position for the particle vector position = cells[1]; // dummy vectors vector surfNorm = vector::zero;//initialize a zero vector scalar wallDist = 1000000.0; // default an ungodly high distance forAll(cells,k) { position = cells[k]; surfNorm = vector::zero; wallDist = 1000000.0; //loop over all patches forAll(patches, patchI) { if ( patches[patchI].name() == electrode || patches[patchI].name() == body )//loop through only the desired walls { polyPatch patch = patches[patchI];//get the patch id vectorField n = patch.faceNormals();//get the face normal vectors vectorField C = patch.faceCentres();//get the face center forAll(C, i)// for a location (k) loop through the boundary face centers { vector normal = n[i]/mag(n[i]);//get the unit normal vector of the face scalar di = mag( position - C[i] );//calculate the distance if (di < wallDist) { wallDist = di;//save smallest distance (as it loops you obtain the smallest distance for the given boundary if(patches[patchI].name() == electrode) { surfNorm = -1*normal;//save the normal -negative applied here due to my case (orientation of cell norm of "electrode" was opposite to E) } else { surfNorm = normal; }; }; }; }; }; refVector[k]=surfNorm; distToWall[k]=wallDist; }; my current code has an issue writing at the patch face. (it instead retains intial values) below is a failed attempt to fix that. Code: //attempt (not working) get the write on patch faces (otherwise it retains initialization) forAll(refVector.boundaryField(), patchi) { refVector.boundaryField()[patchi]= mesh.Sf()/mesh.magSf(); } Finally my calculation for Enorm and Etan. (note I do not like the dimenson issue I am havng to go around to calculate Etan. Code:  volVectorField EnormVol ( IOobject ( "EnormVol", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), (E & refVector) * refVector ); const dimensionedScalar changeDim = dimensionedScalar("changeDim", dimensionSet(0, -2, 0, 0, 0, 0, 0), scalar(1.0)); //dimension issue - Perhaps a source of error? volVectorField EtanVol ( IOobject ( "EtanVol", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), E - EnormVol*changeDim ); This snippet was created with reference to the following thread: Find a particle's distance to a wall I hope this helps and I still look forward to assistance. Tony Last edited by AnthonyP; January 17, 2017 at 16:16. Reason: I was told to shorten my snippets if I want help...

 January 18, 2017, 09:14 #11 New Member     Tony Join Date: May 2016 Posts: 22 Rep Power: 10 Grid refinement studies improved the accuracy the E and En, however Etan is far off. What I am trying to do is simply E-En=Et. If it helps, Etan is approximately 25x larger than it should be. Code:  volVectorField EnormVol ( IOobject ( "EnormVol", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), (E & refVector) * refVector ); volVectorField EtanVol ( IOobject ( "EtanVol", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), E - EnormVol ); *Note: refVector is the stored surface unit vector of the closest face on the chosen wall patch. Dimension issues resolved: RefVector had a unit assignment in the 0/ folder (stupid mistake) Think anyone can figure out why I have a dimension change when calculating Etan? I think this is the remaining error. I had to do the dimension change for compilation, is there an area I must divide by to proper perform my calculation? Last edited by AnthonyP; January 18, 2017 at 19:07. Reason: Revised snippet (found dimension issue, but still Etan is off)

 January 18, 2017, 19:04 #12 New Member     Tony Join Date: May 2016 Posts: 22 Rep Power: 10 I am finding something odd... For some reason the calculation in the IOobject E-Enorm produces off results Etan.PNG Despite E and Enorm being valid. E.PNG and En.jpg If I try to do this E-Enorm in Paraview using the calculator filter i get the same result. BUT if I export from paraview to csv. Then do E-Enorm the result is fine. (Am I missing something here?) Etan_post.PNG I am using OF2.1.1 and Paraview 5.1.8 (also used 5.2.0)

 July 22, 2020, 09:01 #13 Member   X Join Date: Jan 2019 Posts: 63 Rep Power: 7 Hey Anthony, Thank you for the post! I am working on electrohydrodynamic and wanted to find tangential and normal component of E at the fluid interface. I know you indicated that you were able to get correct E and Enorm, but were you able to get correct Etan?