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

Question

---

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
1. Some pre-loading of variables

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;

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

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

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

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

Answer 1

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.