Solving Diffusion Equations Research Articles

Positive and free energy satisfying schemes for diffusion with interaction potentials

In this paper, we design and analyze second order positive and free energy satisfying schemes for solving diffusion equations with interaction potentials. The semi-discrete scheme is shown to conserve mass, preserve solution positivity, and satisfy a discrete free energy dissipation law for nonuniform meshes. These properties for the fully-discrete scheme (first order in time) remain preserved without a strict restriction on time steps. For the fully second order (in both time and space) scheme, a local scaling limiter is introduced to restore solution positivity when necessary. It is proved that such limiter does not destroy the second order accuracy. In addition, these schemes are easy to implement, and efficient in simulations. Both one and two dimensional numerical examples are presented to demonstrate the performance of these schemes.

Changing Diffusion–Convection Modes in Ternary Mixtures with a Diluent Gas

The characteristic features of the diffusion of helium and nitrous oxide diluted by propane to the same extent are experimentally studied at different pressures and concentrations of the diluent gas in the mixture. It is found that, at a certain pressure, convection currents appear in the system which intensify the partial mass transfer, which uncommon for diffusion. An increase in pressure intensifies the mixing of the components of the mixture. An increase in the concentration of the diluent gas in the mixture leads to the stabilization of convection modes. Solving diffusion equations shows that, in systems with a ballast gas, conditions for the nonlinear distribution of the concentrations of the components of the mixture along the length of the channel are implemented which can lead to the appearance of an extremum of the density of the mixture (i.e., inversion of the density gradient of the mixture), which induces the appearance of convection currents. The stability of the mechanical equilibrium of the ternary mixture in a vertical cylindrical channel under the action of gravity and at the predefined concentration gradients of the components is numerically studied. The boundary of the diffusion–convection kinetic transition is determined. The experimental and numerical results are compared in the coordinates of the Rayleigh number and pressure.

Charge Transfer through Redox Molecular Junctions in Nonequilibrated Solvents.

Molecular conduction operating in dielectric solvent environments is often described using kinetic rates based on the Marcus theory of electron transfer at a molecule-metal electrode interface. However, the successive nature of charge transfer in such a system implies that the solvent does not necessarily reach equilibrium in such processes. Here we generalize the theory to account for solvent nonequilibrium and consider a molecular junction consisting of an electronic donor-acceptor system coupled to two metallic electrodes and placed in a polarizable solvent. We determine the nonequilbrium distribution of the solvent by solving diffusion equations in the strong- and weak-friction limits and calculate the charge current and its fluctuating behavior. In extreme limits, the absence of the solvent or fast solvent relaxation, the charge-transfer statistics is Poissonian, while it becomes correlated by the dynamic solvent between these limits. A Kramers-like turnover of the nonequilibrium current as a function of the solvent damping is found. Finally, we propose a way to tune the solvent-induced damping using geometrical control of the solvent dielectric response in nanostructured solvent channels.

The finite volume scheme preserving maximum principle for diffusion equations with discontinuous coefficient

We construct a new scheme whose primary unknowns include both cell-centered unknowns and edge unknowns on discontinuous line to solve diffusion equations with discontinuous coefficient. First, two linear fluxes are given on both sides of the cell edge, respectively. In order to deal with the defect of the existing scheme preserving maximum principle for solving diffusion problem with discontinuous coefficient, in addition to cell-centered unknowns, we also introduce cell-edge unknowns on discontinuous line as basic unknowns. Second, the conservative flux is constructed by using nonlinear weighted combination of these two linear fluxes. For the cell-edge unknowns on discontinuous line, we add an equation by using the continuity of normal flux. Compared to the classical cell-centered nonlinear finite volume scheme, the introduction of cell-edge unknowns on discontinuous edge is the key point for our scheme to solve diffusion equations with discontinuous coefficient. Then we prove that the scheme satisfies the discrete maximum principle. Based on this, the existence of a solution for the scheme is also obtained. Numerical results are presented to show that our scheme obtains almost second order accuracy for solution on random meshes, preserves discrete maximum principle, and is superior to the existing scheme (in Sheng and Yuan (2011)) preserving the maximum principle on dealing with the problems with strong anisotropic discontinuous coefficient on distorted meshes.

Neumann-Weber integral transform for complex indices

Diffusion and heat equations are commonly investigated in mathematical physics and are solvable for potentials in polar coordinates with a separation into a radial and an angular equation. While the angular equation can be solved easily, a common method for solving the radial part consists in the application of the Neumann-Weber integral transform. The Neumann-Weber integral transform, however, has only been shown to be valid for real indices of Bessel functions. In this work, we generalize the Neumann-Weber transform to complex Bessel indices. The back transform then becomes dependent on zeros of Hankel functions, and we provide useful information for its numerical implementation. The results are relevant for solving diffusion equations and heat equations around cylindrical objects.

Finite-time fluctuation theorem for diffusion-influenced surface reactions on spherical and Janus catalytic particles

A finite-time fluctuation theorem for the diffusion-influenced surface reaction is investigated for spherical and Janus catalytic particles. The finite-time rates and thermodynamic force are analytically calculated by solving diffusion equations with the special boundary conditions of the finite-time fluctuation theorem. Theory is compared with numerical simulations carried out with two different methods: a random walk algorithm and multiparticle collision dynamics.

An agent-based model for drug-radiation interactions in the tumour microenvironment: Hypoxia-activated prodrug SN30000 in multicellular tumour spheroids.

Multicellular tumour spheroids capture many characteristics of human tumour microenvironments, including hypoxia, and represent an experimentally tractable in vitro model for studying interactions between radiotherapy and anticancer drugs. However, interpreting spheroid data is challenging because of limited ability to observe cell fate within spheroids dynamically. To overcome this limitation, we have developed a hybrid continuum/agent-based model (ABM) for HCT116 tumour spheroids, parameterised using experimental models (monolayers and multilayers) in which reaction and diffusion can be measured directly. In the ABM, cell fate is simulated as a function of local oxygen, glucose and drug concentrations, determined by solving diffusion equations and intracellular reactions. The model is lattice-based, with cells occupying discrete locations on a 3D grid embedded within a coarser grid that encompasses the culture medium; separate solvers are employed for each grid. The generated concentration fields account for depletion in the medium and specify concentration-time profiles within the spheroid. Cell growth and survival are determined by intracellular oxygen and glucose concentrations, the latter based on direct measurement of glucose diffusion/reaction (in multilayers) for the first time. The ABM reproduces known features of spheroids including overall growth rate, its oxygen and glucose dependence, peripheral cell proliferation, central hypoxia and necrosis. We extended the ABM to describe in detail the hypoxia-dependent interaction between ionising radiation and a hypoxia-activated prodrug (SN30000), again using experimentally determined parameters; the model accurately simulated clonogenic cell killing in spheroids, while inclusion of reversible cell cycle delay was required to account for the marked spheroid growth delay after combined radiation and SN30000. This ABM of spheroid growth and response exemplifies the utility of integrating computational and experimental tools for investigating radiation/drug interactions, and highlights the critical importance of understanding oxygen, glucose and drug concentration gradients in interpreting activity of therapeutic agents in spheroid models.

Simulation of the ion-induced shock waves effects on the transport of chemically reactive species in ion tracks

The passage of energetic ions through tissue initiates a series of physico-chemical events which leads to biodamage. The study of such ion-induced biodamage using a multiscale approach to the physics of radiation damage with ions has led to the prediction of shock waves being initiated by concentrated energy deposition along the ion track. In this work the radial energy deposition around carbon ion paths, calculated solving diffusion equations, is augmented with the inclusion of more energetic δ-electrons. The radial profiles of energy deposition and the induced concentrations of free radicals are used to simulate the shock waves by means of reactive classical molecular dynamics, which predict a characteristic distribution of reactive chemical species which may have an as yet unrecognised contribution to the nascent biodamage. The paper also suggests some experimental methods by which such a shock wave may be detected and the predictions of these simulations verified.Graphical abstract

Battery Performance Analysis Combined with Circuit Simulation and Electrochemical Calculation

Lithium-ion batteries (LIBs) have become one of the popular energy storage devices in various high-power applications, such as power grid systems and electric vehicles. The advancement of simulation techniques is the important task related to the design, performance, and safety of such systems. Currently, simulations of the battery performance using the multiple resistor-capacitor branches equivalent circuit model (Foster-type model) have been carried out. Although this model is easy to calculate the battery voltage and the current through the cell, it is a behavioral model, deviating from the actual electrode structure and the electrochemical reaction. Therefore, it is unlikely that the equivalent circuit represents the electrochemical physics of the battery performance. Our objective here is to develop a unique equivalent circuit model working in conjunction with numerical calculations according to electrochemical theories. Authors demonstrate that details of the electrochemical reaction in the porous electrode, such as distributions of the overvoltage, the reaction current, and the ion concentration, can be calculated on a circuit simulator by using this model. In a simulation technique presented below, both currents through each electric elements and terminal voltages on the equivalent circuit are calculated by the circuit simulator running on MATLAB/Simulink, and these are used as boundary conditions to calculate the Butler-Volmer equation combined with diffusion equations. Thereby a circuit simulation and electrochemical calculations are coupled and the battery performance is calculated according to both theories of the electric circuit and the electrochemistry. The Butler-Volmer equation and diffusion equations are solved by using numerical calculation function implemented in the Simulink. In developing LIBs model, it is very important to express the electrochemical reaction of the porous electrode. This is because the reactive sites spread in the electrode thickness direction, so reaction conditions vary depending on the position of the electrode depth. Therefore, in a equivalent circuit model of LIBs, the porous electrode structure must be reflected. In this simulation technique, a equivalent circuit of a transmission line model (TLM) of the porous electrode is adopted (see Figure 1). The electrolyte resistance of pore inside (Re’), the reaction resistance (Rf’), and the electric double layer capacitance (Cdl’) are given for each segment of electrode pores divided into 100. By using a TLM, it is expected that distributions of the potential and the current in the actual porous electrode can be expressed. Figure 2 shows one segment in a TLM and the method of estimating Rf ' is described below. First, the electric double layer voltage (= overvoltage) is calculated by a circuit simulation, and then both the reaction current and the Li-flux of the intercalation/deintercalation through the electrolyte-electrode interface are calculated by using the Butler-Volmer equation as a function of the overvoltage. Second, the Li ion concentration [Li+] in the electrolyte, the Li concentration [Li] and the vacant site concentration [θ] in the electrode are estimated by solving diffusion equations under the boundary condition of the Li-flux. Moreover, [Li+], [Li], and [θ] are fed back to the Butler-Volmer equation. Third, value of the Rf’ is determined by dividing the overvoltage by the reaction current. Thus, Rf’ is modeled as a variable resistor which changes with overvoltage. Forth, the calculated Rf’ is transferred to the equivalent circuit. Thereby the equivalent circuit is updated from the initial state and the battery current, the electric double layer charging current, the reaction current, and the overvoltage are recalculated on a circuit simulator. The above is the calculation sequence in one calculation step, and this calculation sequence is repeated for each calculation step. The above calculation is carried out for the variable Rf’ element of all segments in Figure 1. Reaction conditions (ex. overvoltage) for each segment is automatically calculated according to a TLM on a circuit simulator. Therefore, Rf’ is a variable resistance which depends on the electrode thickness direction. Electrolyte resistances of the pore outside (Rsol’), Re’, and Cdl’ are estimated from the electrochemical impedance measurement. Diffusion coefficients etc. are determined from literature values. Figure 3 shows the transient current curve at the potential step of 0.5 V obtained with the positive electrode of LIBs. The experimental and the simulated currents were found to be in excellent agreement. In this simulation technique, the behavior of the electric double layer charging current immediately after the potential step and the steady reaction current can be calculated. Moreover, distributions of reaction current and overvoltage in the electrode thickness direction can be estimated (see Figure 4). Figure 1

A Lattice Boltzmann Scheme for Diffusion Equation in Spherical Coordinate

Lattice Boltzmann models for diffusion equation are generally in Cartesian coordinate system. Very few researchers have attempted to solve diffusion equation in spherical coordinate system. In the lattice Boltzmann based diffusion model in spherical coordinate system extra term, which is due to variation of surface area along radial direction, is modeled as source term. In this study diffusion equation in spherical coordinate system is first converted to diffusion equation which is similar to that in Cartesian coordinate system by using proper variable. The diffusion equation is then solved using standard lattice Boltzmann method. The results obtained for the new variable are again converted to the actual variable. The numerical scheme is verified by comparing the results of the simulation study with analytical solution. A good agreement between the two results is established.

Phantom Preparation and Optical Property Determination

Tissue-like optical phantoms are important in testing new imaging algorithms. Homogeneous optical phantoms with determined optical properties are the first step of making a proper heterogeneous phantom for multi-modality imaging. Typical recipes for such phantoms consist of epoxy resin, hardener, India ink and titanium oxide. By altering the concentration of India ink and titanium oxide, we are able to get multiple homogeneous phantoms with different absorption and scattering coefficients by carefully mixing all the ingredients. After fabricating the phantoms, we need to find their individual optical properties including the absorption and scattering coefficients. This is achieved by solving diffusion equation of each phantom as a homogeneous slab under canonical illumination. We solve the diffusion equation of homogeneous slab in frequency domain and get the formula for theoretical measurements. Under our steady-state diffused optical tomography (DOT) imaging system, we are able to obtain the real distribution of the incident light produced by a laser. With this source distribution we got and the formula we derived, numerical experiments show how measurements change while varying the value of absorption and scattering coefficients. Then we notice that the measurements alone will not be enough for us to get unique optical properties for steady-state DOT problem. Thus in order to determine the optical properties of a homogeneous slab we want to fix one of the coefficients first and use optimization methods to find another one. Then by assemble multiple homogeneous slab phantoms with different optical properties, we are able to obtain a heterogeneous phantom suitable for testing multi-modality imaging algorithms. In this paper, we describe how to make phantoms, derive a formula to solve the diffusion equation, demonstrate the non-uniqueness of steady-state DOT problem by analysing some numerical results of our formula, and finally propose a possible way to determine optical properties for homogeneous slab for our future work.

Searching for an optimal shape parameter for solving a partial differential equation with the radial basis functions method

This article presents a procedure for searching for an optimal shape parameter for the solution of partial differential equations with the corresponding initial and boundary conditions, where the solution of the problem is unknown. In recent years, radial basis function methods have emerged as alternative computing methods in the scientific computing community.The numerical solution of partial differential equations has usually been obtained by using finite difference methods, finite element methods (FEMs), boundary elements methods or finite volume methods. In our case, we use the multiquadric radial basis function, Gershgorin’s theorem and the Newton method for searching an optimal shape parameter for solving diffusion equations. More cases are presented, the results of which are compared with the results obtained by the FEM.

Greedy algorithms for diffuse optical tomography reconstruction

Diffuse optical tomography (DOT) is a noninvasive imaging modality that reconstructs the optical parameters of a highly scattering medium. However, the inverse problem of DOT is ill-posed and highly nonlinear due to the zig-zag propagation of photons that diffuses through the cross section of tissue. The conventional DOT imaging methods iteratively compute the solution of forward diffusion equation solver which makes the problem computationally expensive. Also, these methods fail when the geometry is complex. Recently, the theory of compressive sensing (CS) has received considerable attention because of its efficient use in biomedical imaging applications. The objective of this paper is to solve a given DOT inverse problem by using compressive sensing framework and various Greedy algorithms such as orthogonal matching pursuit (OMP), compressive sampling matching pursuit (CoSaMP), and stagewise orthogonal matching pursuit (StOMP), regularized orthogonal matching pursuit (ROMP) and simultaneous orthogonal matching pursuit (S-OMP) have been studied to reconstruct the change in the absorption parameter i.e, Δα from the boundary data. Also, the Greedy algorithms have been validated experimentally on a paraffin wax rectangular phantom through a well designed experimental set up. We also have studied the conventional DOT methods like least square method and truncated singular value decomposition (TSVD) for comparison. One of the main features of this work is the usage of less number of source–detector pairs, which can facilitate the use of DOT in routine applications of screening. The performance metrics such as mean square error (MSE), normalized mean square error (NMSE), structural similarity index (SSIM), and peak signal to noise ratio (PSNR) have been used to evaluate the performance of the algorithms mentioned in this paper. Extensive simulation results confirm that CS based DOT reconstruction outperforms the conventional DOT imaging methods in terms of computational efficiency. The main advantage of this study is that the forward diffusion equation solver need not be repeatedly solved.

A multidimensional unified gas-kinetic scheme for radiative transfer equations on unstructured mesh

In order to extend the unified gas kinetic scheme (UGKS) to solve radiative transfer equations in a complex geometry, a multidimensional asymptotic preserving implicit method on unstructured mesh is constructed in this paper. With an implicit formulation, the CFL condition for the determination of the time step in UGKS can be much relaxed, and a large time step is used in simulations. Differently from previous direction-by-direction UGKS on orthogonal structured mesh, on unstructured mesh the interface flux transport takes into account multi-dimensional effect, where gradients of radiation intensity and material temperature in both normal and tangential directions of a cell interface are included in the flux evaluation. The multiple scale nature makes the UGKS be able to capture the solutions in both optically thin and thick regions seamlessly. In the optically thick region the condition of cell size being less than photon's mean free path is fully removed, and the UGKS recovers a solver for diffusion equation in such a limit on unstructured mesh. For a distorted quadrilateral mesh, the UGKS goes to a nine-point scheme for the diffusion equation, and it naturally reduces to the standard five-point scheme for a orthogonal quadrilateral mesh. Numerical computations covering a wide range of transport regimes on unstructured and distorted quadrilateral meshes will be presented to validate the current approach.

Fourth-order numerical solutions of diffusion equation by using SOR method with Crank-Nicolson approach

The aim of this paper is to investigate the effectiveness of the Successive Over-Relaxation (SOR) iterative method by using the fourth-order Crank-Nicolson (CN) discretization scheme to derive a five-point Crank-Nicolson approximation equation in order to solve diffusion equation. From this approximation equation, clearly, it can be shown that corresponding system of five-point approximation equations can be generated and then solved iteratively. In order to access the performance results of the proposed iterative method with the fourth-order CN scheme, another point iterative method which is Gauss-Seidel (GS), also presented as a reference method. Finally the numerical results obtained from the use of the fourth-order CN discretization scheme, it can be pointed out that the SOR iterative method is superior in terms of number of iterations, execution time, and maximum absolute error.

Three-dimensional numerical simulation of solid oxide fuel cell cathode based on lattice Boltzmann method with sub-grid scale models

In the present study, three-dimensional numerical simulations are carried out to predict the electrochemical characteristics of solid oxide fuel cell (SOFC) cathodes based on a coupled method using sub-grid scale (SGS) model in a lattice Boltzmann method (LBM). Lattice Boltzmann method is used to solve the governing equations. In each coarse computational grid, local gas diffusivities, ion and electron conductivities are obtained by directly solving diffusion equations using the original fine sub-grid scale information. Two types of SGS models are introduced in this study, i.e. isotropic and anisotropic local diffusivities and conductivities. Proposed SGS local diffusivities and conductivities can maintain detailed information of the original fine microstructure even when coarser grid is used in the calculation. In the anisotropic SGS model, weighted summation of particle distribution function is applied in the LBM to maintain the invariances of local concentrations and potentials. From the tortuosity factor and overpotential calculation results, it is concluded that the proposed SGS models can drastically improve the computational accuracy without costing additional calculation time. Moreover, SGS model considering geometric anisotropies inside the calculation grid presents more precise results than the isotropic SGS model.

Abstract 4523: Agent-based models for prediction of tumor spheroid response to radiation and the hypoxia-activated prodrug SN30000

Abstract Introduction: Hypoxia-activated prodrugs (HAPs) are a promising strategy to specifically target radioresistant hypoxic cells in tumors. We have used multicellular layers (MCLs) and steady state Green’s function pharmacokinetic/pharmacodynamic (PK/PD) modelling in microvascular networks to select SN30000 as a tirapazamine analogue with optimal spatial complementarity with radiation. However these models have several limitations such as lack of time dependent PK. Multicellular tumor spheroids are a more tractable model than MCLs for investigating PD parameters. However, the challenges of interpreting spheroid experiments in terms of spatial complementarity limits their application. We have developed an in silico agent-based (AB) tumor spheroid model, in which cell fate can be tracked as a function of time and is determined by spatially-varying concentrations of oxygen, glucose and therapeutic agents (radiation and HAP). The hypothesis is that our AB model can predict tumor spheroid response to radiation and SN30000, both alone and in combination. Methods: In the AB model, spheroid growth is simulated as a function of oxygen, glucose and drug concentrations determined by solving diffusion equations using oxygen consumption measured in HCT116 monolayers and glucose and SN30000 diffusion and metabolism measured in HCT116 MCLs and stirred cell suspensions, linked with intracellular reaction equations. Cell fate is determined by oxygen, glucose and drug concentrations. Spheroid diameter, cell survival and hypoxic fraction of HCT116 spheroids grown at a range of oxygen and initial glucose concentrations were compared to model predictions, as was glucose depletion from the culture medium. Radiation (cobalt-60) and SN30000 induced clonogenic cell killing under oxia and anoxia was determined using HCT116 monolayers. Model simulations of HCT116 spheroid response to radiation and SN30000 were compared to experimental results. Results: The model reproduced known features of spheroids including oxygen and glucose gradients, rapid cell proliferation at the periphery and central hypoxia and necrosis. Measured HCT116 spheroid diameters increased linearly with time, at rates that decreased at lower ambient O2 concentrations. Spheroid growth and glucose depletion was well-fitted by the AB model. Good agreement was found between simulated and measured clonogenic cell killing by radiation and SN30000, while the model underestimated spheroid growth delay by radiation and SN30000 than experimental observation. Discussion: An oxygen and glucose dependent AB model for growth and PD response (cell killing and growth delay) of HCT116 spheroids has been developed and calibrated using experimentally determined parameters. Further development of this model will link cell fates (including cell cycle arrest) explicitly to DNA damage responses and O2 and glucose-dependence of cell proliferation and survival. Citation Format: Xinjian Mao, Gib Bogle, Sarah McManaway, William Wilson, Kevin Hicks. Agent-based models for prediction of tumor spheroid response to radiation and the hypoxia-activated prodrug SN30000 [abstract]. In: Proceedings of the American Association for Cancer Research Annual Meeting 2017; 2017 Apr 1-5; Washington, DC. Philadelphia (PA): AACR; Cancer Res 2017;77(13 Suppl):Abstract nr 4523. doi:10.1158/1538-7445.AM2017-4523

Using Newton basis functions for solving diffusion equation with nonlocal boundary condition

In this paper we introduce a meshless method based on spatial trial spaces spanned by the Newton basis functions. The Newton basis arising from a pivoted Cholesky factorization turns out to be stable and computationally cheap while being orthonormal in the “native” Hilbert space of the kernel. This method can overcome on some of challenges of the full meshless kernel-based methods such as ill-conditioning associated with the dense linear systems that arise. This approach is used to obtain the numerical solution of time-dependent partial differential equations (PDEs), diffusion equation with two types of nonlocal boundary condition. The time-dependent PDE is discretized in space by the collocation method based on Newton basis functions as trial functions and it is converted to a system of ordinary differential equations in time. Stability of the scheme is discussed. Numerical examples are provided to show the accuracy and easy implementation of the novel proposed method.

A comparison study of solving diffusion equations with different algorithm methods

A comparison study for solving diffusion equations with different algorithm methods is studied to understand the oxygen vacancy defect transport under the electric field. We compare computational efficiency and numerical accuracy with different algorithm methods, including finite difference, finite element (COMSOL), and Fourier-Chebysev spectral methods. All the results of oxygen vacancy distribution under an electric field from different algorithm methods are compared with the analytical solution results. Two kinds of boundary conditions are used in solving diffusion equations and the absolute error of different methods are discussed. The main purpose of these results is to provide guidance for studying the role of point defect transport in the degradation and breakdown of devices.

Discrete Conservation Laws on Evolving Surfaces

In this paper we shall introduce a simple effective numerical method for solving diffusion equations on evolving surfaces. Instead of using the finite element method, we shall follow the traditional way to approximate the regular evolving surfaces under consideration by evolving triangular meshes. The key idea of our two-step algorithm is to develop an intrinsic method to compute directly the partial derivatives of functions and indirectly the mean curvatures of the evolving triangular meshes. Most importantly, conservational laws on evolving triangular meshes for diffusion equations are fulfilled.