Multiscale Modeling and Simulation of Organic Solar Cells
Abstract.
In this article, we continue our mathematical study of organic solar cells (OSCs) and propose a twoscale (micro and macroscale) model of heterojunction OSCs with interface geometries characterized by an arbitrarily complex morphology. The microscale model consists of a system of partial and ordinary differential equations in an heterogeneous domain, that provides a full description of excitation/transport phenomena occurring in the bulk regions and dissociation/recombination processes occurring in a thin material slab across the interface. The macroscale model is obtained by a microtomacro scale transition that consists of averaging the mass balance equations in the normal direction across the interface thickness, giving rise to nonlinear transmission conditions that are parametrized by the interfacial width. These conditions account in a lumped manner for the volumetric dissociation/recombination phenomena occurring in the thin slab and depend locally on the electric field magnitude and orientation. Using the macroscale model in two spatial dimensions, device structures with complex interface morphologies, for which existing data are available, are numerically investigated showing that, if the electric field orientation relative to the interface is taken into due account, the device performance is determined not only by the total interface length but also by its shape.
Keywords: Organic solar cell; nonlinear reactiondiffusion system with electrostatic convection; scale transition; multiscale analysis; numerical simulation; finite element method.
1. Introduction and Motivation
Research on photovoltaic energy conversion has recently received great impulse due to the growing demand for low carbon dioxide emission energy sources. In particular, the high manufacturing cost of crystalline silicon and the latest advancements on semiconducting polymer design and synthesis in recent years have directed the attention of the scientific community towards Organic Solar Cells (OSCs), i.e. solar cells based on organic materials [12, 20, 21, 32, 33, 38, 42], especially because of the very limited thermal budget required for the production of such materials and of their amenability to be deposited on large areas, which is fundamental in light harvesting applications. One of the main peculiarities of OSCs is that most physical phenomena that are critical for charge photogeneration occur at the interface between the two materials that constitute the active layer of such devices. In order to increase cell efficiencies, currently of the order of about 10% [22], the optimization of the morphology of such interface is considered by device designers to be an issue at least as important as the optimization of the donor and acceptor optoelectronic characteristics [40].
For this reason, in this article we continue our mathematical study of organic photovoltaic device models started in [17], and we focus on the accurate and computationally efficient modeling of the main dissociation/recombination processes occurring in a thin material slab across the material interface and evaluating their impact on device photoconversion performance. With this aim, we consider a twoscale approach to OSC simulation that is intermediate between a continuum model [6] and a full microscopic model [30, 11], and represents an extension to the case of arbitrary interface geometries of the onedimensional model for bilayer OSC devices proposed in [4].
The approach is based on the introduction of two distinct levels of description of the physical system at hand, a micro and a macro scale, and of two corresponding mathematical models based on classical mass balance conservation laws. At the microscale, a system of PDEs in a heterogeneous domain provides a full description of the excitation/transport phenomena occurring in the bulk regions and of the dissociation/recombination processes occurring in a thin material slab across the interface. The numerical treatment of the microscale model presents several difficulties related to the wide difference in size between the bulk regions and the interfacial width . As a matter of fact, as polaron dissociation is assumed to occur in the first layer of polymer chains on either side of the interface, the length scale can be taken to be that of the average separation between polymer chains which is typically more than two orders of magnitude smaller than the size of the bulk regions [4, 8, 44].
Therefore, to obtain a computationally efficient model, we carry out a microtomacro scale transition that somewhat resembles modelreduction techniques used for porous media with thin fractures [31], for reaction problems with moving reaction fronts [29] and for electrochemical transport across biological membranes [35], and relies on a systematic averaging of the mass balance equations in the normal direction across the interface thickness. The resulting macroscale model is a system of incompletely parabolic PDEs to describe mass transport in the materials, nonlinearly coupled with ODEs and transmission conditions localized at the heterojunction parametrized by the interfacial width . These conditions account in a lumped manner for the volumetric dissociation/recombination phenomena occurring in the interfacial thin slab. The fact that in the macroscale model the interface is reduced to a zero width surface is further exploited to account for the local dependence of the polaron dissociation rate on the electric field orientation, which is the main advantage –together with the computational cost reduction– of our approach, as compared to previous multidimensional models [8, 44, 26].
An outline of the article is as follows. In Sect. 2 we illustrate the sequence of physical phenomena that lead from photon absorbtion to current harvesting in an OSC. Sect. 3 is devoted to characterizing the mathematical model of an OSC. In Sect. 3.1 we describe the geometrical heterogeneous structure of the device, while in Sect. 3.2 we introduce the basic modeling assumptions on the dependent variables of the problem. Then, in Sect. 3.3 we present the microscale PDE/ODE model system with the initial and boundary conditions, while in Sect. 3.4 we describe in detail the scale transition procedure that leads from the microscale model to the macroscale equation system. We complete the mathematical picture of a bilayer OSC by illustrating in Sect. 3.5 a novel model that we have devised for including the dependence of the polaron dissociation rate constant on the local electric field and on the morphology of the material interface. In Sect. 4 we briefly comment about the numerical methods used for discretizing the macroscale model, while Sect. 5 is devoted to presenting and discussing numerical results. In particular, in Sect. 5.1 we successfully perform the validation of the accuracy of the macroscale model with respect to the microscale model through the numerical simulation of a one–dimensional OSC under different working conditions. Extensive simulations of two–dimensional OSC structures are instead reported in Sects. 5.2 to 5.4 in order to both validate the proposed macroscale model with respect to previously available numerical results and to analyze its effectiveness in the study of complex interface morphologies. Finally, in Sect. 6 we draw some conclusions and sketch possible directions for further research in the area of OSC modeling and simulation.
2. Basic Principles of Photocurrent Generation in OSCs
In this section, we describe the basic principles of photocurrent generation in OSCs only to the extent strictly needed for understanding the naming conventions adopted in the following sections. For a more thorough introduction to the subject, we refer the interested reader to [20, 38, 32]. The typical structure of an OSC is constituted by a thin film a crosssection of which is schematically represented in Fig. 1.
The photoactive layer of the device consists of two materials, one with higher electron affinity (the “acceptor”, for example F8BT, P3HT) and one with lower electron affinity (the “donor”, for example PFB, PCBM), sandwiched between two electrodes, one of which is transparent to allow light to enter the photoactive layer while the other is reflecting in order to increase the light path through the device.
The sequence of physical phenomena that leads from photon absorption to current harvesting at the device contacts is represented in Fig. 2.
Absorption of a photon in either material produces an electronhole pair, usually referred to as an exciton whose binding energy is of the order of about . Excitons may diffuse through the device until they either recombine or reach the interface between the donor and acceptor phases. If this latter event occurs, the exciton may get trapped at the interface in such a way that its electron component lays in the high electron affinity region while the hole component lays in the low electron affinity region. Such a trapped excited state is referred to as a polaron pair or geminate pair [4, 44, 36, 34] and has a lower binding energy compared to that of the exciton state, as the Coulomb attraction between the electron and hole is reduced by the chemical potential drop between the two materials. The polaron binding force may be overcome by the electric field induced by the small builtin voltage between the metallic contacts thus leading to the formation of two independent charged particles (one electron and one hole), otherwise the polaron pair may return to the untrapped exciton state or recombine. Free charge carriers move by drift and diffusion mechanisms and, unless they are captured along their path by the coulombic attraction of an oppositely charged particle and recombine at the interface to form a new polaron pair, they eventually reach the contacts thus producing a measurable external current.
3. Mathematical Model
In this section, we propose a PDE/ODE model of photoconversion and charge transport mechanisms in an OSC. The model relies on a twoscale approach that is based on the introduction of two distinct levels of description of the physical system, namely, a micro and a macro scale, and of two corresponding mathematical equation systems based on classical mass balance conservation laws. The construction of the model proceeds through four steps. In Sect. 3.1, we describe the geometrical and heterogeneous structure of the device which consists of two bulk regions (the acceptor and donor phases) separated by an interface region of (finite) thickness , while in Sect. 3.2, we introduce the basic modeling assumptions on the dependent variables of the problem. In Sect. 3.3, we introduce the microscale PDE/ODE model system of conservation laws that governs transport of the various species throughout the device, together with its initial and boundary conditions and the generation/recombination mechanisms that occur in each subdomain of the heterogeneous device. In Sect. 3.4, we describe in detail the scale transition procedure that is applied to the microscale model in order to obtain the macroscale equation system. This latter system basically consists of the same equations as in the microscale model, but satisfied in the separate acceptor and donor phases, coupled through a set of flux transmission conditions across the material interface that synthetize in a “lumped” manner the dissociation and recombination mechanisms that actually occur in the thin volumetric slab of width surrounding the interface itself. The resulting macroscale model system is a compromise between a continuum model and a full microscopic model, and represents a consistent mathematical rationale and generalization of the various models proposed in [4, 43, 44]. We conclude our mathematical picture of the OSC by illustrating in Sect. 3.5 a novel model of the polaron dissociation rate properly devised for including the dependence on the local electric field and on the morphology of the material interface.
3.1. Geometry of the Heterogeneous Computational Domain
A schematic 3D picture of the OSC is illustrated in Fig. 3(a). The device structure is a parallelepipedshaped open subset of divided into two open disjoint subregions, (acceptor) and (donor), separated by a regular oriented surface [14] on which, for each , we can define a unit normal vector directed from into . The top and bottom surfaces of the structure are mathematical representations of the cell electrodes, cathode and anode, denoted as and , respectively, in such a way that and (see Fig. 3(b)). We also denote by the unit outward normal vector over the cell boundary .
Following [4, 43, 44], it is convenient, for modeling purposes, to associate with the interface the subregion depicted in Fig. 4(a) and defined as follows. For each point , let be the “thickness” associated with . Then, set
(1) 
The subregion is thus a 3D thin layer of thickness surrounding which represents the device volumetric portion where the dissociation and recombination mechanisms of Sect. 2 are assumed to occur. It is worth noting that the width is, in general, an unknown of the physical problem. As such, it depends on and , it may assume different values in the two material phases of the photoactive layer and might in principle locally depend on the electric field. According to the data provided in [4, 43, 44], we assume for simplicity to be a constant quantity. Based on the definition (1), we can introduce the two portions and (see Fig. 4(b)). Consistently, we also introduce the boundary portions and and set , in such a way that , and , where . Notice that, unlike , the surfaces and can be regarded as “mathematical” interfaces.
3.2. Modeling Assumptions
Let us denote by , , and the volumetric densities of (singlet) excitons, polaron pairs, electrons and holes, respectively, and by , , and the associated particle fluxes. Based on the physical working principles of an OSC illustrated in Sect. 2, on the heterogeneous geometrical decomposition of the device introduced in Sect. 3.1 and the extensive numerical simulations reported in [8], we make the following modeling and geometrical assumptions:
 A.1:

excitons can be generated at any position in the cell, so that is a nonnegative function over all the cell domain ;
 A.2:

electrons (holes) are not able to penetrate the donor (acceptor) material beyond the interface layer , so that electrons (holes) are nonnegative functions over (), and are identically equal to zero in ();
 A.3:

polarons are trapped and immobile in the interface region , so that is a nonnegative function over and identically equal to zero in ;
 A.4:

the OSC is in the “off” state at (that is, before illumination), so that the initial condition for all the involved densities is for all ;
 A.5:

the geometry of the device is an infinite periodic repetition of the computational domain of Fig. 3, so that periodic boundary conditions are enforced for all variables on the lateral boundary of .
3.3. Microscale Model
In this section, we illustrate the microscale model we advocate in this work to be a mathematical representation of the functioning of an OSC.
We take excitons to obey:  
(2a)  
where we use assumptions A.1 and A.3 to define  
(2b)  
and  
(2c)  
The superscripts and represent the fact that the corresponding volumetric production terms are defined in the bulk and interface regions, respectively. The term denotes the rate at which excitons are generated by photon absorption and is henceforth assumed to be a nonnegative given function of time and position while is the exciton lifetime in the bulk materials. In the interface region additional dissociation and recombination mechanisms are taken into account and and represent the rate constants for the transition of excitons to the polaron state and that of polarons back to the exciton state, respectively. In particular denotes the total rate of polaron recombination events and the fraction of such events which produce a singlet exciton. As excitons have zero net charge, their flux is driven by diffusion forces only, i.e. the flux density may be expressed as  
(2d)  
being the exciton diffusion coefficient. At the contacts we assume perfect exciton quenching [37] so that  
(2e) 
Because of assumption A.2, the following equations hold for electrons:  
(3a)  
and for holes:  
(3b)  
where the term is defined as  
(3c)  
Notice that is identically zero in the bulk region as electrons and holes can only recombine with each other, so no recombination occurs where either of the two species is missing. In the interface region both electrons ad holes exist so the terms in take into account for polaron pair dissociation with rate constant (see Sect. 3.5 for the model) and bimolecular recombination with rate constant . As electrons and holes each bear a nonzero net charge, their flux is driven by both diffusion and electric drift forces [25], therefore:  
(3d)  
and  
(3e)  
where is the electric field while , and , are the diffusion coefficient and mobility for electrons and holes, respectively.  
Because of assumption A.2, the following boundary conditions hold at the artificial interfaces separating the donor and acceptor bulk phases from the thin slab region :  
(3f)  
and  
(3g)  
At the contacts we impose the same Robintype boundary conditions as described in [10, 17]:  
(3h)  
and  
(3i)  
where , , , , are nonnegative coefficients. 
The electric field in (3d) and (3e) is connected to the electric potential by the quasistatic approximation  
(4a)  
and satisfies the Poisson equation  
(4b)  
where is the space charge density in the device. Using assumption A.2, the piecewise smooth definition of turns out to be:  
(4c)  
denoting the quantum of charge. The electric permittivity is equal to , and being the relative material and vacuum permittivities, respectively, with in the acceptor phase and in the donor phase, so that may be discontinuous across the interface . Dirichlet boundary conditions for the electric potential are set at the contacts and , as follows  
(4d)  
and  
(4e)  
where is the builtin voltage of the cell, and are the contact metal work functions while is the externally applied voltage. 
Because of assumption A.3, the flux is identically equal to zero in all and for all , and polarons satisfy the following ODE in the interface region  
(5a)  
while their density is identically zero in the bulk  
(5b) 
3.4. MicrotoMacro Scale Transition
The microscale model of a bilayer OSC described in Sect. 3.3 can be subdivided into three distinct groups of equations:
 1):

parabolic PDEs enforcing mass conservation of excitons, electrons and holes;
 2):

an ODE describing the kinetics of photogenerated polaron pairs;
 3):

an elliptic constraint enforcing Gauss theorem in differential form to be satisfied at each time throughout the whole cell domain.
The markedly spatially heterogeneous nature of the problem may be quite impractical for numerical simulation, in particular when devices with complex interface morphology in multiple spatial dimensions are considered. For this reason, in this section we propose a scale transition procedure which allows us to derive a macroscale model that is more amenable to numerical treatment. Other examples of multiscale mathematical approaches that are based on the scale separation concept and scale transition can be found in [31, 29, 35].
To construct our multiscale model of an OSC, we abandon the perspective focused at the nanoscopic characteristic level adopted so far, and prefer to look at the cell from a “larger” distance. By doing so, necessarily, we loose control of the details (i.e, we cannot distinguish the region from the two bulk regions and ), but, at the same time, we gain the advantage of not needing to resolve the interfacial bulk region across . The resulting macroscale problem is thus posed in the partitioned domain (as a matter of fact, we are still able to neatly distinguish the interface separating the two material phases!) without including the interfacial production terms in the mass balance and kinetics equations 1) and 2) introduced above.
Of course, we cannot simply limit ourselves to neglecting these latter terms, rather, we do need to incorporate their effects, in the macroscale model, in an alternative way. For this, the simplest approach to microtomacro scale transition consists of replacing , at each point of and for each time level, with its average across the thickness of in the normal direction, and then, of using as a source term for suitable flux transmission conditions, to be enforced on the interface in the case of mass balance equations. In the case of polaron pair equation, the averaging procedure automatically transforms the volumetric kinetics balance within into a surface kinetics balance over . In any case, the scale transition results in the introduction of suitable interfacial terms that replace in a “lumped” manner the volumetric dissociation/generation phenomena microscopically occurring in . Having characterized the averaging procedure for equations 1) and 2), the (macroscale) differential Gauss theorem 3) remains automatically (formally) unchanged and is expressed in terms of the (macroscale) space charge density as in Eqns. (4).
3.4.1. Derivation of the Macroscale Equations
The macroscale model for excitons reads:  
(6a)  
with  
(6b)  
and subject to the interface conditions:  
(6c)  
where denotes for any function the jump of across the interface , and being the traces on of the restrictions of from and , respectively. The continuity of at the interface is a requirement consistent with the elliptic regularity of both microscale and macroscale problems.  
The interfacial source term is defined as  
(6d)  
In the above relation, is the (single–valued) trace of over , while is the areal density of the bonded pairs, defined as  
(6e)  
and the midpoint quadrature rule is used for approximating the third integral in (6d). The macroscale model for excitons is completed by the constitutive relation (2d) for exciton flux density and by the perfect exciton quenching boundary conditions (2e). 
The macroscale model for electrons reads:  
(7a)  
subject to the interface/boundary condition  
(7b)  
The interfacial source term is defined as  
(7c)  
where definition (6e) is used in the first integral while the midpoint quadrature rule is again used to approximate the third integral in (7c). The macroscale model for electrons is completed by the constitutive relation (3d) for electron flux density and by the Robintype boundary condition (3h). 
Proceeding in a completely analogous manner as done with electrons, the macroscale model for holes reads:  
(8a)  
subject to the interface/boundary condition  
(8b)  
The macroscale model for holes is completed by the constitutive relation (3e) for hole flux density and by the Robintype boundary condition (3i).  
The conditions (7b) and (8b) assume an interesting physical meaning upon introducing the electron and hole current densities, defined respectively as and , and the total (conduction) current density . Recalling that () in (), we have:  
from which we get  
(8c)  
that expresses the property of current conservation across the interface . 
Integration of (5a) across the interface thickness yields the following macroscale model for the areal density of polaron pairs  
(9a)  
where  
(9b) 
The macroscale model for the differential Gauss theorem is expressed by the following Poisson problem in heterogeneous form:  
(10a)  
with  
(10b)  
and subject to the interface conditions:  
(10c)  
Two remarks are in order with system (10). First, we notice that the Gauss theorem in differential form (10a) looks formally identical to the corresponding microscale formulation (4b), the difference between the two methodologies being in the definition of the space charge density (compare (4c) with (10b)). Second, the transmission conditions (10c) express the physical fact that the normal component of the electric displacement vector and the electric potential do not experience any discontinuity at the material interface, as is the case of the microscale formulation. 
3.4.2. Summary of the Macroscale Model
For sake of convenience, we summarize below the macroscale model of an OSC written in primal form:
(11a)  
(11b)  
(11c)  
(11d)  
(11e) 
System (11) is completed by periodic boundary conditions on , as stated in assumption A.5. For the physical models of the coefficients in system (11) we refer to [4, 19, 23], except for the description of the polaron dissociation rate constant which is addressed in detail in Sect. 3.5. In particular, for the carrier mobilities, we neglect the effect of energetic disorder, so that they can be assumed to depend only on the electric field magnitude, according to the PooleFrenkel model. As for diffusivities, in the computations of Sect. 5, Einstein relations
(12) 
are assumed to hold, although the proposed multiscale formulation remains unchanged if such an assumption is removed. In (12), is Boltzmann’s constant and is the absolute temperature. Finally, for the bimolecular recombination rate constant a Langevintype relation is used [4].
3.5. Model for the Polaron Dissociation Rate
Numerical simulations as those reported in Sect. 5 show that the polaron dissociation rate has a significant impact on the cell photoconversion efficiency, for this reason we devote this entire section to modeling the dependence of on the electric field and on the morphology of the material interface. A commonly used polaron dissociation rate model is the BraunOnsager model [7] which is derived assuming the OSC bulk to be a homogeneous medium and takes into account only the magnitude of the electric field. In [4] the authors propose a model for tailored for bilayer devices which is derived by performing an average over the admissible range of the escape angle relative to the electric field direction. In this latter model, the electric field is assumed to be always directed orthogonally to the interface consistently with the planar geometry of the device considered therein. The authors of [44] apply the dissociation rate model of [4] to more complex geometries by performing an average along the interface of the field component normal to the contacts which amounts to neglecting the effect of local electric field orientation.
To construct a novel model which also takes into account this latter effect we repeat the derivation of [4] with two differences. The first difference is that of removing the assumption that the field is normal to the interface. The second difference is that of considering a limited range of admissible escape directions to account for the physical fact that polaron pairs tend to be aligned with the gradient of the electron affinity due to the different materials in the two device subregions.
Referring to Fig. 5 for the geometrical notation, we let
(13) 
where is the zerofield dissociation rate constant, is the escape direction of the electron part of the polaron at the point , is a nonnegative weight representing the probability distribution of admissible escape directions, and such that , and is an enhancement/suppression factor given by the PooleFrenkel formula:
(14) 
having set . The product can be expressed in terms of the normal component and the tangential component of the electric field as
To specify an expression for we assume an escape direction to be admissible only when the angle it forms with respect to the normal unit vector is not too large. Indicating by the maximum admissible value for and allowing all admissible values to be equally likely, we obtain:
Two limits are of particular interest, and .
In the first case, Eq. (13) can be checked to yield
(15) 
This corresponds to assuming that all geminate pairs are exactly aligned with the interface normal unit vector, thus neglecting any possible variability in their orientation due to, e.g. interface surface roughness and/or thermal vibrations.
In the second case, Eq. (13) becomes
(16) 
which, in the special case where , coincides with Eqs. (17)(21) of [4]. Notice that if the choice may overestimate the effective dissociation rate as it corresponds to completely neglecting the alignement of the geminate pairs with the electron affinity gradient. This is observed to give rise to nonphysical effects as shown by the simulations of Sect. 5.2. Therefore, for practical purposes, the quantity should be used as a fitting parameter to be calibrated on experimental data.
Fig. 6 shows the dissociation rate constant (normalized to ) computed by model (15) (left) and (16) (right) for several values of the angle between and and having set and .
We notice that the dissociation rate computed by model (16) has a significantly smaller range of variability than predicted by model (15). A possible explanation to this difference is related to the smoothing operated by the integral in (16). The higher variability of the dissociation rate translates into an higher sensitivity of model (15) to the inclination of the electric field with respect to the interface normal as will be further discussed in the numerical results section when commenting Fig. 14(b). A discussion of the impact of (15) and (16) on the model predictions will be carried out in Sect. 5.2.
4. Numerical Approximation
In this section we describe the numerical techniques used to solve the mathematical models introduced in Sects. 3.3 and 3.4. The full details of the discrete system of linear algebraic equations resulting from problem approximation are postponed to A. As for the simulation of the model in the transient regime carried out in Sect. 5.1, we have adapted to the case at hand the numerical method described in [17] based on Rothe’s method and on the use of adaptive Backward Differentiation Formulas (BDF). In the steadystate simulations illustrated in Sects. 5.2, 5.3 and 5.4, all partial derivatives with respect to time have been dropped out in system (11) in such a way that Eq. (11b) reduces to an algebraic constraint.
The numerical strategy adopted in the present paper is basically composed of three steps:

Linearization

Spatial discretization

Solution of the linear algebraic system
Step (1)
For model linearization, we adopt a quasiNewton approach similar to
that used in [17], where, in the computation of Jacobian matrix
entries, the dependence of mobilities and polaron pair dissociation rate on the solution
is neglected.
Step (2)
Similarly to [17], for the spatial discretization of the sequence of linear systems of PDEs stemming from Step (1) we adopt the GalerkinFinite Element Method (GFEM) stabilized by means on an Exponential Fitting technique [3, 18, 45, 28]
in order to deal with possibly dominating drift terms in the continuity equations.
A peculiarity of the heterojunction
model (11) as compared to the homogenized model of [17] is the presence of nontrivial interface conditions at the donoracceptor interface, which are taken care of by means of the substructuring techniques described, e.g.,
in [39, 24] which turn out to be of straightforward implementation in the adopted GFEM method.
Step (3)
To solve the linear algebraic systems arising from problem
discretization, we employ the Unsymmetric Multi Frontal method
implemented in the UMFPACK library [15]
as on current hardware architectures memory constraints
are not the main limiting factor and the use of
a direct sparse solver has the advantage of being
more robust than iterative approaches
with respect to coefficient matrix conditioning
5. Simulation Results
In this section we carry out an extensive computational study of the micro and macroscale models introduced in Sects. 3.3 and 3.4. In Sect. 5.1, onedimensional transient simulations under different working conditions are carried out to verify the accuracy of the macroscale model with respect to the microscale system. In Sects. 5.2, 5.3 and 5.4, computations are performed in steadystate conditions, in order to validate the macroscale model by comparison with available results in the literature. The numerical schemes of Sect. 4 have been implemented in Octave using the OctaveForge package bim [16] for matrix assembly.
5.1. Numerical Validation of the Accuracy of the Macroscale Model
In Sects. 3.4 and 3.3, we have illustrated two models of the operating principles of bilayer OSCs at two increasing levels of detail, corresponding to the macro and micro scales, respectively. The two modeling descriptions are expected to provide a correspondingly more refined level of quality in the representation of the principal physical phenomena that govern the functioning of an OSC, at the price, however, of a substantial increase in implementation complexity and computational effort, especially in the case of multidimensional simulations (mesh generation of complex interface morphologies, solution of large algebraic systems with possibly badlybalanced matrices). The natural question that arises at this point of the discussion is whether the macroscale formulation of Sect. 3.4.2 is capable of returning an output picture of the performance of a bilayer OSC with sufficient accuracy compared to that of the microscale formulation of Sect. 3.3. By construction of the two models, the extent by which the word “accuracy” is mathematically identified cannot certainly refer to a pointwise comparison between micro and macroscale solutions (they will certainly look different!), rather, it should be concerned with average quantities that best represent the overall performance of the device. In this respect, the verification test we are going to carry out later on, is a check of the total current per unit area predicted by the micro and macro formulations, where , being either of the two contacts or . The choice of for model validation is due to the fact that the total output current density is an easily accessible quantity in experiments, and thus represents the most significant parameter for assessing the photoconversion performance of a solar cell [32].
For sake of computational simplicity, we consider a biplanar OSC, so that the resulting spatial geometrical description can be reduced to a one–dimensional model. The total length of the device is equal to nm, with the two regions and occupying each one half of the cell. All model coefficients are assumed to be constant quantities, and their values are listed in Tab. 1.
Parameter  Symbol  Numerical value 

Acceptor relative dielectric constant  2.5  
Donor relative dielectric constant  2.5  
Builtin voltage  0.6V  
Temperature  298K  
Electron mobility  
Hole mobility  
Exciton diffusion coefficient  
Exciton lifetime  
Exciton dissociation time  
Polaron pair recombination rate constant  
Singlet exciton recombination fraction  
Polaron pair dissociation rate constant ( V)  
Polaron pair dissociation rate constant ( V)  
Bimolecular recombination rate constant  
Boundary condition parameters for electrons  
Boundary condition parameters for holes  
We start by simulating the cell turnon transient at short circuit condition (), which corresponds to computing the device response to an abrupt variation at in the photon absorption rate from zero to . The simulation time interval is taken wide enough for the device to reach stationary conditions. In Fig. 7(a) the relative discrepancy between the stationary values of computed with the two methods is reported for several values of the interface width parameter . Results allow us to conclude that in the physically relevant range of variation of , the relative discrepancy between the micro and macroscale models remains consistently below 10% and that, as expected, the predictions of the two models tend to become undistinguishable as tends to zero. In Fig. 7(b) we set and we show the time evolution of the total current per unit area in the two biasing conditions and this latter being the “flat band” voltage. Accordingly to Fig. 7(a), in both regimes the curves almost coincide over the whole simulation time interval.
5.2. Model Validation through Comparison with Existing Simulation Data
In this section, we aim to compare the predictions of our macroscale model to those of [43, 44] and to investigate the impact of the model for proposed in Sect. 3.5 on the simulated device performance. We consider the same device as in [43, 44, 13] where the acceptor and donor materials are F8BT and PFB, respectively. The values of the model parameters are listed in Table 2.
The device morphology, shown in Fig. 8, is an interpenetrating rodshaped structure of donor and acceptor materials with , , and . Throughout this section, we denote by the direction between the two electrodes and .
In [43, 44] an optical model has been used to determine the exciton generation term . Here, instead, we follow a simpler approach by considering to be constant in the entire device structure and equal to the value obtained averaging the result in [43, 44]. The choice corresponds to enforcing Dirichlet boundary conditions for the carrier densities at the contacts and amounts to neglecting the dependence of charge injection on the electric field and assuming an infinite recombination rate at the contacts.
Fig. 9(a) shows the current densityvoltage characteristics in the case of an exciton generation rate .
Parameter  Symbol  Numerical value 

Acceptor relative dielectric constant  4  
Donor relative dielectric constant  4  
Temperature  298K  
PooleFrenkel mobility model parameters for electrons [4]  
PooleFrenkel mobility model parameters for holes [4]  
Exciton diffusion coefficient  
Exciton lifetime  
Exciton dissociation time  
Polaron pair recombination rate constant  
Singlet exciton recombination fraction  
Polaron pair zerofield dissociation rate constant  
Interface width  
Boundary condition parameters for electrons [17]  
Boundary condition parameters for holes [17]  
The three curves correspond to the use of three different expressions for the polaron dissociation rate constant , identified as follows: (A) the model proposed in [43, 44] with as the driving parameter for polaron pair dissociation (solid line); (B) the model (16) (dashdotted line); (C) the model (15) (dashed line). The result computed using model (A) is in excellent agreement with that of Fig. 7(right) in [44] despite the above mentioned modeling differences. Model (A) does not account for the orientation of the electric field with respect to the donoracceptor interface and is expected to overestimate dissociation in the case where . This is confirmed by the curve for model (B). As a matter of fact, in this case all the dissociation directions are assumed to be equally likely and the computed output current density before flatband condition occurs () is smaller than predicted by the solid line curve. For the computed currentvoltage characteristic exhibits a nonmonotonic behavior. This latter behavior is not observed in any of the experimental measurements we are aware of, and is most probably to be ascribed to a too important contribution of the tangential component of the electric field that leads (16) to overestimate polaron dissociation at the material interface. If, instead, model (C) is used, the obtained output current density characteristics is the dashed line in Fig. 9(a). We observe a smoother trend than in previous cases for all applied voltages, and close to short circuit condition we note that the current density is further reduced since dissociation is assumed to occur only in the normal direction and on a significant portion of the interface is almost vanishing. In all the considered cases, the nonsmooth behavior at flat band condition () is to be ascribed to the discontinuity of at in (14).
Fig. 9(b) shows the results of the same analysis as above in the case of an exciton generation rate . The shape of the characteristics is very similar to those with low light up to a scaling factor of about 100, this suggesting a linearity between the output current density and the illumination intensity. Notice the absence of the bump for in the case of model (B). This is a consequence of the increased magnitude of the charge carrier densities compared to the previously considered illumination that in turn determines stronger Coulomb attraction forces and hence more recombination phenomena. With reduced attractions, instead, charge carriers have more chances to escape from the interface following concentration gradients.
In Fig. 10 we show the charge carrier densities in a device with geometrical data set to , , and , at short circuit condition with exciton generation rate . We first observe that computed charge carrier distributions in Fig. 10(left) are in very good agreement with those of Fig. 3(i) in [43]
and show the same peaks close to the vertical sides of the donoracceptor interface. It is interesting to notice that the total number of holes in the donor material is higher than the number of electrons in the acceptor material because of the significantly different values of their respective mobilities. Negative charges can move through the device faster to be finally extracted at the cathode so that an overall positive charge buildsup in the device. The charge densities computed using models (B) and (C) exhibit a qualitatively similar profile with a gradual reduction of the magnitude compared to the result of model (A). This behavior is completely consistent with the previous analysis of the currentvoltage characteristics predicted by the three models of .
We conclude this preliminary validation of model (11) by illustrating in Fig. 11 the open circuit voltage and short circuit current density of a device with the same characteristics as in the previous set of simulations for values of exciton generation rate in the range from to . Fig. 11(a) is in excellent agreement with Fig. 6(right) of [44], and indicates that models (A), (B) and (C) predict a linear behavior of with respect to the logarithm of the exciton generation rate, as already pointed out in [43, 44]. Fig. 11(b) illustrates the current density that can be extracted from the device at short circuit condition. The logscale plot indicates that increases linearly in a wide range of illumination regimes until values of the order of . With more intense irradiation the increase becomes sublinear, suggesting that saturation of the device occurs due to more relevant excitonic and electronhole recombination phenomena which in turn are a consequence of the increased densities.
5.3. The Role of Interface Morphology
In this section, we aim to investigate the role of interface configuration in affecting the OSC performance. Referring to Fig. 8, we set and , and we analyze the importance of interfacial length by considering devices with an increasing density of interpenetrating structures, starting from a biplanar device and then taking decreasing values for the rod width . Model parameters are the same as in the previous simulations and the exciton generation rate is .
Fig. 12 illustrates the computed short circuit current density as a function of the interfacial length for the various polaron pair dissociation rate models we previously considered in this section. In all cases, current saturation is predicted for high densities of nanostructures due to the depletion of excitons in the interface area that in turn is a consequence of the abundance of dissociation sites. Computed saturation levels greatly differ among the three choices of the model for , in accordance with the analysis of Sect. 5.2. Fig. 12 also shows that when a biplanar device is considered, using model (C) a higher short circuit current density is obtained compared to the other approaches. An explanation of this result is that the electric field in this case is actually vertically directed and this fact, combined with the assumption that dissociation occurs only in the normal direction, brings to overestimate its rate (cf. the solid lines in Fig. 6). Qualitatively similar results have been obtained in [8, 43, 44].
Also the orientation of the interface is expected to play a role in determining device operation and the following set of simulations aims to investigate this issue. This is a distinctive feature of our model that, to our knowledge, has not been treated in previous works. For a proper analysis, we allow the orientation of the donoracceptor interface to change while its overall length remains almost constant, in order to single out the effect of the former and analyze it.
The considered device geometry is shown in Fig. 13, where the number of rods is kept constant to four for each material but the incidence angle is varied in a range from to . The geometric data are , and .
Since the changes in the amplitude for are small, the interface length does not vary significantly (as demonstrated by Fig. 14(a)) and we expect model (A) to be quite insensitive to such small modifications since mainly depends on the potential drop across the electrodes. Concerning with model (B), we again do not expect a relevant sensitivity to such variations in the interface morphology since the changes in and should balance in the overall contribution. We instead expect model (C) to be most sensitive since the normal field that is screened at the interface may experience significant variations as a function of the angle .
Our expectations are confirmed by the results in Fig. 14(b), showing that the performance of the device in terms of computed short circuit current density does not vary with when models (A) and (B) are considered, while if model (C) is used, an increase of the short circuit current density is observed as soon as the inclination of the nanorod structure is modified with respect to the initial configuration. This behavior can be explained as follows. The choice of model (15) predicts an increase of the dissociation for negative values of the normal electric field that is higher than the reduction for positive values of . Since at short circuit the electric field can be reasonably assumed to be directed along the axis (i.e., from the cathode to the anode), the sides of each rod experience opposite normal fields. As a result, the overall effect is dominated by the contribution of the sides with negative fields and dissociation is enhanced.
5.4. The Case of a Complex Interface Morphology
In this concluding section, we test the versatility of the model proposed in the present article in dealing with a very complex internal morphology as that shown in Fig. 15. In this regard, it is important to notice that the use of the microscale model (2)(5) would require an extremely fine grid resolution to accurately describe the volumetric terms in the active layer around the donoracceptor interface, while the use of the macroscale model (11) has the twofold advantage of considerably simplifying the design of the computational mesh and reducing the size of the nonlinear algebraic system to be solved.
Fig. 16(a) illustrates the computed charge carrier density at short circuit () for the geometry of Fig. 15, where the domain is a square sided, with exciton generation rate and using model (B) for . Notice in particular that the densities assume much higher magnitudes compared to those of Fig. 10. This is a consequence of the complexity of the geometry, where donor and acceptor form deadend areas in which the charges are trapped and experience recombination. In Fig. 16(b) we show again a comparison of the currentvoltage characteristics obtained using the three different polaron dissociation rate models. The differences among the obtained characteristic lines are reduced with respect to the previous simulated cases. In particular the computed short circuit current densities attain closer values with respect to more regular morphologies, such as that of Fig. 8, for comparable values of the interface length (approximately ), see Fig. 12. This is probably to be ascribed to the tortuosity of device internal morphology which makes interface recombination effects more significant than in the case of a more regular internal structure.
6. Concluding Remarks and Future Perspectives
The research activity object of the present article is a continuation of the mathematical study of organic photovoltaic devices started in [17] and is focused to:
 :

the accurate and computationally efficient modeling of photoconversion mechanisms occurring at the interface separating the acceptor and donor layers;
 :

the investigation of the impact of the interface morphology and of polaron pair dissociation on device performance.
With this aim, we propose a twoscale (micro and macroscale) multidimensional model for organic solar cell devices with arbitrary interface geometries. The microscale model is a system of incompletely parabolic nonlinear PDEs in driftdiffusion form set in a heterogeneous domain. The macroscale model is obtained via a microtomacro scale transition that consists of averaging the mass balance equations in the direction normal to the interface, giving rise to nonlinear transmission conditions parametrized by the interfacial width. This averaging procedure is similar to modelreduction techniques used in porous media with thin fractures [31], in reaction problems with moving reaction fronts [29] and in electrochemical transport across biological membranes [35]. The fact that in the macroscale model the interface is reduced to a zero width surface is further exploited to account for the local dependence of the polaron dissociation rate on the electric field orientation, which is the main advantage –together with the computational cost reduction– of our approach, as compared to previous multidimensional models [8, 44, 26].
Extensive numerical simulations of realistic device structures are carried out to study the performance of the proposed models and the impact of the lumping procedure. First, onedimensional transient simulations under different working conditions are carried out to verify the accuracy of the macroscale model with respect to the microscale system. Results indicate that in the physically reasonable range of values for the parameter the relative discrepancy between the micro and macroscale formulations is consistently below 10%. Twodimensional realistic device structures with various interface morphologies are then numerically investigated to assess the impact of our novel model for on the main device properties (short circuit current and open circuit voltage). Simulation results indicate that, if the electric field orientation relative to the interface is taken into due account, the device performance is determined not only by the total interface length but also by its shape.
Research topics currently under scrutiny include:
 :

application of the proposed computational model to the study of more complex threedimensional morphologies, as considered in [27];
 :
 :

extension of the model to the general case where and are free boundaries to be determined for each and at each time level ;
 :
7. Acknowledgements
The authors wish to thank Dr. Dario Natali from Dipartimento di Elettronica e Informazione, Politecnico di Milano and Dr. Mosè Casalegno from Dipartimento di Chimica “Giulio Natta”, Politecnico di Milano for many fruitful and stimulating discussions. The third author was supported by the M.U.R.S.T. grant nr. 200834WK7H005 Adattività Numerica e di Modello per Problemi alle Derivate Parziali (20102012).
Appendix A Finite Element Discretization
In this appendix, for sake of completeness, we provide more detail about the Finite Element (FE) discretization of the linear problem resulting from the application of time semidiscretization and linearization to the equation system (11) as schematically described in Sect.4. Before we proceed we need to introduce some notation.
The time semidiscretization consists of approximating the time derivative of the generic quantity (representing any of the unknowns in (11)) as
(17) 
where is the index of the current time step and is the order of the adopted BDF formula. The notation is used to group together the terms that depend only on results from past time steps and is therefore a known quantity at the th time integration level.
To treat the spatial discretization of the problem, we assume only for ease of presentation that is a rectangular domain, as depicted in Fig. 3(b), but the approach remains completely valid also in the threedimensional case, provided to replace “triangle” by “tethrahedron” and “edge” by “face”. Let denote a conformal partition into open triangles of the computational domain , being the maximum diameter over all triangles, and let and denote the finite element partitions of the subregions and , such that is their separating interface, consisting of the union of a set of edges of .
We introduce the following finite dimensional spaces of FE functions:
(18a)  
(18b)  
(18c)  
(18d)  
(18e) 
Let be such that and . Then, we denote by
(19) 