[chriss85]

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

[RodriguezFatz]

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?

[chriss85]

I've made some progress myself, but I need to look into this further.
Indeed I want to set the potentials at the electrodes. Alternatively I could also use a current driven method, but this would require a more elaborate treatment of the boundaries, it would require an area with high conductivity (representing the electrodes), so the current distribution can develop properly.
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? solving laplace(sigma, phi)=0 should be equal to solving div(j)=0 because of j=-sigma*grad(phi). This is why I would expect that the current field should be divergence free, atleast from an analytical perspective. I'm relatively inexperienced with numerical methods, is there any reason this is not conservative?

a) With smeared out I mean a discretization error that happens at jumps in conductivities between two cells. Unfortunately I can't upload an image right now due to company policies . However, if you assume that the conductivity is a step function, you get a potential with two different gradients. So far, so good. However, when calculating the electric field with -fvc:.grad(phi), you don't reproduce the step function but a "smeared out" step, with two values at the step being wrong. This is what I believe is due to central differencing gradient method. I'm not sure if there is actually a better method there. The problem comes up when I now calculate the current with j=sigma * E. E is "smeared out", while sigma isn't, so the multiplication produces inconstant values at the step which is unphysical.

Would it be preferable to calculate snGrad instead and interpolate the conductivities to the cell faces? Or is it an error I can't avoid? Is it actually producing significant problems when the grid is small enough?

b) Should already be answered in a), everything is calculated at the cell centers right now.

c) The ZeroGradient BC was applied to the conductivity. I believed this would yield better results at the boundaries because it should avoid this step problem I explained in a), but it didn't converge to a reasonable result. The method I linked to in the first post works though.

[RodriguezFatz]

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.

[chriss85]

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.

[RodriguezFatz]

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?

[chriss85]

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

[RodriguezFatz]

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?

[chriss85]

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.

[chriss85]

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?

[RodriguezFatz]

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?

[chriss85]

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

[RodriguezFatz]

to 2) Can you try "Gauss linear" and "linear"?

[RodriguezFatz]

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?

[chriss85]

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...