|
[Sponsors] |
Calculation Tangential & Normal gradients of a Scalar near a surface. |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
September 9, 2016, 01: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:
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]; } For this calculation, I have been using the following application for a basis: Code:
$FOAM_APP/utilities/postProcessing/wall/wallGradU/wallGradU.C http://www.cfd-online.com/Forums/ope...tml#post186010 http://www.cfd-online.com/Forums/ope...hi-patchi.html |
|
September 11, 2016, 13:24 |
|
#2 |
Senior Member
Bobby
Join Date: Oct 2012
Location: Michigan
Posts: 454
Rep Power: 16 |
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, 17: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]; // } (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<Foam::Vector<double> >’ and ‘Foam::tmp<Foam::Field<double> >’) Enorm.boundaryField()[patchi] = 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, 18:37 |
|
#4 |
Senior Member
Bobby
Join Date: Oct 2012
Location: Michigan
Posts: 454
Rep Power: 16 |
Hey Tony
what I know is that: within a patch, the face normals are given by: Code:
patch.nf() Code:
mesh.Sf().boundaryField()[patchID] / mesh.magSf().boundaryField()[patchID] 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] } } |
|
September 11, 2016, 19: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]; } } Code:
ElectricComponents.H: In function ‘int main(int, char**)’: ElectricComponents.H:55:59: error: invalid types ‘<unresolved overloaded function type>[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]; ^ |
|
September 12, 2016, 00: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]; } Code:
$FOAM_SRC/fvMotionSolver/motionDiffusivity/directional/directionalDiffusivity.C 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 |
|
September 12, 2016, 10: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; } 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 14:50. Reason: Additional thread found. |
|
October 2, 2016, 15: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; } |
|
October 6, 2016, 13: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. |
|
January 17, 2017, 03: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 ); 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; }; 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(); } 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 ); 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 17:16. Reason: I was told to shorten my snippets if I want help... |
|
January 18, 2017, 10: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 ); 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 20:07. Reason: Revised snippet (found dimension issue, but still Etan is off) |
|
January 18, 2017, 20: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, 10: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? |
|
Tags |
calculation, gradient, patch surface |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
How is the normal unit vector calculated for 2D airfoil surface gradients? | Ry10 | SU2 | 5 | February 27, 2018 07:28 |
surface tension as function of normal gradient | tom | FLUENT | 2 | September 7, 2016 04:48 |
dieselFoam problem!! trying to introduce a new heat transfer model | vivek070176 | OpenFOAM Programming & Development | 10 | December 24, 2014 00:48 |
[snappyHexMesh] Layers don't fully surround surface | EVBUCF | OpenFOAM Meshing & Mesh Conversion | 14 | August 20, 2012 05:31 |
[CAD formats] my stl surface is seen as just a line | rcastilla | OpenFOAM Meshing & Mesh Conversion | 2 | January 6, 2010 02:30 |