Riccardo Sacco

Carlo de Falco, Riccardo Sacco, Maurizio Verri

This article is an attempt to provide a self consistent picture, including existence analysis and numerical solution algorithms, of the mathematical problems arising from modeling photocurrent transients in Organic-polymer Solar Cells (OSCs). The mathematical model for OSCs consists of a system of nonlinear diffusion-reaction partial differential equations (PDEs) with electrostatic convection, coupled to a kinetic ordinary differential equation (ODE). We propose a suitable reformulation of the model that allows us to prove the existence of a solution in both stationary and transient conditions and to better highlight the role of exciton dynamics in determining the device turn-on time. For the numerical treatment of the problem, we carry out a temporal semi-discretization using an implicit adaptive method, and the resulting sequence of differential subproblems is linearized using the Newton-Raphson method with inexact evaluation of the Jacobian. Then, we use exponentially fitted finite elements for the spatial discretization, and we carry out a thorough validation of the computational model by extensively investigating the impact of the model parameters on photocurrent transient times.

Carlo de Falco, Matteo Porro, Riccardo Sacco et al.

In this article, we continue our mathematical study of organic solar cells (OSCs) and propose a two-scale (micro- and macro-scale) 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 micro-to-macro 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.

Chiara Lelli, Riccardo Sacco, Paola Causin et al.

An adequate control of cell response in tissue engineering applications is of utmost importance to obtain products suitable to clinical practice. This paper is the first part of a series of two connected publications in which we study via mathematical tools the cultivation in bioreactors of articular chondrocytes. The proposed model combines poroelastic theory of mixtures and cellular population models into a framework including stress state and oxygen tension as main determinants of engineered culture evolution. The special mechanosensitivity of articular chondrocytes to the surrounding environment is accounted for in the model through the novel concept of "force isotropy" acting on the cell which is assumed as the promoting factor of the production of new cells or extracellular matrix.

Emanuela Abbate, Matteo Porro, Thierry Nieus et al.

In this article we propose and investigate a hierarchy of mathematical models based on partial differential equations (PDE) and ordinary differential equations (ODE) for the simulation of the biophysical phenomena occurring in the electrolyte fluid that connects a biological component (a single cell or a system of cells) and a solid-state device (a single silicon transistor or an array of transistors). The three members of the hierarchy, ordered by decreasing complexity, are: (i) a 3D Poisson-Nernst-Planck (PNP) PDE system for ion concentrations and electric potential; (ii) a 2D reduced PNP system for the same dependent variables as in (i); (iii) a 2D area-contact PDE system for electric potential coupled with a system of ODEs for ion concentrations. The backward Euler method is adopted for temporal semi-discretization and a fixed-point iteration based on Gummel's map is used to decouple system equations. Spatial discretization is performed using piecewise linear triangular finite elements stabilized via edge-based exponential fitting. Extensively conducted simulation results are in excellent agreement with existing analytical solutions of the PNP problem in radial coordinates and experimental and simulated data using simplified lumped parameter models.

Chiara Lelli, Riccardo Sacco, Paola Causin et al.

In Part I of this article we have developed a novel mechanobiological model of a Tissue Engineering process that accounts for the mechanisms through which an isotropic or anisotropic adherence condition regulates the active functions of the cells in the construct. The model expresses mass balance and force equilibrium balance for a multi-phase mixture in a 3D computational domain and in time dependent conditions. In the present Part II, we study the mechanobiological model in a simplified 1D geometrical setting with the purpose of highlighting the ability of the formulation to represent the influence of force isotropy and nutrient availability on the growth of the tissue construct. In particular, an example of isotropy estimator is proposed and coded within a fixed-point solution map that is used at each discrete time level for system linearization and subsequent finite element approximation of the linearized equations. Extensively conducted simulations show that: 1) the spatial and temporal evolution of the cellular populations are in good agrement with the local growth/production conditions predicted by the mechanobiological stress-dependent model; and 2) the isotropy indicator and all model variables are strongly influenced by both maximum cell specific growth rate and mechanical boundary conditions enforced at the interface between the biomass construct and the interstitial fluid.

Riccardo Sacco, Aurelio Giancarlo Mauri, Giovanna Guidoboni

This work focuses on a class of elliptic boundary value problems with diffusive, advective and reactive terms, motivated by the study of three-dimensional heterogeneous physical systems composed of two or more media separated by a selective interface. We propose a novel approach for the numerical approximation of such heterogeneous systems combining, for the first time: (1) a dual mixed hybrid (DMH) finite element method (FEM) based on the lowest order Raviart-Thomas space (RT0); (2) a Three-Field (3F) formulation; and (3) a Streamline Upwind/Petrov-Galerkin (SUPG) stabilization method. Using the abstract theory for generalized saddle-point problems and their approximation, we show that the weak formulation of the proposed method and its numerical counterpart are both uniquely solvable and that the resulting finite element scheme enjoys optimal convergence properties with respect to the discretization parameter. In addition, an efficient implementation of the proposed formulation is presented. The implementation is based on a systematic use of static condensation which reduces the method to a nonconforming finite element approach on a grid made by three-dimensional simplices. Extensive computational tests demonstrate the theoretical conclusions and indicate that the proposed DMH-RT0 FEM scheme is accurate and stable even in the presence of marked interface jump discontinuities in the solution and its associated normal flux. Results also show that in the case of strongly dominating advective terms, the proposed method with the SUPG stabilization is capable of resolving accurately steep boundary and/or interior layers without introducing spurious unphysical oscillations or excessive smearing of the solution front.

Riccardo Sacco, Paolo Airoldi, Aurelio G. Mauri et al.

In this article we address the three-dimensional modeling and simulation of biological ion channels using a continuum-based approach. Our multi-physics formulation self-consistently combines, to the best of our knowledge for the first time, ion electrodiffusion, channel fluid motion, thermal self-heating and mechanical deformation. The resulting system of nonlinearly coupled partial differential equations in conservation form is discretized using the Galerkin Finite Element Method. The validation of the proposed computational model is carried out with the simulation of a cylindrical voltage operated ion nanochannel with K+ and Na+ ions. We first investigate the coupling between electrochemical and fluid-dynamical effects. Then, we enrich the modeling picture by investigating the influence of a thermal gradient. Finally, we add a mechanical stress responsible for channel deformation and investigate its effect on the functional response of the channel. Results show that fluid and thermal fields have no influence in absence of mechanical deformation whereas ion distributions and channel functional response are significantly modified if mechanical stress is included in the model. These predictions agree with biophysical conjectures on the importance of protein conformation in the modulation of channel electrochemical properties.

Riccardo Sacco, Fabio Manganini, Joseph W. Jerome

In this article we address the study of ion charge transport in the biological channels separating the intra and extracellular regions of a cell. The focus of the investigation is devoted to including thermal driving forces in the well-known velocity-extended Poisson-Nernst-Planck (vPNP) electrodiffusion model. Two extensions of the vPNP system are proposed: the velocity-extended Thermo-Hydrodynamic model (vTHD) and the velocity-extended Electro-Thermal model (vET). Both formulations are based on the principles of conservation of mass, momentum and energy, and collapse into the vPNP model under thermodynamical equilibrium conditions. Upon introducing a suitable one-dimensional geometrical representation of the channel, we discuss appropriate boundary conditions that depend only on effectively accessible measurable quantities. Then, we describe the novel models, the solution map used to iteratively solve them, and the mixed-hybrid flux-conservative stabilized finite element scheme used to discretize the linearized equations. Finally, we successfully apply our computational algorithms to the simulation of two different realistic biological channels: 1) the Gramicidin-A channel considered in~\cite{JeromeBPJ}; and 2) the bipolar nanofluidic diode considered in~\cite{Siwy7}.

Aurelio Mauri, Andrea Bortolossi, Giovanni Novielli et al.

In this article we propose two novel 3D finite element models, denoted method A and B, for electron and hole Drift-Diffusion (DD) current densities. Method A is based on a primal-mixed formulation of the DD model as a function of the quasi-Fermi potential gradient, while method B is a modification of the standard DD formula based on the introduction of an artificial diffusion matrix. Both methods are genuine 3D extensions of the classic 1D Scharfetter-Gummel difference formula. The proposed methods are compared in the 3D simulation of a p-n junction diode and of a p-MOS transistor in the on-state regime. Results show that method A provides the best performance in terms of physical accuracy and numerical stability. Method A is then used in the 3D simulation of a n-MOS transistor in the off-state regime including the impact ionization generation mechanism. Results demonstrate that the model is able to accurately compute the I-V characteristic of the device until drain-to-bulk junction breakdown.