# For a determined (known) Space charge density, what are the conditions to obtain the Electric potential/field distribution? (COMSOL, MATLAB)

Theoretic part

From the theory, in Electrostatics inside a real dielectric material between real conductors, in a simple 1D plane geometry between points $$P1$$ and $$P2$$, according to the current density continuity equation, when an inequality occurs between the flow of charges into a region and the flow of charges outside then charge accumulates in that region.
According to Gauss's law, a space charge field $$E\rho$$ can be associated to that charge distribution,

$$\nabla(\epsilon \cdot \boldsymbol{E}\rho)=\rho$$ ..... (Eq. 1)

Therefore, the electric field $$Etot$$ within the material, in the presence of space charge $$\rho$$ is given by the sum of two contributions: the space charge field $$E\rho$$ (generated by the existence of $$\rho$$) and the external field $$E0$$ (also called the Laplacian field, when $$\rho = 0$$), which is induced by the applied voltage (or electric potential difference at the boundaries):

$$Etot = E\rho + E0$$ ..... (Eq. 2)

* Formulas were adapted from the theory described in paragraph 5.1 (page 4/73) and paragraph 5.6 (page 22/73) of this link: HERE .

Conclusion: If one has a value for the $$\rho(x)$$ distribution and knows the $$\epsilon$$ (assume its constant) he can determine the electric field $$E\rho$$ from eq: $$\nabla(\epsilon \cdot \boldsymbol{E}\rho)=\rho$$ , as

$$E\rho(p)=(\int_{P1}^{p}{\rho(x) dx})/\epsilon$$ ..... (Eq. 3)

where $$p$$ is a point between $$P1$$ and $$P2$$.

If he finds out $$E\rho$$ and adds $$E0$$ to it, he can determine the total field $$Etot$$ (lets call it simply $$E$$), and then proceed to find out the electric potential distribution as

$$V(p)=-(\int_{P1}^{p}{E(x) dx})$$ (from $$E = -\nabla V$$ ) ..... (Eq. 4)

What I don't understand is: If you know the value for $$\rho$$ and $$\epsilon$$ then you automatically know the distribution of the electric field $$E\rho$$, or of the total field $$Etot$$?
-- think of it like this: if you know the value of $$\rho$$, then the contribution of the external field (the surface charge distributions at the boundaries) is/are already quantified inside of $$\rho$$? and the $$E\rho$$ term in (Eq. 1) is actually $$Etot$$ directly, so by solving (Eq. 3) I find directly the electric field of interest?.

Practical example: I have two different $$\rho(x)$$ distribution (Pastebin Data 1 and Pastebin Data 2) obtained for a material of 0.5 mm in length (from $$P1$$ at 0mm to $$P2$$ at 0.5mm), with $$\epsilon$$ = 2.85*$$\epsilon0$$, for an applied voltage of 5000 V
(at $$P1$$ the potential is $$Vapp$$ = 5000 and at $$P2$$ the potential is 0).

I note x0 as the first column from Pastebin Data 1, y0 as second column from Pastebin Data 1. Same for x60 and y60, from Pastebin Data 2.

HOWEVER, things get a bit more complicated when I test this out, because I get different results when solving it numerically and via simulation. WHY? Attached results for $$E$$ can be found HERE.

A) Numerical testing in MATLAB

data1 = importdata('data1.txt');
data2 = importdata('data2.txt');
x0 = data1(:,1);
y0 = data1(:,2);
x60 = data2(:,1);
y60 = data2(:,2);
Vapp = 5000;
P1 = 0;
P2 = 0.5e-3;

1. I can determine the Laplacian contribution to the field ($$E0$$) (the expected 10 MV/m) as:

E0 = Vapp/(P2-P1);

and the Laplacian potential as:

for k = 1:length(x0)
V0(k) = Vapp.*(x0(k)-P2)./(P1-P2);
end

1. For the 1st $$\rho(x)$$ distribution, I determine $$E\rho$$ with (Eq. 3) and then the potential $$V$$ numerically with (Eq. 4) where $$E$$ is just $$E\rho$$ (so no contribution of $$E0$$), between points $$P1$$ and $$P2$$. I check the potential on $$P1$$ and $$P2$$ and its 0 V and approx. -1900 V.

2. For the 2nd $$\rho(x)$$ distribution, I determine $$E\rho$$ and then the potential $$V$$ numerically in same conditions. I check the potential on $$P1$$ and $$P2$$ and its 0 V and approx. +1000 V.

3. Figures with the distribution of $$\rho(x)$$, $$E\rho$$ and $$V$$ are uploaded here. The code I used in Matlab is:

% to find out E
for k = 2:length(x0)
field0(k) = trapz(x0(1:k),y0(1:k))./(2.85.*8.854187e-12);
end
for k = 2:length(x60)
field60(k) = trapz(x60(1:k),y60(1:k))./(2.85.*8.854187e-12);
end

% to find out V
for k = 2:length(x0)
voltage0(k) = -trapz(x0(1:k),field0(1:k));
end
for k = 2:length(x60)
voltage60(k) = -trapz(x60(1:k),field60(1:k));
end

To these two, I also have from above the $$E0$$ and $$V0$$.

B) Simulation testing in COMSOL
1. I assume 1D geometry and a domain from P1 (0) to P2 (0.5 mm) on $$x$$, representing material M that has, in Electrostatics, the constant relative permittivity of 2.85.
2. I then add in the physics a Space charge density node as I have a distribution of $$\rho(x)$$ that I determined experimentally for that applied voltage of 5000 V.
(I can know $$\rho$$ by doing some experimental PEA measurements and I can import the $$\rho$$ data as an Interpolation function that varies with $$x$$.)
3. In COMSOL I now need to consider the boundaries. If I add a Floating potential boundary at point $$P1$$ then the results I get in COMSOL are those for $$E\rho$$. If I add the Electric potential boundary with the potential equal to 5000 V, I do get a solution for $$E$$ that is very similar to what I should get for $$Etot$$, with a mean electric field of 10 MV/m. But it always has a mean of 10 MV/m, so does this mean that the influence of $$E\rho$$ on $$Etot$$ is always 0, since the mean of $$E\rho$$ is always 0?

*1) If I add a Floating potential boundary at point P1

If I add a Floating potential boundary at point P1 then the field distribution is the one obtained only from the term $$E\rho$$ from (Eq. 2), and is identical to the one computed numerically from $$E\rho(x)=(\int_{P1}^{p}{\rho(x) dx})/\epsilon$$ . However, its mean is not 10 MV/m and the potential is not 5000 V at point P1 (it can be whatever else). But the potential at the boundaries should be 5000 and 0 V for my case in order to have an exact/correct solution, no?

*2) If I add an Electric potential boundary of V = 5000 V

If I add an Electric potential boundary at point P1 (of V = 5000 V) then the field distribution is probably given by the term of $$Etot$$ from (Eq. 2).
So it is NO LONGER just $$E\rho$$ computed numerically above, BUT IT IS ALSO NOT $$Etot$$ from (Eq. 2), as $$E0$$ + $$E\rho$$ distributions (that were both computed numerically above).
The electric potential is fixed at 5000 and 0 V at the boundaries BUT IT IS NOT the sum $$V0$$ + $$V$$ (both were computed numerically above).

*3) The resulting figures for $$E$$ and $$V$$ from the simulations of *1) and *2) are uploaded here.

What I don't understand What I get in COMSOL appears to be correct, but it always have the mean of $$E$$ as the value of $$E0$$. Does this mean that the influence of $$E\rho$$ is zero? Also, why am I not able to get to similar results via the numerical approach? Not even if I add $$E\rho$$ and $$E0$$ to get to the $$Etot$$, still not same result as the one I get from *2).

TLDR:

1. Is there such a thing as what's described in (Eq. 2)? $$Etot$$ = $$E\rho$$ + $$E0$$ if I have a known $$\rho$$ distribution? Or by simply knowing $$\rho$$ I directly know $$Etot$$, since the contribution of $$E0$$ (the Laplacian contribution of the surface charges at the boundaries) is already included in $$\rho$$.

2. Why do I not get the same result when I solve the equation numerically and by simulation, if I know $$\rho(x)$$ and $$\epsilon$$ ? What more should I do to my numerical approach to yield the same results from *2) (the ones that I expect to have) and not those from *1) (the ones I actually get)? Even if I add the Laplacian contribution to my space charge contribution I still do not get the same. (SEE HERE.)

3. Are the results I got from *2) actually correct, considering that the mean of $$E$$ is always $$E0$$ no matter what $$\rho$$ distribution I add (so a different $$E\rho$$ but zero influence on $$Etot$$)?

4. Should I solve this numerically via another way (not via simple numerical integration with the trapezoidal rule)? I assumed that if I solve it via trapz then I'm not missing any extra constant of integration, but I am unsure now about the whole thing.

• It's really difficult to understand what you really want to accomplish here. Is it possible for you to edit your question and make it more clear? – Alone Programmer Sep 29 '19 at 0:26
• Thank you for the feedback. I edited the question, please check it again. – Lucian Viorel Sep 29 '19 at 9:56
• It's not even more difficult to understand what you are asking... Please make it really brief and objective and before that please do some research about basic of Maxwell equations and numerical methods to solve them. – Alone Programmer Sep 30 '19 at 12:55

I think you have misunderstanding about what the Gauss's law and subsequently Poisson equation means. You have these relations (Gauss's units):

$$\nabla \cdot \mathbf{E} = 4 \pi \rho$$

$$\mathbf{E} = -\nabla \phi$$

$$\nabla ^{2} \phi = -4 \pi \rho$$

$$\phi(\mathbf{r}) = \int_{\Omega} \frac{\rho(\mathbf{r^{'}})}{|\mathbf{r}-\mathbf{r^{'}}|} d^{3} \mathbf{r^{'}}$$

Where $$\mathbf{E}$$ is the electric field, $$\rho$$ is electric charge density, and finally $$\phi$$ is electric potential and $$\Omega$$ is your domain (in your case 1D line).

There is no such things as electric potential caused by Laplacian term and electric potential caused by charge density! There are some possibilities for the fact that when $$\mathbf{E}$$ is not zero and you have electric field such as but not limited to:

1. Volumetric charge density ($$\rho$$) inside a domain ($$\Omega$$) and zero electric potential at the boundaries ($$\partial \Omega$$) as $$\phi = 0$$.

2. Zero charge density ($$\rho = 0$$) inside a domain ($$\Omega$$) and non-zero electric potential at the boundaries ($$\partial \Omega$$) as $$\phi = \phi_{0}$$.

Both of these two situations still could be solved easily by Poisson equation and you will get one $$\mathbf{E}$$ as electric field and one $$\phi$$ as electric potential. Also, you should understand that you should choose the boundary conditions based on your physical situation. I hope it helps you to understand what Gauss's law means and there is no such conception as $$\mathbf{E}_{\rho}$$ caused by electric charge density and $$\mathbf{E}_{\nabla^{2}}$$ caused by Laplacian term.

• Could you please explain your rationale in contrast to what is described in paragraph 5.1 (page 4/73) and paragraph 5.6 (page 22/73) of this link here: ocw.mit.edu/resources/… ? This is where I got the notion of electric potential caused by Laplacian term and electric potential caused by charge density - which you noted that it doesn't exist as described. Am I missing something? Please help. – Lucian Viorel Oct 2 '19 at 11:26
• That paragraph 5.1 in page 4 is talking about the fact that you can add any function that satisfies the Laplace equation and boundary conditions to the solution of Poisson equation cause: Let's say you have this Poisson equation $\nabla^{2} \phi = -4\pi \rho$ and a function that satisfies the Laplace equation: $\nabla ^{2} \psi = 0$, if I add these two equations: $\nabla ^{2} (\phi + \psi) = -4\pi \rho$ and I can call: $\phi^{'} = \phi + \psi$ and finally again I will have Poisson equation: $\nabla^{2} \phi^{'} = -4\pi \rho$. But, I'm not sure what do you want to extract from this fact... – Alone Programmer Oct 2 '19 at 13:11
• That makes more sense now. As to your last paragraph, I am trying to backtrack my experimental results and see if they are actually correct. For a material and an applied electric field I determined the $\rho$ distribution for it at that electric field. So I have a $\rho$ distribution that I know was determined under particular conditions (known voltage at boundaries, known applied field etc), and now I want to do the backward and answer the question ... for this $\rho$ distribution that I have, if I solve the equation, do I get this field back? Can't seem to, and I do not understand why. – Lucian Viorel Oct 2 '19 at 13:37
• "For a material and an applied electric field I determined the $\rho$ distribution for it at that electric field." how did you do that? did you extract the $\rho$ from the Gauss's law? If yes and if you did the math correct, now if you put the $\rho$ back into your Gauss's law you MUST get the $\mathbf{E}$ back. – Alone Programmer Oct 2 '19 at 13:42
• I extract $\rho$ by actual measurements called PEA (Pulsed electroacoustic), used non-destructively to study space charge ($\rho$). It is not exact and it adds noise, but the results are accepted if obtained correctly. So I must get E- and I actually do, if I solve via Simulation: I import the $\rho$ distribution in the software and add potentials at boundaries. However, if I try to solve numerically via trapezoidal integration of $\rho(x)$ (Eq. 3), I get the same shape but with an offset. I assume is because of that $ψ$ that I don't account for. See imgur.com/a/JA8gWPw – Lucian Viorel Oct 2 '19 at 14:20