An Efficient Numerical Technique for the Simulation of Charge Transport in Polymeric Dielectrics

In recent years, the use of HVDC cables has grown exponentially. One of the main challenges that remains concerns the space charge accumulation inside the insulating materials. A better understanding of the mechanisms governing this phenomenon is essential to improve the performance of HVDC systems. Numerical simulations are often employed to achieve this goal. For this reason, it is important to perform them in an efficient way. In this work, we test several numerical techniques, aiming to assess which one is the best to use for fast and reliable simulations. We consider a well-known bipolar dynamic model from the literature for our simulations. The model considers a single level of deep traps and is implemented in a one-dimensional Cartesian coordinate system, considering a thin specimen of polymeric material. We compare three different time discretization methods: a fully explicit, a semi-implicit, and a fully implicit approaches. For the advective flux discretization, we compare the first-order upwind scheme (FOU) with a second-order upwind scheme coupled with the Koren flux limiter (SOU/KL). Regarding the computation of the polarization current, we introduce a simple approach using Sato’s equation and compare it with the well-established approach based on the total current density.


I. INTRODUCTION
In the last decades, high voltage direct current (HVDC) cables have become increasingly popular.There are many factors that contribute to the exponential growth in the use of this kind of technology.In general, these systems allow for efficient and reliable long-distance energy transmission.For this reason, these are often used as cross-country interconnections, among many other applications.However, a number of challenges related to HVDC systems still exist and need to be addressed [1].One of these critical issues is the space-charge accumulation within the extruded cable insulation [2].The dielectric materials used to insulate HVDC systems are usually solid polymers, such as low-density polyethylene (LDPE) or cross-linked polyethylene (XLPE).The charge trapped inside these materials can significantly The associate editor coordinating the review of this manuscript and approving it for publication was Roberto C. Ambrosio .alter the electric field distribution within the cable insulation.This may result, in turn, in increased electrical stresses on the dielectric, leading to accelerated aging and premature system failures [3].For this reason, a deeper understanding of charge transport and accumulation processes inside solid polymeric dielectric materials is needed to improve the performance of HVDC cables [4], and allow higher design temperatures and electric fields, while maintaining the appropriate reliability levels.Experimental measurements are a fundamental tool to achieve this goal.The pulsed electro-acoustic (PEA) method is one of the most widely used techniques to measure the space-charge distribution inside a dielectric sample [5].Such measurements are performed by applying voltage pulses to the material.This creates detectable acoustic signals due to the internal space-charge motion which, when measured, provide information on the space-charge distribution.This technique has been applied for decades on flat specimens of dielectric materials [6], [7], [8].Only in recent years PEA measurements started to be performed on mini and full-size cables [9], [10].Another well-known diagnostic technique is the polarization/depolarization current measurement [11].One drawback of measurements on dielectric materials is that these can be expensive and time-consuming.For this reason, the development of physical models is a fundamental tool to complement experiments.These models can be grouped into two categories, i.e., analytical and numerical models.Examples of the first category can be found in the works of Chen and Xu [12] on trapping and detrapping mechanisms, and in the work of Mazzanti et al. [13] on assessing apparent trap-controlled mobility and trap distribution.However, these models are based on several simplifying assumptions; for example, [12] considers volume-averaged quantities, whereas the model in [13] is directly based on measurements.Models aiming to provide a space-and time-resolved description of the charged species involve the solution of systems of partial differential equations (PDEs), which in general do not admit analytical solutions in closed form for the quantities of interest.In this case, it is necessary to set up numerical simulations to compute the solution of the discretized charge generation and transport equations.Regarding this case, one of the first bipolar charge transport (i.e., considering both electrons and holes) models for dielectric materials was proposed by Alison and Hill [14].They considered a constant charge injection from both electrodes and, regarding the extraction, charge was free to leave the material unimpeded.Constant trapping and recombination coefficients were used and the detrapping process was neglected.Fukuma et al. added a Schottky emission law for the charge injection and considered an energetic barrier for charge extraction.They also considered the effect of temperature changes over time [15].Later on, Kaneko et al. proposed a model to reproduce the packet-like space charge in XLPE [16].They considered an electric field dependent mobility.The governing equations were solved by means of a finitedifference scheme, whereas the discretization technique is not explicitly discussed in other referenced works.Le Roy et al. presented a model to describe bipolar transport in LDPE, using an effective mobility and a single trap level [17], [18].Again, the detrapping process was not considered.In these works, they focused also on the polarization current that was computed using the total current density.The physical contributions from transport and species generation/depletion in the continuity equation are accounted for with an operatorsplitting technique.The detrapping process was added in a later work by Le Roy et al., where the developed model was used to study electroluminescence in polyethylene-based materials [19].In this work, an initial number density of carriers, uniformly distributed throughout the domain, was considered to obtain the best fit between experiments and simulations.Subsequently, Baudoin et al. used a similar model to study the current versus applied electric field characteristic at steady state, using an under-relaxation technique to avoid divergence in the iterative solution of strongly nonlinear equations.[20].Later on, Le Roy et al. presented a model to study charge accumulation inside electron-beam irradiated low density polyethylene [21].They considered the extraction of charges possible at the grounded electrode through an ohmic law.In all the above works by Le Roy and colleagues (except [20], where a steady-state analysis is performed), the continuity equation is integrated over time by means of a fully explicit method (Euler's method).The Poisson equation, used to compute the electric potential, is solved with the boundary element method (BEM).The transport equation is discretized using the finite-volume (FV) approach.In particular, the drift term is approximated by a third-order upwind scheme (QUICKEST), coupled to a flux limiter (ULTIMATE) to avoid negative number densities values caused by numerical oscillations.Other discretization techniques have also been successfully applied to the problem of modeling bipolar charge transport.For example, in [22], [23], and [24] the finite element method (FEM) was used for spatial discretization, coupled with a Runge-Kutta time integration scheme.Meftali et al. focused on the effect of the injection barrier height on the results of the bipolar transport model using finite differences and an explicit time integrator [25].Recently, Doedens et al. proposed a comprehensive model to describe the effect of electrodes surface roughness on the space charge accumulation phenomenon inside XLPE [26], [27].They considered also Fowler-Nordheim injection processes in addition to the well-established Schottky law.The diffusion of charge carriers was also added to the model.The governing equations are solved using COMSOL Multiphysics.A similar analysis was conducted in [28], where the effect of nanometric scale processes at play at the electrode interfaces was studied.In this work, 2D simulations were carried out using COMSOL Multiphysics.Regarding the charge extraction process, all the listed models assume non-blocking electrodes, meaning that the charge is free to leave the dielectric under the electric field effect.Charge accumulation is often modeled using a single level of traps although there are works where multiple levels of traps were considered, and also continuous trap distributions [29], [30].These models were often employed for simulations in 1D or 2D Cartesian coordinates, but also for cable geometries [31], [32] in cylindrical coordinates.The models described so far are fluid, but other approaches are possible.For example, Cambareri recently developed a circuital model to describe polarization and depolarization currents [33].His model employs the least possible number of parameters to fit experimental measurements and considers a single carrier type with negative charge.
In this work, we compare the results yielded by a Cartesian one-dimensional fluid model when using different numerical techniques to solve the governing equations.The analysis is focused (but not limited) on polarization current results.The goal of this work is to determine which numerical techniques should be employed for fast simulations of charge transport and accumulation inside dielectric materials, since the charge transport problem is not always studied by numerical specialists.Similarly to what was done in [34], we propose an analysis of two different discretization schemes for the advective fluxes.In our case, the schemes are the first-order upwind scheme (FOU) and a second-order upwind scheme with the Koren flux limiter (SOU/KL) [35].When comparing the accuracy and computational time of the schemes in the context of a specific physical case study, our findings diverge from those reported by Le Roy et al. in [34].We also add an analysis of three different time-discretization schemes: a full explicit method, a semi-implicit method, and a full implicit method.In addition, we show an alternative expression to evaluate the polarization current, which allows to avoid the displacement current calculation (required if the total current density approach is used) and that is based on Sato's equation [36].The code used to produce the results in this work is open-source and available online at the link https://github.com/PTL-Unibo /CALLIOPE.The developed model is presented in Section II.The numerical techniques that will be compared are discussed in Section III.In Section IV, one of the methods for the post-processing computation of the polarization current, Sato's formula, is described.The comparison between the results obtained from the simulations is presented in Section V.

II. PHYSICAL MODEL
In this work we developed and implemented a fluid model for the bipolar charge transport in solid dielectric materials under the effect of electric fields.This model shares many features with those described in Section I and, in particular, with those developed by the group of Le Roy and colleagues and discussed in several works, including [17], [19], [21].The approach is based on the band theory model for dielectric materials, which describes the conduction processes by means of a valence band and a conduction band [37]; for a dielectric material, these bands are separated by a forbidden region (band gap), characterized by a wide energy gap.Due to physical and chemical imperfections, energy states are introduced in the band gap [38], [39].An electron takes part in the conduction process when located in the conduction band, whereas a hole is available for the conduction when located in the valence band.The energy states in the band gap are usually referred to as traps, since electrons and holes that occupy them are excluded from the conduction mechanism.We made the hypothesis that a single level of deep traps is present, for both electrons and holes, and therefore the model described in this work considers four types of charge carriers: free holes h µ ; free electrons e µ ; trapped holes (h t ); trapped electrons (e t ).The fast dynamics related to shallow traps are modeled with a macroscopic transport parameter, the mobility [27] .The species' number density time evolution is described by means of a continuity equation: In (1), n s , s and s are the number density, the flux density and the source term of the given species s; the flux density is given by the sum of two contributions: the first one is due to spatial concentration gradients, whereas the second is due to the drift induced by the electric field.The expression used to evaluate the flux density is, therefore: In (2), K s indicates the diffusion coefficient, q s /|q s | is the sign of the species' charge, µ s is the electrical mobility and E is the electric field.The diffusion coefficient is evaluated using Einstein's relation: In (3), T s is the species' temperature, k B is the Boltzmann constant and e is the elementary charge.The source term in (1) represents the net number of charge carriers that are generated or lost per unit-volume and per unit-time due to elementary processes.Indeed, a free carrier can be trapped and removed from the conduction process; a previously trapped charge can become free due to thermal activation, and recombination can occur between charges with opposite sign.
For each of the four carriers, the source term s is computed using the following equations:  [27].In this way it is possible to better represent the non-linear behavior of the material at higher electric fields.Still, it should be noted that even with constant rate coefficients the model is already non-linear due to the presence of a product between the number density and the electric field in (2) and the presence of products between number density in the source terms (4), ( 5), (6), and (7).A modified Schottky emission law is employed to compute the fluxes of electrons and holes injected from the electrodes into the dielectric: In (8), a is the Richardson constant, divided by the elementary charge, and φ is the magnitude of the energetic barrier located at the interface between the electrode and the dielectric; β E is a constant that can be expressed as: where ε 0 is the vacuum permittivity and ε r is the relative permittivity.The electric potential is determined from the Poisson equation where ϕ and ρ are the electric potential and the charge density, respectively.The model considers four species, hence the net charge density is given by: The model couples the transport equations with the Poisson equation.Therefore, it is capable of reproducing local enhancements of the electric field due to space charge accumulation.This can be used to assess physical conditions that may lead to local breakdown.

III. NUMERICAL MODEL
The physical model described in the previous section has been applied to a one-dimensional geometry, using the finitevolume method.The domain is subdivided into a finite number of cells; inside each cell there is a node (domain point), where the physical quantities are evaluated under the assumption that nodes store cell-averaged values.When a non-uniform cell spacing is considered, interfaces between adjacent control volumes are placed at the midpoint between nodes.This choice is not the only possible, as one could chose instead the position of interfaces such that each node is exactly at the midpoint of its control volume.We chose the first approach over the second to maintain second-order accuracy when centered finite difference formulas are used to estimate physical quantities [40].For this work, a constant and spatially uniform mobility is considered.The diffusion coefficient, calculated with (3), has the same properties; the trapping, detrapping, and recombination parameters and the deep trap density are also considered constant and uniform.We implemented three different versions of the model, corresponding to three different time-discretization schemes: a full explicit method, a semi-implicit method, and a full implicit method.In addition, two different discretization schemes for the fluxes at the interfaces between adjacent cells were tested: the FOU and the SOU/KL [35].
Equation ( 1) is integrated over the i-th cell of the domain, giving the following: In ( 12), V i is the volume of the i-th cell of the domain and the symbol S V i is used to denote the boundary surface of that volume.We are considering a one-dimensional geometry, thus the surface integral appearing in ( 12) can be evaluated considering only the left and right boundary surfaces, respectively , obtaining: The subscript i + 1 2 in the expression is referring to the interface between the adjacent cells i and i + 1.In the case of a uniformly spaced domain, ( 13) may be rewritten as where is the uniform spacing between domain points.The first and last domain interfaces represent the electrodes.Holes are injected into the dielectric material from the positive electrode, whereas electrons are emitted from the negative electrode.The fluxes of injected carriers are computed using (8).Regarding the extraction of carriers, blocking electrodes are considered; this implies that charge carriers cannot escape from within the dielectric but can only accumulate at the electrode interface.This is different from what is usually done in the literature, but in the simulations we conducted, considering blocking electrodes produced results very similar to the ones obtained when non-blocking electrodes were considered.The electric potential is computed solving the discretized version of Poisson's equation: A. FLUX DISCRETIZATION In this section we introduce two different strategies for the numerical discretization of advective fluxes, that will be compared in Section V.

1) FOU -FIRST-ORDER UPWIND SCHEME
Using the FOU, the discretized flux density takes the form: ) In ( 16), we introduced the notation: 12548 VOLUME 12, 2024 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The symbol u i+ 1 2 ,s indicates the velocity of a carrier for the s species at the interface i + 1 2 due to the electric field drift: As is well known, the FOU is affected by considerable numerical diffusivity.To overcome this issue, one may employ higher-order schemes.This strategy, helpful in lowering numerical diffusivity, may result in oscillations in domain regions with strong gradients, possibly leading to negative number densities [41].A well-known technique to improve on this behavior is the use of a flux limiter.There is a vast literature on flux limiter types and their properties, see, e.g., [42].For this work, we tested the SOU/KL since it has been successfully used in state-of-the-art codes such as the Afivo software for streamer propagation [43].

2) SOU/KL -SECOND-ORDER UPWIND SCHEME WITH THE KOREN FLUX LIMITER
Convective flux densities computed with a flux limiter at the interface between adjacent cells are usually expressed as: is the flux density evaluated using a low-order scheme and γ (high) i+ 1   2   represents the flux density evaluated using a high-order scheme.The function and the expressions for γ (low) i+ 1   2   and γ (high) i+ 1   2   depend on the specific flux limiter considered.The Koren flux limiter function is defined as [35]: In (21), r is the ratio between adjacent gradients of the advected physical property (in our case the number density) and depends on the velocity direction.The above expression is also equivalent to the piecewise-defined function: The flux limiter function KN (r) is shown in Fig. 1.
Piecewise definition was used in the numerical implementation in order to avoid the costly evaluation of the max and min functions and improve code performance.20) is evaluated with a low-order scheme, i.e., a firstorder upwind.γ (high) is instead evaluated with a second-order  upwind scheme: Note that ( 24) is only valid in the case of uniformly spaced domain points.A schematic representation of the considered domain is shown in Fig. 2. Combining ( 20), ( 23), (24) and considering the diffusion, the discretized flux density is: From ( 26) and ( 27) it is possible to see that the ratio r i+ 1 2 has different expressions depending on the sign of the advection velocity.In general, when evaluating the flux for the interface i + 1  2 , the denominator is given by the difference between the two nodal number densities adjacent to the interface upstream of i + 1 2 .The numerator is, instead, always equal to the difference between the two nodal number densities adjacent to i + 1 2 .Clearly, the evaluation of (25) requires more computational resources than (16), since the Koren flux limiter function, , needs to be computed.In problems dominated by strong advection, this additional cost is often justified by the improved accuracy with respect to the results obtained when using the FOU [41].The code implementation of ( 25) is achieved efficiently using the following vectorized Matlab function: This function is adapted from the Fortran implementation in [44], where loops are employed.Vectorization of the Matlab code allows the computational time required to be significantly reduced.However, the instructions to obtain ii and jj are computationally burdensome.Such instructions are not required in the code implementation of ( 16), making it faster.

B. TIME DISCRETIZATION
In this section we discuss and compare the performance of three different time-discretization schemes for the advectiondiffusion-reaction equations describing the evolution of charged species.

1) FULLY EXPLICIT METHOD
For the fully explicit integration method, the time derivative in the first member of ( 14) is approximated with a finite forward difference using the forward Euler method.The discretized expression reads as and thus it is possible to compute the number density at the next time instant explicitly as: The time step duration t can't be arbitrary, as the Euler method is conditionally stable [40].t should indeed be evaluated at each time step to avoid numerical instabilities.A Von Neumann stability analysis is required to rigorously assess the right value of the time step needed to avoid instabilities.For a convective-diffusive problem, with constant convective velocity u, the CFL stability condition is with C m = 1.This condition can also be used to estimate the time step for the problem described by (29), using C m < 1.The appropriate value for C m is often determined by trial and error.This approach is not rigorous but allows to simplify the computation and to avoid repeating the Von Neumann stability analysis, for additional terms or different discretization schemes.For this work, we propose an alternative method to evaluate the time step, which takes into account source terms.If the flux density is discretized using the FOU, we can write (29) as: The coefficients α, β, γ and δ depend on the particular kind of species that is being considered.In the case of free holes, we have: Since α, γ , δ are always greater than zero, the stability condition that we propose is determined by imposing β > 0, resulting in Note that (36) provides a time step for each cell of the domain.Similar expressions can be derived for the other model species, also yielding a time step for each point of the domain.For free electrons: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
For trapped holes: For trapped electrons: In total, we obtain a time step for each domain cell and for each species.The minimum of these values is the one actually used in the algorithm.This approach is more accurate than the one proposed in (30).The two methods were implemented and tested using the same initial conditions.We used a value of C m = 0.5 for the simpler stability condition.The computational time required by the simulation using our stability condition is approximately 33% lower than the time required by the simulation using the simpler stability condition.

2) SEMI-IMPLICIT METHOD
In order to develop a semi-implicit approach, the contribution due to the source term is split from the one due to the drift and diffusion mechanisms (operator splitting technique).
We used the superscripts and to indicate the number density values obtained considering only the drift and diffusion contribution, and the source term contribution, respectively.
a: DRIFT-DIFFUSION CONTRIBUTION Equation ( 40) may be rewritten in matrix notation as We will use square brackets and curly brackets to indicate respectively matrices and arrays from here on in this section.
The fluxes have been rewritten isolating the number densities for the diffusive flux contribution and the electric potential for the advective flux contribution, respectively.Note that the cell volumes are incorporated into the matrices K D s and K µ s .We apply the Crank-Nicholson integration scheme, which is second-order accurate in time and is unconditionally stable [41], for the diffusive operator.Similarly to what has been done in [45] and [46] we evaluate the drift term using a semi-implicit scheme.This has the advantage of making the method robust with respect to rapid electric field variations.In contrast to [45] and [46], where the electric potential was evaluated at the time-instant k + 1, we evaluate it at k + 1 2 ({ϕ} 2 ) for coherence with the Crank-Nicholson scheme.
The resulting numerical scheme is then: The above matrix equation represents a linear system that can be solved at each time-instant to compute the number density of the species at the next one.
To estimate {ϕ} k+ 1   2 , we combine the matrix form of ( 15) and ( 11) and we apply a forward finite difference to ( 42) Note that the approach used in ( 45) is not fully rigorous since {n s } and {ϕ} in the right-hand-side are not evaluated at the same time-instant.Isolating the term {n s } k+ 1   2 from ( 45), and substituting it in (44), yields: The above system yields the electric potential at the half step required to solve (43).Note that the drift and diffusion contributions regard only mobile species.The described procedure, therefore, is carried out only for mobile holes and mobile electrons.

b: SOURCE TERM CONTRIBUTION
The second member of ( 41) is expanded using Taylor series, neglecting higher-order terms: The chain rule is applied to obtain Equation ( 48) can be rewritten as where we denoted Rewriting the first member of ( 49) with a centered finite difference we get Rearranging (51) yields: Using the above expression, also employed in [47], one may compute the number density variation of each species due to the effect of the source term alone.This process needs to be repeated for each domain point.The time step t for the source term integration is determined by imposing the condition and therefore: A schematic representation of the calculation algorithm is shown in Fig. 3. Starting with known number densities throughout the domain at the generic time instant k, the electric potential is computed with (46).The number densities at the next time instant k + 1 due to drift and diffusion fluxes are evaluated with (43).Then, the number density variation of each species for each cell of the domain due to the contribution of the source term alone is computed with (52).n and n are added together to yield the updated total number densities at the next time instant.The process is then iterated to advance in time.

3) IMPLICIT METHOD
We based our implementation of an implicit solver on the implicit adaptive scheme TR-BDF2 (trapezoidal rule / backward differentiation formula) [48], [49] provided by the ode23tb function from Matlab.This function is capable of solving stiff systems of differential equations that can be written in the form where {y} represents the state vector.This vector uniquely determines all the relevant quantities at a given time instant.
The user needs to provide f , a function that takes in input the time instant and the state vector: the number densities of all species over the entire domain.All the relevant physical quantities, e.g.charge density, electric field, flux density at domain interfaces, can be computed from the state vector.The output of this function is the time derivative of the state vector at the considered time instant, i.e., ∂n/∂t in (1).Rearranging (14) to isolate the time derivative of the number density, one gets the following: The right-hand-side of ( 56) represents the computation that the function f needs to perform.Since the JIT Matlab compiler is not always able to handle loops with complex expressions efficiently [50], [51], the implementation of the function f was extensively vectorized, aiming at reducing computational time.A drawback of using the implicit approach is that the output of the ode23tb function is the state vector at every required time instant.Therefore, the computation of non-state variables, such as electric field or flux density at domain interfaces, needs to be performed again after the simulation ends if physical quantities other than the number densities are needed.For this reason, an additional post-processing phase is required in the implicit implementation.This, however, does not affect significantly the computational efficiency of our code thanks to the aforementioned vectorization process.The ode23tb function requires the estimation of a Jacobian matrix and therefore increasing the number of domain points may result in a significant increment in the computational time required by the simulations.This aspect represents a potential limitation of the proposed approach, but for domains with approximately one-hundred points the computational time was deemed acceptable.

C. SOFTWARE STRUCTURE
The basic steps of a simulation are illustrated in Fig. 4. The first thing to do is to create a parameter structure P.This structure contains the details regarding the geometry, the material and the value of electron and hole mobility, as well as the parameters appearing in ( 4), ( 5), ( 6), (7).The time instants at which the solution is sought can be specified using the vector time_instants.At this point the user can provide the specific options for the simulation, for example it is possible to set the flux discretization scheme or the time integration scheme.At this point the simulation can be launched using the Run function.When using the implicit time-integration method an additional post processing phase is required to compute the evolution in time of the charge density, the polarization current and the electric field.The final output is a structure out containing all the relevant quantities ready to be plotted.

IV. SATO'S FORMULA
The performance of the three time-discretization schemes described above will be benchmarked based on the results of the simulated polarization current in Section V. To compute the polarization current, one must first obtain the total current density J tot at each control volume interface at each time instant.We start by recalling the common technique to compute J tot and then introduce a simplified approach based on Sato's formula, originally developed for gas discharge simulations [36].The standard approach to evaluate J tot is to compute the total current density at each domain interface and at a given time instant.The polarization current value at the considered time instant is then obtained by averaging J tot over the entire domain.In our simulations we also computed the polarization current using Sato's formula [36]: In (58), I p indicates the polarization current, V a is the applied voltage and E s is the static electric field.The static electric field, also known as Laplacian, is obtained considering only the applied voltage and neglecting the charge density inside the volume.In a one-dimensional domain, the equation becomes where J p is the polarization current density and L is the length of the specimen.It should be noted that the above equations are valid only in the case of a constant applied voltage between the electrodes.Equation (59) states that the polarization current can be computed as the mean over the specimen of the conduction current, J cond .The advantage of Sato's formula is that the evaluation of the displacement current is not required.This results in a shorter computational time and a simpler implementation.More details on Sato's equation can be found in Appendix.

V. RESULT COMPARISON
In this section we compare the performances of the three time-discretization schemes (Sections V-A,V-B), two different flux discretization schemes (Section V-C) and of the two above methodologies for the current density computation (Section V-D).The physical parameters for all the simulations conducted are taken from the work in [18], where a similar model was used for a charge transport study on polyethylene.We refer the reader to [18] for a tabular description of the parameters.A sample with thickness of 0.40 mm and subjected to a voltage difference of 4.00 kV was considered.We simulated a time interval of 1.00 • 10 5 s with the three solvers, focusing on the polarization current results.

A. IMPLICIT VS. SEMI-IMPLICIT TIME DISCRETIZATION SCHEME
Here we compare the accuracy obtained when simulating the same problem with the same flux discretization scheme (FOU, ( 16)), using the implicit and the semi-implicit timediscretization schemes.Accuracy-wise, throughout this study the results yielded by the fully implicit solver were taken as a reference for the other two schemes for two main reasons: first, using an implicit approach, there is no loss of accuracy since no operator splitting is performed.In other words, the flux and source terms are integrated simultaneously without the loss of accuracy inherently bound with fractional time-stepping techniques where different numerical schemes are employed to discretize the two mentioned contributions.Second, the ode23tb uses an adaptive time-step length, which allows to control the absolute and relative accuracy throughout the simulation.The results are shown in Fig. 5.It is possible to observe that there is a good agreement between the results obtained using the implicit integration scheme and the semi-implicit method.In particular, the error is negligible in the first part of the simulation.After approximately 30.00 s from the beginning of the simulation, the percentage error starts to increase and remains approximately constant until the end of the simulation.The maximum percentage difference between the computed current values is 23.54% at t = 7.36 • 10 3 s.We deem this value acceptable considering that it is found for polarization currents below 1.00 • 10 9 A. In addition, only the implicit solver has an adaptive time-step control.The code using the semi-implicit time-discretization scheme was written in Fortran 90, a language known for its high efficiency, whereas the other two using the implicit and the explicit scheme were written in Matlab.For this reason, Comparison between the polarization current computed using the implicit method and using the semi-implicit method.

FIGURE 6.
Comparison between the polarization current computed using the implicit method and using the explicit method.
we did not compare the computational time required by the semi-implicit scheme with the other ones.

B. EXPLICIT VS. IMPLICIT TIME DISCRETIZATION SCHEME
The results obtained from two simulations are shown in Fig. 6.The first simulation was carried out using the explicit integration scheme with C m = 0.8, see (30).The second simulation was conducted using Matlab's ode23tb function (implicit method).For both simulations, the flux density was discretized by means of a FOU scheme for the advective part and a centered finite difference formula for the diffusion part.At first glance the performance of the explicit scheme may seem superior to the one of the semi-implicit one show in Fig. 5.In reality, the maximum percentage error with respect to the implicit scheme was of 23.94%, which is larger than the one of the discussed semi-implicit scheme.However, this time the maximum error was found in an earlier stage of the simulated time, at t = 21.21s, where the current density is undergoing a rapid decrease.If the time step is reduced in the explicit scheme (using lower values of C m ), the maximum percentage error drops as expected.However, this comes at the cost of a larger computational time.The two simulations were timed using the Matlab timeit function.The one performed using the explicit time stepping took about 9.20 s, whereas the one with the implicit method based on ode23tb took only about 1.30 s, making it ∼ 7 times faster.Note that, for what concerns the implicit scheme results, the reported time measurement includes the time needed to complete the post processing phase to get the polarization current from the number densities.In reality, the most efficient solver depends on the time span covered by the simulation.Indeed, as shown in Fig. 7, the explicit integration scheme is the most computationally efficient for simulated time spans shorter than 1.00 • 10 4 s.Nevertheless, for the study of solid polymeric dielectric materials, it is often necessary to perform studies over time spans considerably longer than 1.00•10 4 s, where the implicit integration scheme is clearly the fastest.This difference can be quite relevant if many simulations need to be performed, e.g., in the case of an optimization problem aiming to find a best fit for physical parameters, such as those in [52], [53], [54], [55], and [56].Such optimization problems are of interest e.g. for diagnostic purposes.Even-though the identification of an appropriate parameter set could be non-trivial it could provide important insight about the condition of the material.In general, the results in Fig. 7 show that as the simulated time interval increases, the advantages of employing the implicit method become more pronounced.We conclude that using an implicit integration scheme is far more efficient than using the Euler explicit method when simulating solid polymeric dielectric materials over long time spans.The implicit scheme based on Matlab's ode23tb was used to obtain all the results that will be discussed in the rest of this work.

C. SOU/KL VS. FOU FLUX DISCRETIZATION SCHEME 1) RESULT COMPARISON
In this section we study the performance of two different schemes for the discretization of advective fluxes stemming  from the drift term in (2), i.e., the hyperbolic part of the drift-diffusion equation describing the number density evolution over time.We once again perform two simulations where all the physical parameters are set according to the ones in [18], as in the previous section, and only the scheme used for the drift term is varied.The first simulation employs a first-order upwind scheme (FOU), the second instead is performed using a second-order upwind scheme with the Koren flux limiter (SOU/KL).The results in Fig. 8 show the polarization current, computed by averaging the total current density across the domain at each time instant.The results are close, despite the higher order of accuracy of the SOU/KL with respect to the FOU.Throughout the simulation, the maximum percentage difference in the polarization currents computed with the two approaches was observed at t = 7.54• 10 3 s and has a value of 10.96%.For the same simulation, Fig. 9 shows the charge density distribution yielded by the two schemes over the whole domain at the last time instant of the simulation.The maximum absolute difference is observed  The four figures show only a fraction of the domain, closer to the left electrode.From these figures it is possible to observe a significant difference in the species number density computed using the FOU and the SOU/KL.The discrepancies are concentrated in the region between the left electrode and x = 4.00 • 10 −5 m.However, these differences have a small impact on the polarization current and on the charge density profile, at least for the physical conditions that we explored.On the other hand, the computational time variation between the two methods is not negligible.The simulation employing the FOU scheme was timed with the Matlab timeit function and took an average of 1.30 s.The simulation performed with the SOU/KL scheme was timed in the same way and took an average of 11.00 s.This corresponds to a ratio between the two computational  times greater than 8.In previous works by Le Roy and colleagues, such as [17], [18], [19], and [21], a third-order upwind scheme (QUICKEST) coupled with the ULTIMATE flux limiter was often employed to solve the transport equation.However, our analysis shows that the significant computational time increase due to the usage of a high-order scheme coupled with a flux limiter does not correspond to a relevant difference in the polarization current and charge density results.For this reason, we suggest to employ a highaccuracy scheme, such as the SOU/KL or the QUICKEST-ULTIMATE, only if the knowledge of the species number density with great detail is desired.In general, we suggest employing FOU when the polarization current or the charge density profile is computed.the species number density in the vicinity of the electrodes.The most accurate results are the ones obtained with the SOU/KL scheme, which is between first and second order accurate in space.On the other hand, the FOU scheme is only first order accurate.In order to better show the effect of this difference in accuracy a comparison between the two numerical schemes was performed considering simple test cases for which analytical solutions exist.As a first test case we solved the steady-state equation where A one-dimensional domain with length equal to L = 1.00 m, and discretized in 100 cells was considered.The diffusion coefficient K was set equal to 1.00 • 10 −1 m 2 s −1 and the velocity U was set to 5.00 m s −1 .Dirichlet conditions were considered for the number density n at both domain boundaries, in particular n W = 5.00 • 10 9 m −3 is the number density fixed at the western domain boundary and n E = 2.00 • 10 10 m −3 is the number density fixed at the eastern domain boundary.In this case an analytical solution exists for the distribution of the number density: In Fig. 14 the analytical solution is compared with the results yielded by the two schemes for the case presented above.In Fig. 15 the absolute value of the percentage error of the two schemes with respect to the analytical solution is shown.
From the two figures it is possible to see that the SOU/KL scheme produces results closer to the analytical solution with respect to the FOU scheme.The maximum percentage error  is always less than 3% for the SOU/KL scheme, while the FOU scheme reaches values of almost 10%.The second test case considered was a purely convective time dependent problem: = n U .
This time, a one-dimensional periodic domain with length equal to L = 1.00 m, and discretized in 1600 cells was considered.There is no need for boundary conditions being the domain periodic.The initial distribution for the number density consists of a triangular shape, a rectangular shape and a Gaussian.The velocity U was set equal to 1.00 m s −1 and a time span of 1.00 s was simulated employing the two flux schemes.In Fig. 16 the results at the end of the simulated time interval are shown.Again the results yielded by the SOU/KL are more accurate than the results obtained with the FOU.From Fig. 16 it is clearly visible the numerical diffusion associated with the FOU scheme [41].The test cases presented highlight the improved accuracy of the SOU/KL scheme with respect to the FOU scheme, in particular in those regions where steep gradients of the advected property occur.For this reason, in Fig. 10, 11, 12, and 13 the more accurate results are the one obtained using the SOU/KL scheme.

D. SATO'S FORMULA VS. TOTAL CURRENT DENSITY
A comparison between the polarization current obtained using two different techniques is shown in Fig. 17.The current is calculated starting from the output (i.e., the computed number densities over time) of the same simulation for the sake of consistency.In this way, it is possible to isolate the effect of the scheme employed to compute the polarization current.The continuous blue line in Fig. 17 shows the current directly computed from the number densities, the mobility, and the static electric field using Sato's equation.The orange dots show the current that was computed by averaging the total current density over the domain.The figure shows that the two methods are equivalent.The advantage of using Sato's equation is that the computation of the displacement current -which requires using data from at least two instants -is not needed.This results in a simpler implementation and a shorter computational time compared to the total current density method.

VI. CONCLUSION
In this work we compared several numerical techniques for the computation of charge transport and accumulation inside solid polymeric dielectric materials.The employed physical model is based on the well-known drift-diffusion equation for charged species and considers a constant mobility, a single level of deep traps and constant recombination, trapping, and detrapping coefficients.The model has a total of four types of carriers and the charge injection from the electrodes is evaluated by means of a Schottky emission law.We tested three different time discretization algorithms: a full explicit method, a full implicit method (implemented using Matlab's ode23tb function) and a semi-implicit technique.The full implicit method was identified as the most accurate and computationally efficient for long simulation times.Two different cell-centered finite volume schemes for the flux density discretization were compared: the first-order upwind scheme (FOU) and the second-order upwind scheme with the Koren flux limiter (SOU/KL).We observed that the polarization current and the charge density profiles are weakly dependent on the employed flux scheme.For this reason, since the code using the FOU was about 10 times faster than the (fully vectorized) implementation of the SOU/KL, we suggest to use the FOU unless a high accuracy in the species number density is required.
Concerning the computation of the polarization current, we highlighted the benefits of using Sato's equation instead of the more common approach based on the total current density.Indeed, we have shown that Sato's equation yields the same results, allowing for a simpler implementation and to avoid the computation of the displacement current at each time instant.

APPENDIX SATO'S FORMULA DERIVATION
A complete step-by-step derivation of Sato's formula is presented in this section.A generic 3D domain with two electrodes is considered.A generic charge distribution is let free to evolve in time inside the domain.The total electric field inside the domain may be expressed as the sum of two terms: the static (or Laplacian) electric field E s , which is dependent only on the applied voltage, and the electric field due to the charge distribution inside the domain, E ρ : We can write: where D ρ is the displacement field due to the charge distribution and D s is the static displacement field.We can express E ρ as the gradient of an electric potential ϕ ρ : The assumption that a constant voltage is applied between the electrodes is made, and thus: An extended version of Sato's equation where time-dependent electric fields are considered can be found in [57].The generic continuity equation is written: In (70), S = E × H is the Poyinting vector, B is the magnetic flux density, and H is the magnetic field.
If we consider a constant B inside the domain and a linear and isotropic dielectric material, with constant dielectric permittivity ε, (70) can be rewritten as: Equation ( 71) is integrated throughout the domain, producing the following.
The divergence theorem is applied to the first member of (72): Assuming that the normal component of the displacement current is equal to zero at the domain boundary (without the need to include the electrodes), we obtain: Using ( 65) and (69) it is possible to write: Employing (68) we obtain: We recall the identity valid for a generic vector V and a scalar f : Applying (77) to the second integral on the right-hand side of (76) yields: The divergence theorem is then applied, obtaining: 12558 VOLUME 12, 2024 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The surface integral in (79) is equal to zero since the electric potential ϕ ρ is zero at the electrodes and the normal component of the conduction current may be assumed to be zero at the remaining portion of the boundary surface.ϕ ρ is zero since at the electrodes the electric potential is imposed and therefore the contribution due to internal charges must be zero in those regions.In addition, we recall that and using (66) we can rewrite (79) as: Applying (77) to the second integral in the right hand side of (81) and using (68), we get: The divergence theorem is applied, yielding: The surface integral in (83) is equal to zero since ϕ ρ is zero at the electrodes and the normal component of the displacement current is zero elsewhere at the boundary surface.Using (65) and recalling that the material is linear and isotropic, we get the following.
Applying (68) we get: Employing (77) yields: Recalling (67) and applying the divergence theorem, (86) can be rewritten as: The surface integral in (87) is equal to zero if we assume that the normal component of the static displacement field (D s • n) is zero at the domain boundary, without necessarily including the electrodes, where the time derivative of ϕ ρ is equal to zero.
In this case, we obtain: Considering the following expression for the conduction current density we obtain the same expression as (58).

FIGURE 1 .
FIGURE 1. Koren flux limiter function for different values of r , the ratio between advected property gradients, defined in (26) and (27).

FIGURE 3 .
FIGURE 3. Schematic representation of the algorithm used for the semi-implicit method.

FIGURE 4 .
FIGURE 4. Structure of a typical simulation.

FIGURE 5 .
FIGURE 5.Comparison between the polarization current computed using the implicit method and using the semi-implicit method.

FIGURE 7 .
FIGURE 7. Wall clock time versus simulation time for the explicit and implicit time schemes.

FIGURE 8 .
FIGURE 8. Comparison between the polarization current computed using the FOU and the SOU/KL schemes.

FIGURE 9 .
FIGURE 9.Comparison between the charge density computed using the FOU and the SOU/KL schemes.

FIGURE 10 .
FIGURE 10.Comparison between free electron number density computed using the FOU and the SOU/KL schemes in a domain fraction close to the left electrode.

FIGURE 11 .
FIGURE 11.Comparison between free hole number density computed using the FOU and the SOU/KL schemes in a domain fraction close to the left electrode.

FIGURE 12 .
FIGURE 12.Comparison between trapped electron number density computed using the FOU and the SOU/KL schemes in a domain fraction close to the left electrode.

FIGURE 13 .
FIGURE 13.Comparison between trapped hole number density computed using the FOU and the SOU/KL schemes in a domain fraction close to the left electrode.

2 )
ACCURACY COMPARISON USING ANALYTICAL SOLUTIONSFrom Fig.10, 11, 12, and 13 it is possible to see that the two flux schemes yield greatly different results regarding

FIGURE 14 .
FIGURE 14.Comparison between the results yielded by the FOU and SOU/KL schemes and the analytical solution of a stedy-state problem.

FIGURE 15 .
FIGURE 15.Absolute value of the percentage error of the two schemes with respect to the analytical solution for the steady-state case.

FIGURE 16 .
FIGURE 16.Comparison between the results yielded by the SOU/KL and FOU schemes and the analytical solution after 1.00 s.

FIGURE 17 .
FIGURE 17.Comparison between the polarization current computed using Sato's formula and the total current density.
S 2 n h µ n e t − S 3 n h µ n e µ ;