# Calculating gradients & discretization error

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

 December 6, 2013, 08:10 Calculating gradients & discretization error #1 Senior Member   Join Date: Oct 2013 Posts: 379 Rep Power: 7 Hello Foamers, I'm having an issue with an error coming from the discretization & gradient calculation. I have calculated an electric potential (using laplace(sigma, phi)=0), from which I then calculate the gradient to get the electric field (using -fvc::grad(phi)). The field is smeared out, likely due to central differencing being used to calculate the gradient. This is probably unavoidable and doesn't pose such a huge problem. The issue comes up when I want to calculate the electric current from the field and the conductivity volScalarField sigma. The conductivity is not smeared out.This leads to a current field which is not divergence-free as it should be. My question now is how can I avoid this problem and get a divergence-free current field? Is it possible to smooth the conductivity field to correct this error? By the way, and just for reference for others, this problem also comes up at the boundaries, because the conductivity set at the boundaries differs from the conductivity in the field. ZeroGradient BC didn't work for me, it didn't produce any proper solution. Any idea why? However, this issue can be fixed by explicitly setting the boundary field values in the solver to the values of the neighbouring cells. Example code can be found in this post: http://www.cfd-online.com/Forums/ope...tml#post407396

 December 6, 2013, 10:49 #2 Senior Member     Philipp Join Date: Jun 2011 Location: Germany Posts: 1,297 Rep Power: 19 The current isn't divergence free, because the differential euqation you solve for isn't a current-conserving equation. It is just the laplace equation. What you need to do to get a current coserving equation is solving div(J/sigma) = 0, with J = sigma*E (with J the current density vector and E the electric field vector). If you solve this equation with finite volumes in conservative form the result will be divergence free, also the current density. Unfortunately you will need a boundary condition for E. I guess you want to set the voltage "phi" on the electrodes? I think this doesn't work and you need to solve laplace equation as you already did. 1) What do you mean by "smeared out"? 2) Where do you calculate the electric field? At the cell centers or at the cell faces? 3) Why / Where did you try ZeroGradient BC? For the potential? __________________ The skeleton ran out of shampoo in the shower.

December 9, 2013, 04:44
#4
Senior Member

Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 19
Quote:
 Originally Posted by chriss85 I'm not sure that div(J/sigma)=0 is correct, because this would effectively mean that J/sigma (i.e. E) is divergence-free, and not J?
Chris, but this is Maxwell's equation (at least if the charge density is zero).
Edit: Generally it should be div(J/sigma) = rho / epsilon.

But: Do you solve the euquations in regions with different / varying conductivity? If so, you will get charge density or surface charges. You can not just solve d2 Phi / dx^2 = 0 then.
__________________
The skeleton ran out of shampoo in the shower.

 December 9, 2013, 05:02 #5 Senior Member   Join Date: Oct 2013 Posts: 379 Rep Power: 7 I'm fine with div(J/sigma) = rho / epsilon of course. I'm solving a part of the charge conservation equation in the steady limit ( -> drho/dt = 0): drho/dt + div(j) = 0 div(sigma*E) = 0 div(-sigma*grad(phi)) = 0 laplacian(sigma, phi) = 0 I do have varying conductivity, this is why the conductivity is a field in my application. I believe that the last equation above should produce correct results, and it appears to do so in simple cases (except at the interfaces between different conductivities due to this numerical problem I described above), but please correct me if this is wrong. Do I need to include time-dependence and consider charge distribution? If so, which equations and BCs are suited for this? I've seen that the rodFoam solver uses a similar approach to mine.

 December 9, 2013, 05:26 #6 Senior Member     Philipp Join Date: Jun 2011 Location: Germany Posts: 1,297 Rep Power: 19 Ok, sorry, you are right. I didn't read you initial post properly. I thought you were trying to solve Maxwell's equation instead of charge conservation. Everything is fine, no dt. This "div(J)" problem, where do you see it? Only at the boundaries or also inside the volume? __________________ The skeleton ran out of shampoo in the shower.

 December 9, 2013, 05:33 #7 Senior Member   Join Date: Oct 2013 Posts: 379 Rep Power: 7 It happens at every two cells where there is a step in conductivity. I need to try to get a picture uploaded, I have drawn some graphs that show this problem quite well I think. If you can pm me your mail address I can send it to you, otherwise I'll have to see how to get it online

 December 9, 2013, 05:51 #8 Senior Member     Philipp Join Date: Jun 2011 Location: Germany Posts: 1,297 Rep Power: 19 Ok, sure - at every step in sigma. I think I understand the problem now. What your solver (laplacian(sigma,phi)=0) does, is to solve the charge conservation for each cell. That means, every computational cell is charge conserving by means of what flows into it through the cell's faces. This should work pretty fine in your case. What you do afterwards is calculating the current divergence in the cell centers. This introduces new numerical errors, since up to then, only the divergence at the cell faces was calculated. But the current divergence in the center of each cell is some new artifical value the solver didn't try to reduce. I think you can not avoid this problem. But - as I understand it - this error "isn't really there". It just comes up during evaluation of your results. To get what you actually want, I think you need to sum up the current density of each cell's faces. But this should be similar to the local residuum, right? Do you understand what I am trying to say? __________________ The skeleton ran out of shampoo in the shower.

 December 9, 2013, 05:51 #9 Senior Member   Join Date: Oct 2013 Posts: 379 Rep Power: 7 Here you go: https://drive.google.com/file/d/0B2I...it?usp=sharing I'm not explicitly calculating current divergence, but current density varies where it shouldn't. I'm using a simple "cable" conductor case, a rectangular long cable with two different conductivities. Please see the drawn graphs as illustration. I have thought about looking at the cell faces before, maybe similar to how it is done in the electrostaticFoam solver. I'm just not sure yet how to properly work with surfaceScalarFields when it comes to displaying them in paraView. I also need to get the net currents through the electrodes, right now I'm using patchIntegrate for this. I believe it should also work with surfaceScalarFields, but I need to try this.

 December 10, 2013, 04:19 #10 Senior Member   Join Date: Oct 2013 Posts: 379 Rep Power: 7 I tried using surfaceScalarFields now instead of volVectorFields for electric field and current density, using EField = -fvc::snGrad(phi) and j = fvc::interpolate(sigma) * EField * mesh.magSf(). I'm not quite sure how to properly visualize surfaceScalarFields, so I used fvc::reconstruct after the calculation to get back volVectorFields. The electric field looks very similar if not same to the field I had before, the current density looks a bit different but still wrong. Do you think the calculations are more exact this way? In the end I'm mostly interested in a correct calculation and in the current through the electrodes which I can get with patchIntegrate. If I lose some precision during visualization then so be it, but I'd like to have correct current fields. The problem which arises here is that I have a joule heating term (sigma * EČ) which needs to be reconstructed on the cell centers, so I think I'd still have the error I'm currently seeing?

 December 10, 2013, 04:54 #11 Senior Member     Philipp Join Date: Jun 2011 Location: Germany Posts: 1,297 Rep Power: 19 I don't know much about these different fields and visualization, sorry. I think you deal with the "normal" discretization error. Since you are aware of it and know where it comes from, I would just do it like you did. If someone asks you about that, you will still "look good" (? how do you say that in english ?), because you can explain it and you can also explain why you can't avoid it. 1) Did you try if the current density wobble decreases when the mesh becomes finer? 2) I don't understand, why the current density still looks wrong. For the "fvc::snGrad(phi)" and "fvc::interpolate(sigma)", do you use the same interpolation method as during the solution of your differential equation? __________________ The skeleton ran out of shampoo in the shower.

 December 10, 2013, 05:18 #12 Senior Member   Join Date: Oct 2013 Posts: 379 Rep Power: 7 It's clear that this is due to discretization, I'm just wondering if it actually produces problems and if I can avoid it somehow. 1) I didn't try yet, but I wouldn't expect it to, atleast not directly. I believe that the gradients would result in the same values on finer meshes, so the smeared out step in the electric field would become a bit sharper, but the smeared out values would still be there, same for the current. However, on finer meshes the conductivity will have smaller steps and this should indeed result in a smaller error. Not sure if it sums up again over the larger amount of cells though... 2) I've been using linear interpolation everywhere as far as I can tell. However, different gradient methods for grad and snGrad were used (leastSquares and corrected). I'm not sure if this is the source of the error, I can probably try later on, right now I've reversed my changes again

 December 10, 2013, 05:24 #13 Senior Member     Philipp Join Date: Jun 2011 Location: Germany Posts: 1,297 Rep Power: 19 to 2) Can you try "Gauss linear" and "linear"? __________________ The skeleton ran out of shampoo in the shower.

 December 10, 2013, 05:29 #14 Senior Member     Philipp Join Date: Jun 2011 Location: Germany Posts: 1,297 Rep Power: 19 But wait... even then: The gradient(phi) for the differential equation will be calculated by using (in 1D) the two adjacent cell values. The gradient after solving (snGrad) will interpolate the two cell values of the gradient to the surface. You will still get the wobble. Is there any way to calculate the gradients exactly the same way? __________________ The skeleton ran out of shampoo in the shower.

 December 10, 2013, 06:31 #15 Senior Member   Join Date: Oct 2013 Posts: 379 Rep Power: 7 Not sure, I'm not that experienced with the different algorithms yet. However, I noticed that this seems to produce quite a large error in my application, in a real case I get a difference of about 10% in electrode currents. This is a significant problem...

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post igo OpenFOAM 1 November 20, 2012 21:16 Yogini Fluent UDF and Scheme Programming 7 October 3, 2012 07:24 vivek070176 OpenFOAM Installation 22 June 1, 2010 12:34 phsieh2005 OpenFOAM Running, Solving & CFD 2 March 21, 2009 06:18 cfduser CFX 0 April 29, 2006 10:58

All times are GMT -4. The time now is 21:03.