Sepic Converter Simulation Matlab

SpinDoctor: a MATLAB toolbox for diffusion MRI simulation Jing-Rebecca Lia,∗, Van-Dang Nguyenb, Try Nguyen Trana, Jan Valdmanc, Cong-Bang Tranga, Khieu Van Nguyena, Duc Thach Son Vua, Hoang An Trana, Hoang Trong An Trana, Thi Minh

Phuong Nguyena

aINRIA Saclay, Equipe DEFI, CMAP, Ecole Polytechnique, Route de Saclay, 91128 Palaiseau Cedex, France Information Theory and Automation of the ASCR, Prague, Czech Republic

Abstract

The complex transverse water proton magnetization subject to diffusion-encoding magnetic field gradient pulses in a heterogeneous medium can be modeled by the multiple compartment Bloch- Torrey partial differential equation. Under the assumption of negligible water exchange between compartments, the time-dependent apparent diffusion coefficient can be directly computed from the solution of a diffusion equation subject to a time-dependent Neumann boundary condition.

This paper describes a publicly available MATLAB toolbox called SpinDoctor that can be used 1) to solve the Bloch-Torrey partial differential equation in order to simulate the diffusion magnetic resonance imaging signal; 2) to solve a diffusion partial differential equation to obtain directly the apparent diffusion coefficient; 3) to compare the simulated apparent diffusion coefficient with a short-time approximation formula.

The partial differential equations are solved by P1 finite elements combined with built-in MATLAB routines for solving ordinary differential equations. The finite element mesh generation is performed using an external package called Tetgen.

SpinDoctor provides built-in options of including 1) spherical cells with a nucleus; 2) cylindrical cells with a myelin layer; 3) an extra-cellular space enclosed either a) in a box or b) in a tight wrapping around the cells; 4) deformation of canonical cells by bending and twisting; 5) permeable membranes; Built-in diffusion-encoding pulse sequences include the Pulsed Gradient Spin Echo and the Oscillating Gradient Spin Echo.

We describe in detail how to use the SpinDoctor toolbox. We validate SpinDoctor simulations using reference signals computed by the Matrix Formalism method. We compare the accuracy and computational time of SpinDoctor simulations with Monte-Carlo simulations and show significant speed-up of SpinDoctor over Monte-Carlo simulations in complex geometries. We also illustrate several extensions of SpinDoctor functionalities, including the incorporation of T2 relaxation, the simulation of non-standard diffusion-encoding sequences, as well as the use of externally generated

Arxiv:1902.01025V2 [Math.Na] 16 Sep 2019

geometrical meshes.

Keywords:

Bloch-Torrey equation, diffusion magnetic resonance imaging, finite elements, simulation, apparent diffusion coefficient.

1. Introduction

Diffusion magnetic resonance imaging is an imaging modality that can be used to probe the tissue micro-structure by encoding the incohorent motion of water molecules with magnetic field gradient pulses. This motion during the diffusion-encoding time causes a signal attenuation from which the apparent diffusion coefficient , (and possibly higher order diffusion terms, can be calculated .

For unrestricted diffusion, the root of the mean squared displacement of molecules is given by ¯x = √2 dim σ0t, where dim is the spatial dimension, σ0 is the intrinsic diffusion coefficient, and t is the diffusion time. In biological tissue, the diffusion is usually hindered or restricted (for example, by cell membranes) and the mean square displacement is smaller than in the case of unrestricted diffusion. This deviation from unrestricted diffusion can be used to infer information about the tissue micro-structure. The experimental parameters that can be varied include 1. the diffusion time (one can choose the parameters of the diffusion-encoding sequence, such as Pulsed Gradient Spin Echo and Oscillating Gradient ).

2. the magnitude of the diffusion-encoding gradient (when the magnetic resonance imaging sig- nal is acquired at low gradient magnitudes, the signal contains only information about the apparent diffusion coefficient, at higher values, Kurtosis imaging becomes possible); 3. the direction of the diffusion-encoding gradient (many directions may be probed, as in high angular resolution diffusion imaging ).

Using diffusion magnetic resonance imaging to get tissue structural information in the mamalian brain has been the focus of much experimental and modeling work in recent years . The pre- dominant approach up to now has been adding the diffusion magnetic resonance imaging signal from simple geometrical components and extracting model parameters of interest. Numerous biophysical models subdivide the tissue into compartments described by spheres, ellipsoids, cylinders, and the extra-cellular space [7–9, 11, 12, 15–19]. Some model parameters of interest include axon diameter and orientation, neurite density, dendrite structure, the volume fraction and size distribution of cylinder and sphere components and the effective diffusion coefficient or tensor of the extra-cellular space.

Numerical simulations can help deepen the understanding of the relationship between the cellular structure and the diffusion magnetic resonance imaging signal and lead to the formulation of appro- priate models. They can be also used to investigate the effect of different pulse sequences and tissue features on the measured signal which can be used for the development, testing, and optimization of novel diffusion magnetic resonance imaging pulse sequences .

Two main groups of approaches to the numerical simulation of diffusion magnetic resonance imaging are 1) using random walkers to mimic the diffusion process in a geometrical configuration; 2) solving

2

the Bloch-Torrey PDE , which describes the evolution of the complex transverse water proton magnetization under the influence of diffusion-encoding magnetic field gradients pulses. The first group is referred to as Monte-Carlo simulations in the literature and previous works include [13, 24–27]. A GPU-based acceleration of Monte-Carlo simulation was proposed in [28, 29]. Some

Software Packages Using This Approach Include

1. Camino Diffusion MRI Toolkit developed at UCL (http://cmic.cs.ucl.ac.uk/camino/); 2. DIFSIM developed at UC San Diego (http://csci.ucsd.edu/projects/simulation.html); 3. Diffusion Microscopist Simulator developed at Neurospin, CEA.

The second group relies on solving the Bloch-Torrey PDE in a geometrical configuration. In [30– 32] a simplifying assumption called the narrow pulse approximation was used, where the pulse duration was assumed to be much smaller than the delay between pulses. This assumption allows the solution of the diffusion equation instead of the more complicated Bloch-Torrey PDE. More generally, numerical methods to solve the Bloch-Torrey PDE. with arbitrary temporal profiles have been proposed in . The computational domain is discretized either by a Cartesian grid [33, 34, 37] or finite elements [30–32, 35, 36]. The unstructured mesh of a finite element discretization appeared to be better than a Cartesian grid in both geometry description and signal approximation . For time discretization, both explicit and implicit methods have been used. In a second order implicit time-stepping method called the generalized α−method was used to allow for high frequency energy dissipation. An adaptive explicit Runge-Kutta Chebyshev method of second order was used in [34, 35]. It has been theoretically proven that the Runge-Kutta Chebyshev method allows for a much larger time-step compared to the standard explicit Euler method . There is an example showing that the Runge-Kutta Chebyshev method is faster than the implicit Euler method in . The Crank-Nicolson method was used in to also allow for second order convergence in time. The efficiency of diffusion magnetic resonance imaging simulations is also improved by either a high-performance FEM computing framework [39, 40] for large-scale simulations on supercomputers or a discretization on manifolds for thin-layer and thin-tube media .

In this paper, we present a MATLAB Toolbox called SpinDoctor that is a simulation pipeline going from the definition of a geometrical configuration through the numerical solution of the Bloch- Torrey PDE to the fitting of the apparent diffusion coefficient from the simulated signal. It also includes two other modules for calculating the apparent diffusion coefficient. The first module is a homogenized apparent diffusion coefficient mathematical model, which was obtained recently using homogenization techniques on the Bloch-Torrey PDE. In the homogenized model, the apparent dif- fusion coefficient of a geometrical configuration can be computed after solving a diffusion equation subject to a time-dependent Neumann boundary condition, under the assumption of negligible wa- ter exchange between compartments. The second module computes the short time approximation formula for the apparent diffusion coefficient. The short time approximation implemented in Spin- Doctor includes a recent generalization of this formula to account for finite pulse duration in the pulsed gradient spin echo. Both of these two apparent diffusion coefficient calculations are sensitive to the diffusion-encoding gradient direction, unlike many previous works where the anisotropy is neglected in analytical model development.

In Summary, Spindoctor

1. solves the Bloch-Torrey PDE in three dimensions to obtain the diffusion magnetic resonance

Imaging Signal;

2. robustly fits the diffusion magnetic resonance imaging signal to obtain the apparent diffusion

Coefficient;

3. solves the homogenized apparent diffusion coefficient model in three dimensions to obtain the

Apparent Diffusion Coefficient;

4. computes the short-time approximation of the apparent diffusion coefficient; 5. computes useful geometrical quantities such as the compartment volumes and surface areas; 6. allows permeable membranes for the Bloch-Torrey PDE (the homogenized apparent diffusion coefficient assumes negligible permeabilty).

7. displays the gradient-direction dependent apparent diffusion coefficient; in three dimensions

Using Spherical Harmonics Interpolation;

SpinDoctor provides the following built-in functionalities: 1. placement of non-overlapping spherical cells (with an optional nucleus) of different radii close

To Each Other;

2. placement of non-overlapping cylindrical cells (with an optional myelin layer) of different radii close to each other in a canonical configuration where they are parallel to the z-axis; 3. inclusion of an extra-cellular space that is enclosed either

(B) In A Rectangular Box;

4. deformation of the canonical configuration by bending and twisting; Built-in diffusion-encoding pulse sequences include

1. The Pulsed Gradient Spin Echo ;

2. the Oscillating Gradient Spin Echo (cos- and sin- type gradients).

Spindoctor Uses The Following Methods:

1. it generates a good quality surface triangulation of the user specified geometrical configuration by calling built-in MATLAB computational geometry functions; 2. it creates a good quality tetrehedra finite elements mesh from the above surface triangulation by calling Tetgen , an external package (executable files are included in the Toolbox

4

3. it constructs finite element matrices for linear finite elements on tetrahedra (P1) using routines

From ;

4. it adds additional degrees of freedom on the compartment interfaces to allow permeability conditions for the Bloch-Torrey PDE using the formalism in ; 5. it solves the semi-discretized FEM equations by calling built-in MATLAB routines for solving ordinary differential equations .

The SpinDoctor toolbox has been developed in the MATLAB R2017b and requires no additional MATLAB toolboxes. The toolbox is publicly available at:

2. Theory

Suppose the user would like to simulate a geometrical configuration of cells with an optional myelin layer or a nucleus. If spins will be leaving the cells or if the user wants to simulate the extra-cellular space (ECS), then the ECS will enclose the geometrical shapes. Let Ωe be the ECS, Ωin

I

the cytoplasm (or the myelin layer) of the ith cell. We denote the interface

And Ωe By Σi, Finally The Outside

boundary of the ECS by Ψ.

2.1. Bloch-Torrey Pde

In diffusion MRI, a time-varying magnetic field gradient is applied to the tissue to encode water diffusion. Denoting the effective time profile of the diffusion-encoding magnetic field gradient by f(t), and letting the vector g contain the amplitude and direction information of the magnetic field gradient, the complex transverse water proton magnetization in the rotating frame satisfies the

(3)

where γ = 2.67513×108 rad s−1T−1 is the gyromagnetic ratio of the water proton, I is the imaginary unit, σl is the intrinsic diffusion coefficient in the compartment Ωl

I. The Magnetization Is A Function

of position x and time t, and depends on the diffusion gradient vector g and the time profile f(t). We denote the restriction of the magnetization in Ωin

I

and M e. Some commonly used time profiles (diffusion-encoding sequences) are: 1. The pulsed-gradient spin echo (PGSE) sequence, with two rectangular pulses of duration δ, separated by a time interval ∆−δ, for which the profile f(t) is

(4)

where t1 is the starting time of the first gradient pulse with t1 + ∆> TE/2, TE is the echo time at which the signal is measured. 2. The oscillating gradient spin echo (OGSE) sequence [4, 45] was introduced to reach short diffusion times. An OGSE sequence usually consists of two oscillating pulses of duration T, each containing n periods, hence the frequency is ω = n 2π T , separated by a time interval τ −T.

(5)

where τ = TE/2. The BTPDE needs to be supplemented by interface conditions. We recall the interface between

And Ωe Is Σi, And The Outside Boundary Of The Ecs

is Ψ. The two interface conditions on Γi are the flux continuity and a condition that incorporates

6

where n is the unit outward pointing normal vector. Similarly, between Ωout

,

x ∈Σi.

0 = Σe∇M E(X, T) · Ne,

x ∈Ψ.

(X, 0) = Ρout,

M e(x, 0) = ρe. where ρ is the initial spin density. The dMRI signal is measured at echo time t = TE > ∆+δ for PGSE and TE > 2σ for OGSE. This

, Ωe}

M(x, TE) dx.

(6)

In a dMRI experiment, the pulse sequence (time profile f(t)) is usually fixed, while g is varied in amplitude (and possibly also in direction). S is usually plotted against a quantity called the b-value. The b-value depends on g and f(t) and is defined as

2

.

For Pgse, The B-Value Is :

b(g, δ, ∆) = γ2∥g∥2δ2 (∆−δ/3) .

(7)

For the cosine OGSE with integer number of periods n in each of the two durations σ, the corre-

4N2Π2 = Γ2∥G∥2 Σ

ω2 .

(8)

The reason for these definitions is that in a homogeneous medium, the signal attenuation is e−σb, where σ is the intrinsic diffusion coefficient.

2.2. Fitting The Adc From The Dmri Signal

An important quantity that can be derived from the dMRI signal is the “Apparent Diffusion Co- efficient” (ADC), which gives an indication of the root mean squared distance travelled by water molecules in the gradient direction g/∥g∥, averaged over all starting positions:

B=0

.

Log S(B) = C0 + C1B + · · · + Cnbn,

increasing n from 1 onwards until we get the value of c1 to be stable within a numerical tolerance.

7

2.3.

Hadc Model

In a previous work , a PDE model for the time-dependent ADC was obtained starting from the Bloch-Torrey equation, using homogenization techniques. In the case of negligible water exchange between compartments (low permeability), there is no coupling between the compartments, at least to the quadratic order in g, which is the ADC term. The ADC in compartment Ωis given by

(11)

is a quantity related to the directional gradient of a function ω that is the solution of the homoge- neous diffusion equation with Neumann boundary condition and zero initial condition:

(12)

n being the outward normal and t ∈[0, TE], ug is the unit gradient direction. The above set of equations, (10)-(12), comprise the homogenized model that we call the HADC model.

2.4. Short Diffusion Time Approximation Of The Adc

A well-known formula for the ADC in the short diffusion time regime is the following short time

Where A

V is the surface to volume ratio and σ is the intrinsic diffusivity coefficient. In the above formula the pulse duration δ is assumed to be very small compared to ∆. A recent correction to the above formula , taking into account the finite pulse duration δ and the gradient direction

!

.

√

∆.

3. Method

Below is a chart describing the work flow of SpinDoctor.

Bend And Twist The Fe Mesh Nodes

by analytical transformation.

Plot Hadc

Figure 1: Flow chart describing the work flow of SpinDoctor The physical units of the quantities in the input files for SpinDoctor are shown in Table 1, in particular, the length is in µm and the time is in µs. Below we discuss the various components of SpinDoctor in more detail.

3.1. Read Cells Parameters

The user provides an input file for the cell parameters, in the format described in Table 2.

(Μsµm)−1

Table 1: Physical units of the quantities in the input files for SpinDoctor.

Height Of Cylinders

Table 2: Input file containing cells parameters.

3.2. Create Cells (Canonical Configuration)

SpinDoctor supports the placement of a group of non-overlapping cells in close vicinity to each other. There are two proposed configurations, one composed of spheres, the other composed of cylinders. The algorithm is described in Algorithm 1.

Algorithm 1: Placing ncell non-overlapping cells. Generate a large number of possible cell centers. Compute the minimum distance, dist, between the current center and previously accepted cells.

Find the intersection of [dist −dmax × Rmean, dist −dmin × Rmean] and [Rmin, Rmax],

2

. If the intersection is not empty, then take the middle of the intersection as the new radius and accept the new center. Otherwise, reject the center. Loop through the possible centers until get ncell accepted cells.

3.3. Plot Cells

SpinDoctor provides a routine to plot the cells to see if the configuration is acceptable (see Fig. 2). Figure 2: SpinDoctor plots cells in the canonical configuration.

3.4. Read Simulation Domain Parameters

The user provides an input file for the simulation domain parameters, in the format described in Table 3.

3.5. Create Surface Triangulation

Finite element mesh generation software requires a good surface triangulation. This means the surface triangulation needs to be water-tight and does not self-intersect. How closely these require- ments are met in floating point arithmetic has a direct impact on the quality of the finite element mesh generated.

It is often difficult to produce a good surface triangulation for arbitrary geometries.

Thus, We

restrict the allowed shapes to cylinders and spheres. Below in Algorithms 2 and 3 we describe how to obtain a surface triangulation for spherical cells with nucleus, cylindrical cells with myelin layer, and the ECS (box or tightly wrapped). We describe a canonical configuration where the cylinders are placed parallel to the z-axis. More general shapes are obtained from the canonical configuration by coordinate transformation in a later step.

Path To Tetgen Cmd

Table 3: Input file of simulation domain parameters.

3.6. Plot Surface Triangulation

SpinDoctor provides a routine to plot the surface triangulation (see Fig. 3).

3.7. Finite Element Mesh Generation

SpinDoctor calls Tetgen , an external package (executable files are included in the toolbox package), to create a tetrehedra finite elements mesh from the surface triangulation generated by

12

Algorithm 2: Surface triangulation of spherical cells and ECS. Suppose we have ncell spherical cells with nucleus. Denote a sphere with center c and radius R by S(c, R), we use the built-in functions (convex hull, delaunnay triangulation) in MATLAB to get its surface triangulation, T(c, R). Call the radii of the nucleus r1, · · · , rncell and the radii of the cells R1, · · · , Rncell. Then the boundaries between the

{Σi = T(Ci, Ri)}, I = 1, · · · , Ncell;

For the box ECS, we find the coordinate limits of the set

S(Ci, Ri) ∈[X0, Xf] × [Y0, Yf] × [Z0, Zf]

and add a gap k = ECS gap × max{xf −x0, yf −y0, zf −z0} to make a box B = [x0 −k, xf + k] × [y0 −k, yf + k] × [z0 −k, zf + k]. We put 2 triangles on each face of B to make a surface triangulation Ψ with 12 triangles.

For the tight-wrap ECS, we increase the cell radius by a gap size and take the union

2

. We use the alphaShape function in MATLAB to find a surface triangulation Ψ that contains W. Algorithms 2 and 3. The FE mesh is generated on the canonical configuration. The numbering of the compartments and boundaries used by SpinDoctor are given in Tables 4 and 5. The labels are related to the values of the intrinsic diffusion coefficient, the initial spin density, and the perme- ability requested by the user. Then the FE mesh nodes are deformed analytically by a coordinate transformation, described in Algorithm 4.

3.8. Plot Fe Mesh

SpinDoctor provides a routine to plot the FE mesh (see Fig. 4 for cylinders and ECS that have been bent and twisted).

3.9. Read Experimental Parameters

The user provides an input file for the simulation experimental parameters, in the format described in Table 6.

13

Algorithm 3: Surface triangulation of cylindrical cells and ECS. Suppose we have ncell cylindrical cells with a myelin layer, all with height H. Denote a disk with center c and radius R by D(c, R), and the circle with the same center and radius by C(c, R). Let the radii of the axons be r1, · · · , rncell and the radii of the cells be R1, · · · , Rncell, meaning the thickness of the myelin layer is Ri −ri.

The boundary between the axon and the myelin layer is:

C(Ci, Ri) × [−H/2, H/2]

We discretize C(ci, ri) as a polygon P(ci, ri) and place one at z = −H/2 and one at z = H/2. Then we connect the corresponding vertices of P(ci, ri) × {−H/2} and P(ci, ri) × {H/2} and add a diagonal on each panel to get a surface triangulation Γi.

Between the myelin layer and the ECS we discretize C(ci, Ri) as a polygon and place one at z = −H/2 and one at z = H/2 to get a surface triangulation Σi. For the box ECS, we find the coordinate limits of the union of D(ci, Ri) and add a gap to make a rectangle in two dimensions. Then we place the rectangle at z = −H/2 and at z = H/2 to get a box. Finally, the box is given a surface triangulation with 12 triangles.

For tight-wrap ECS, we increase the cell radius by a gap size and take the union

I

D(ci, Ri + kRmean). We use the alphaShape function in MATLAB to find a two dimensional polygon Q that contains W. We place Q at z = −H/2 and at z = H/2 and connect correponding vertices, adding a diagonal on each panel. Suppose Q is a polygon with n vertices, then the surface triangulation of the side of the ECS will have 2n triangles.

The above procedure produces a surface triangulation for the boundaries that are parallel to z-axis. We now must close the top and bottom. The top and bottom boundaries is just the interior of Q. However, the surface triangulation cannot be done on Q directly. We must cut out D(ci, ri), the disk which touches the axon, and Ai = D(ci, Ri) −D(ci, ri), the annulus which touches the myelin. Then we triangulate Q −S

I Ai Using The

MATLAB built-in function that triangulates a polygon with holes to get the boundary that touches the ECS. The surface triangulation for Ai and D(ci, ri) are straightforward.

3.10. Btpde

The spatial discretization of the BTPDE is based on a finite element method where interface (ghost) elements are used to impose the permeable interface conditions. The time stepping is done using the MATLAB built-in ODE routine ode23t. See Algorithm 5.

3.11.

Hadc Model

Similarly, the DE of the HADC model is discretized by finite elements. See Algorithm 6.

14

Figure 3: SpinDoctor plots the surface triangulation of the canonical configuration. Left: spherical cells with ECS; Right: cylindrical cells with ECS.

2Ncell + 1

Table 4: The labels and numbers of compartments.

3.12. Some Important Output Quantities

In Table 7 we list some useful quantities that are the outputs of SpinDoctor. The braces in the ”Size” column denote MATLAB cell data structure and the brackets denote MATLAB matrix data structure.

4Ncell + 1

Table 5: The labels and numbers of boundaries.

4. Spindoctor Examples

In this section we show some prototypical examples using the available functionalities of SpinDoctor. 4.1. Comparison of BTPDE and HADC with Short Time Approximation In Fig. 5 we show that both BTPDE and HADC solutions match the STA values at short diffusion times for cylindrical cells (compartments 1 to 5). We also show that for the ECS (compartment 6), the STA is too low, because it does not account for the fact that spins in the ECS can diffuse around several cylinders. This also shows that when the interfaces are impermeable, the BTPDE ADC and that from the HADC model are identical. The diffusion-encoding sequence here is cosine OGSE with 6 periods.

4.2. Permeable Membranes

In Fig. 6 we show the effect of permeability: the BTPDE model includes permeable membranes (κ = 1×10−3 m/s) whereas the HADC has impermeable membranes. We see in the permeable case, the ADC in the spheres are higher than in the impermeable case, whereas the ECS show reduced

16

Algorithm 4: Bending and twisting of the FE mesh of the canonical configuration. The external package Tetgen generates the finite element mesh that keeps track of the different compartments and the interfaces between them. The mesh is saved in several text files.

The connectivity matrices of the finite elements and facets are not modified by the coordinates transformation described below. The nodes are transformed by bending and twisting as described next.

The set of FE mesh nodes {xi, yi, zi} are transformed in the following ways: Twisting around the z-axis with a user-chosen twisting parameter αtwist is defined by



. Bending on the x −z plane with a user-chosen bending parameter αbend is defined by



. Given [αbend, αtwist], bending is performed after twisting. ADC because the faster diffusing spins in the ECS are allowed to moved into the slowly diffusing spherical cells. We note that in the permeable case, the ADC in each compartment is obtained by using the fitting formula involving the logarithm of the dMRI signal, and we defined the ”signal” in a compartment as the total magnetization in that compartment at TE, which is just the integral of the solution of the BTPDE in that compartment.

4.3. Myelin Layer

In Fig. 7 we show the diffusion in cylindrical cells, the myelin layer, and the ECS. The ADC is higher in the myelin layer than in the cells, because for spins in the myelin layer diffusion occurs in the tangential direction (around the circle). At longer diffusion times, the ADC of both the myelin layer and the cells becomes very low. The ADC is the highest in the ECS, because the diffusion distance can be longer than the diameter of a cell, since the diffusing spins can move around multiple cells.

4.4. Twisting And Bending

In Fig. 8 we show the effect of bending and twising in cylindrical cells in multiple gradient directions. The HADC is obtained in 20 directions uniformly distributed in the sphere. We used spherical harmonics interpolation to interpolate the HADC in the entire sphere.

Then We Deformed The

radius of the unit sphere to be proportional to the interpolated HADC and plotted the 3D shape. The color axis also indicates the value of the interpolated HADC.

17

Figure 4: FE mesh of cylinders and ECS after bending and twisting. Compartment number is 1 to 8 for the cylinders and 9 for the ECS.

4.5. Timing

In Table 8 we give the average computational times for solving the BTPDE and the HADC. All simulations were performed on a laptop computer with the processor Intel(R) Core(TM) i5-4210U 2 axons and a tight wrap ECS, the simulated sequence is PGSE (δ = 2.5ms, ∆= 5ms).

In

the impermeable case, the compartments are uncoupled, and the computational times are given separately for each compartment. In the permeable membrane case, the compartments are coupled, and the computational time is for the coupled system (relevant to the BTPDE only).

5. Numerical Validation Of Spindoctor

In this section, we validate SpinDoctor by comparing SpinDoctor with the Matrix Formalism method [49, 50] in a simple geometry. The Matrix Formalism method is a closed form representation of the dMRI signal based on the eigenfunctions of the Laplace operator subject to homogeneous Neumann boundary conditions. These eigenfunctions are available in explicit form for elementary

Tributed Uniformly On A Sphere;

if ngdir = 1, take the gradient direction from the

Depending On Line 13;

Table 6: Input file for simulation experiment parameters. geometries such as the line segment, the disk, and the sphere . The dMRI signal obtained using the Matrix Formalism method will be considered the reference solution in this section.

The accuracy of the SpinDoctor simulations can be tuned using three simulation parameters:

1. Htetgen Controls The Finite Element Mesh Size;

(a) Htetgen = −1 means the FE mesh size is determined automatically by the internal algorithm of Tetgen to ensure a good quality mesh (subject to the constraint that the radius to edge ratio of tetrahedra is no larger than 2.0).

(b) Htetgen = h requests a desired FE mesh tetrahedra height of h µm (in later versions of Tetgen, this parameter has been changed to the desired volume of the tetrahedra).

19

Algorithm 5: BTPDE. FE matrices are generated for each compartment by the finite element method with continuous piecewise linear basis functions (known as P1). The basis functions are denoted as ϕk for k = 1, . . , Nv, where Nv denotes the number of mesh nodes (vertices). All matrices are sparse matrices. M and S are known in the FEM literature as mass and

Ω

σi ∇ϕi · ∇ϕj dx. J has a similar form as the mass matrix but it is scaled with the coefficient g · x, we

Ω

g · x ϕiϕj dx.

W Φiϕj Ds

where a scalar function w is used as an interface marker. The matrices are assembled from local element matrices and the assembly process is based on vectorized routines of , which replace expensive loops over elements by operations with 3-dimensional arrays. All local elements matrices in the assembly of S, M, J are evaluated at once and stored in a full matrix of size 4 × 4 × Ne, where Ne denotes the number of tetrahedral elements. The assembly of Q is even simpler; all local matrices are stored in a full matrix of size 3 × 3 × nbe, where nbe denotes the number of boundary triangles.

Double nodes are placed at the interfaces between compartments connected by permeable membrane. Q is used to impose the interface conditions and it is associated with the interface (ghost) elements. Specifically, assume that the double nodes are defined in a pair

−Q¯I¯J

if vertex i and j belong to two different interfaces The fully coupled linear system has the following form

(14)

where ξ is the approximation of the magnetization M. SpinDoctor calls MATLAB built-in ODE routine ode23t to solve the semi-discretized system of equations. 2. rtol controls the accuracy of the ODE solve. It is the relative residual tolerance at all points of the FE mesh at each time step of the ODE solve;

20

Algorithm 6: HADC model. Eq. (12) can be discretized similarly as described for the BTPDE and has the matrix form

(15)

where ζ is the approximation of w and ¯ζi = σi F(t) ug · n(xi). We note that the matrices here are assembled and solved separately for each compartment. SpinDoctor calls MATLAB built-in ODE routine ode23t to solve the semi-discretized equation.

Integral Of Magnetization At Te

in each compartment.

Integral Of Magnetization At Te

summed over all compartments.

[Ncmpt × Nexperi ]

ADC in each compartment.

Adc Accounting For All Compart-

ments.

Adc In Each Compartment In

each direction.

Adc Accounting For All Compart-

ments in each direction. Table 7: Some important SpinDoctor output quantities.

N/A

Table 8: Computational times for solving the BTPDE and the HADC. All simulations were performed on Intel(R) 2 axons and a tight wrap ECS, the simulated sequence is PGSE (δ = 2.5ms, ∆= 5ms). 3. atol controls the accuracy of the ODE solve. It is the absolute residual tolerance at all points of the FE mesh at each time step of the ODE solve; We varied the finite element mesh size and the ODE solve accuracy of SpinDoctor and ran 6 simulations with the following simulation parameters: SpinD Simul 5-1: rtol = 10−3, atol = 10−6, Htetgen = −1;

Sta Experi 1

Figure 5: Geometry: 5 cylinders, tight wrap ECS, ECS gap = 0.2, ug = [1, 1, 1], σout = σecs = 2 × 10−3 mm2/s, κ = 0 m/s, OGSE cosine (δ = 14ms, ∆= 14ms, number of periods = 6). The vertical bars indicate the ADC in each compartment. The ADC in the rightmost position is the ADC that takes into account the diffusion in all the compartments.

SpinD Simul 5-2: rtol = 10−6, atol = 10−9, Htetgen = −1; SpinD Simul 5-3: rtol = 10−9, atol = 10−12, Htetgen = −1; SpinD Simul 5-4: rtol = 10−3, atol = 10−6, Htetgen = 1; SpinD Simul 5-5: rtol = 10−6, atol = 10−9, Htetgen = 1; SpinD Simul 5-6: rtol = 10−9, atol = 10−12, Htetgen = 1;

Hadc Experi 1

Figure 6: Geometry: 3 spheres, tight wrap ECS, ECS gap = 0.3, ug = [1, 1, 0], σin = σecs = 2 × 10−3 mm2/s, κ = 1 × 10−3 m/s (left), κ = 0 m/s (right). PGSE (δ = 5ms, ∆= 5ms). The vertical bars indicate the ADC in each compartment. The ADC in the rightmost position is the ADC that takes into account the diffusion in all the compartments.

• 3LayerCylinder is a 3-layer cylindrical geometry of height 1µm and the layer radii, R1 = 2.5µm, R2 = 5µm and R3 = 10µm. The middle layer is subject to permeable interface conditions on both the interior and the exterior interfaces, with permeability coefficient κ.

The exterior boundary R = R3 is subject to impermeable boundary conditions. The top and bottom boundaries are also subject to impermeable boundary conditions. • For this geometry, Htetgen = −1 gives finite elements mesh size (nnodes = 440, nelem = 1397).

Htetgen = 1 gives finite elements mesh size (nnodes = 718, nelem = 2088). The dMRI experimental parameters are the following:

23

• the diffusion coefficient in all compartments is 2 × 10−3 mm2/s; • the diffusion-encoding sequence is PGSE (δ = 10ms, ∆= 13ms); • 8 b-values: b = {0, 100, 500, 1000, 2000, 3000, 6000, 10000} s/mm2;

• 1 Gradient Direction: [1, 1, 0];

In Figure 9 we show the signal differences (in percent) of the reference Matrix Formalism method and the SpinDoctor simulations, normalized by the reference signal at b = 0:

Smf (B = 0)

× 100.

(16)

We see that the signal difference is less than 0.35% for κ = 10−5 m/s and it is less than 0.25% for κ = 10−4 m/s for all 6 SpinDoctor simulations. The signal difference becomes smaller when the ODE solve tolerances are changed from (rtol = 10−3, atol = 10−6) to (rtol = 10−6, atol = 10−9), but there is no change when the tolerances are further reduced to (rtol = 10−9, atol = 10−12). If we refine the FE mesh, but keep the ODE solve tolerances the same, the signal difference is in fact larger using the refined mesh than using the coarse mesh at the smaller b-values, though this effect disappears at higher b-values and larger permeability. This is probably due to parasitic oscillatory modes on the finer mesh that need smaller time steps to be sufficiently damped.

6. Computational time and comparison with Monte-Carlo simulation In this section, we compare SpinDoctor with Monte-Carlo simulation using the publicly available software package Camino Diffusion MRI Toolkit , downloaded from http://cmic.cs.ucl.ac.

MATLAB R2019a on the same computer. We give SpinDoctor computational times for three relatively complicated geometries. We also give Camino computational times for the first two geometries. We did not use Camino for the third The number of the degrees of freedom in the SpinDoctor simulations is the finite element mesh size (the number of nodes and the number of elements). For Camino it is the number of spins. The time stepping choice of the SpinDoctor simulations is given by the ODE solve tolerances. For Camino it is given by the number of time steps. Camino has an initialization step where it places the spins and we give the time of this initialization step separately from the Camino random walk simulation time.

Given the interest of the dMRI community in the extra-cellular space and neuron simulations,

We Chose The Following Three Geometries:

1. ECS400axons. See Figure 10. This models the extra-cellular space outside of 400 axons. We generated 400 cylinders with height 1µm and radii ranging from 2 −5µm, randomly placed according to Algorithm 1. The small height of the cylinders means that this geometry should

Deployment In Out-Of-Position Situations

D. Bendjaballah1, A. Bouchoucha1, M. L. Sahli1,2* and J-C. Gelin2

Abstract

Side-impact collisions represent the second greatest cause of fatality in motor vehicle accidents. Side-impact airbags have been installed in recent model year vehicle due to its effectiveness in reducing passengers’ injuries and fatality rates. In meeting these requirements, simulations of folding and deploying airbags are very useful and are widely used. The paper presents a simulation method for the deploying airbags using three materials in different working conditions. Finite element analysis is primarily used to evaluate this concept. In these simulations, the gas flow is described by the conservation laws of mass, momentum, and energy. The numerical results indicate that the FE method in this paper is capable of capturing airbag deploying process accurately.

Keywords: Airbag simulations, Out-of-position, Crash, Modeling, Out-of-position

Background

The passive safety of cars has become a very high prior- ity issue for the automotive industry. Today, there are not only one or two airbags in a car; certain models have ten times more than that. With the increasing usage of airbags, the number of accidents where the airbag itself can cause an injury to the occupant also increases

(Augenstein Et Al. 2003; Gabauer And Gabler 2010;

Audrey et al. 2011). As is well known, safety belts are also now devices designed to provide protection to the users of vehicles during crash events, minimizing the loads necessary to adapt their movement to the move- ment of the car (Freesmeier and Butler 1999; Schmitt et al. 1997). In general, the seat belt is designed to restrain the occupant in the vehicle and prevent the

Occupant From Having Harsh Contacts With Interior

surfaces of the vehicles. The airbag acts to cushion any impact with vehicle structure and has positive internal pressure, which can exert distributed restraining forces over the head and face. As a safety component of auto- mobile, an airbag decreases occupants’ injury likelihood effectively in case of an accident (Ruff et al. 2007). These safety elements can reduce the death rates on the roads, and its protection effects have been widely approved (Crandall et al. 2001; Teru and Ishikawa 2003). With computational tools such as finite element methods designed for dynamic contact problems, crashworthiness simulations can now be used with reliable accuracy to evaluate occupant protection in various collision condi- tions with safety metric/parameters such as acceleration, head injury criteria, intrusion distance, intrusion vel- ocity, and neck forces (neck injury risk or whiplash).

Thus, new types of airbag products are being developed to handle different collision scenarios.

Become Standard Equipment On Most New Passenger

vehicles (Braver and Kyrychenko 2004; Teng et al. 2007; Yoganandan et al. 2007). The airbag cushion is com- posed of a woven fabric which is rapidly inflated during a car crash. The airbag dissipates the passenger’s kinetic energy thereby reducing injury through biaxial stretching of the fabric bag and escaping gas through vents. There- fore, the performance of the airbag is greatly influenced by the mechanical properties of the fabric. Generally, air bags are designed to deploy in a crash that is equivalent to a vehicle crashing into a solid wall at 8 to 14 mph.

Air bags most often deploy when a vehicle collides with another vehicle or with a solid object like a tree. There are various types of airbags: frontal, side-impact, and curtain airbags. In general, the passenger side airbags are usually larger than the driver airbags (see Fig. 1).

Besançon, France

© The Author(s). 2017 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Bendjaballah et al. International Journal of Mechanical

Doi 10.1186/S40712-016-0070-2

Extensive studies have shown that the airbag deploy- ment in load cases consists of two occupant loading phases: a punch-out effect where the airbag bursts out of its container with the airbag and airbag module cover accelerating towards the occupant and a second loading phase during which the airbag is taking on its deployed shape and volume (membrane-loading effect). Bankdak et al. (2002) developed an experimental airbag test system to study airbag-occupant interactions during close proximity deployment. The results provided insight for simulating the effect of inflation energy and mass flow on target response. Bedard et al. (2002) found that while left-side (driver-side) impacts accounted for only 13.5% of all crashes, the fatality rate among these

Crashes Was 68.3% In Comparison To Front Impact

(48.3%), right-side impact (31.3%), and rear impact (38.4%). These studies underscore the importance of oc- cupant safety during side-impact collisions. In the last years, the current market requested to reduce the time and cost airbag development. In order to achieve this result, virtual simulations play an important role since they allow to minimize the number of experimental tests (Pei et al. 2013; Cao et al. 2014). Several simulation models of airbag were established (Wang et al. 2007). It is feasible to optimize the parameters of airbag deploy- ment using simulation technology. Experimental and numerical studies have quantified injury risks to close- proximity occupants from deploying side airbags. These studies have focused on the prevention of the most ad- verse effects of airbag deployment (Duma et al. 2003).

Other studies have proposed airbag characteristics to minimize particular biomechanical responses (Haland and Pipkorn 1996). In a more recent study, Marklund and Nilsson (2003) compared deformation patterns with experimental data as well as the computational costs associated with three different airbag deployment simu- lation methods; they concluded that the SPH method is relatively inexpensive and produces incremental deform- ation patterns that compare most closely to the experi- mental results. The process of inflation of an airbag is one of the determining factors in saving lives. The duration from the initial impact of the crash to the full inflation of an airbag is about 40 ms, and during this time, the airbag goes from being in a folded state to a fully inflated state, with a high internal pressure. After achieving this state, the airbag begins to deflate, thus providing a nice cushion for the body impacting it.

Ideally, the person in the crash should come into contact with the airbag at this time. In the present study, a large volume passenger side airbag model is developed to handle different collision scenarios. The main aim is evaluate the performance of deploying of passenger side airbag using finite element methods (FEM).

Materials

The tensile specimens were made in different airbags (P: Peugeot, R: Renault, and VW: Volkswagen) with a length of 200 mm long and a width of 40 mm. Table 1 shows the mechanical properties of the airbag.

Tensile Tests

To determine the mechanical properties of the material of airbag used in the test pieces, tensile tests were performed on Lloyd EZ20 universal testing machine in Constantine. These tests were conducted using rect- angular samples. The axial force and axial displacement acquired during a test are converted into stress and the strain in order to be used for the fabric material model.

The continuous recording of the stress-strain data was performed during both the load and unload phases. A minimum of five samples were made in order to check the repeatability of the measurements. All the data was collected by using a PC-based data acquisition system and analyzed by commercial software. The picture frame test device that is made for this study is shown in Fig. 2.

Fig. 1 a Frontal and side airbags. b Oblique view of facet occupant model in sitting posture following airbag deployment (Lim et al. 2014)

0.150

Bendjaballah et al. International Journal of Mechanical and Materials Engineering (2017) 12:12

Page 2 Of 9

Figure 3 shows the stress-strain relationship of the airbag sample under axial tensile loads. The results are showing a linear increase in extension with the increas- ing stresses. This is an expected output and it confirms with the theoretical behavior of a sample subjected to tensile stress. The rupture strain values for different airbags (R/P/VW) were 0.322, 0.441, and 0.472, respect- ively. The measured elastic parameters (i.e., Young’s modulus E and initial yield strength) and Poisson’s ratio are summarized in Table 2. The tensile tests of the woven fabrics can show differences on mechanical prop- erties because woven fabrics can resist in-plane shear loads once the yarn lock-up angle has been reached. The differences of material property on material direction can affect the shape of fully deployed bag (see Fig. 3b).

Theoretical Background

Numerical simulations of airbags use very complex and techniques such as an orthotropic model to identify the mechanical behaviors during the airbag inflation and the fluid mechanics (gas flow) to describe the inflator gas flow (pressure gradient) and improve the representation of the pressures within the airbag. To model the airbag as an orthotropic model, three material constants have to be provided. Assuming a plane stress condition, the

Ð1Þ

where σ is the normal stress and τ is the shear stress, the subscript refers to the principal material directions, i.e., the fill and warp directions. Also, ε and γ are the strain components. The material elastic constants Qij are

Ð2Þ

where E1 and E2 are the Young’s modulus in the fill and wrap directions and G12 is the shear modulus of the fabric material. νij is the Poisson ratio of the material.

The gas exerts a pressure load on the airbag causing it to expand. This expansion puts the airbag under tensile stress lowering the expansion rate. In this study, heat conduction and heat transfer is not taken into account.

Fig. 2 A photograph of Lloyd EZ20 universal testing Fig. 3 Stress versus strain using Lloyd EZ20 machine for a three different airbags at 0° and 90° and b VW airbag test specimens at

Different Angles

Table 2 Physical and mechanical properties of the airbag

Page 3 Of 9

In the deployment of an airbag, an inflator supplies high velocity gas into an airbag causing it to expand rapidly. The gas inside the airbag is assumed to be ideal, to be of constant entropy, and to satisfy the equation of state:

Ð3Þ

Here p, ρ, and e are respectively the pressure, density, and specific internal energy, and γ is the ratio of the heat capacities of the gas. The gas flow is described by the conservation laws for mass, momentum, and energy that

Ð4Þ

here, V is a volume, A is the boundary of this volume,

N Is The Normal Vector Along The Surface A, And U

denotes the velocity vector in the volume. Applying Bernoulli’s equation in the case of an ideal gas with

Ð5Þ

Here, the subscript ex denotes quantities at the throat of the tube. Furthermore u, p, and ρ denote the quan- tities inside that part of the tube that is supplying mass.

Materials And Boundary Conditions

The airbag system mainly consists of three parts: the airbag itself, the inflator unit, and the crash sensor or diagnostic unit. Thus, to study the behavior of the airbag using FE simulations, we need to have an FE model of the airbag in the folded position. A FE model of the airbag was used to simulate the test condition as shown in Fig. 5. LS-DYNA® material model FABRIC (MAT_34) is used to simulate the airbag material. It is a variation of the layered orthotropic material model. Additionally, in the LS-DYNA® material model, fabric leakage can be accounted for. However, for this CAB material, the leak- age is almost negligible and therefore no leakage is specified. The mechanical properties can be determined from the physical test. Typical material properties for airbag fabrics are taken as given in Chawla et al. (2004a) (Table 3). These properties are used to simulate inflation process of airbag (see Table 1). The car dashboard is modeled as the rectangular thin plate using a MAT_RI-

Gid Material, And The Degrees Of Freedom Are Con-

strained in all the directions. The similar properties of thermoplastic polymer are assigned for contact purposes. The porosity of the fabric is assumed zero. The nitro- gen gas is taken for inflating the airbag. Properties of nitrogen gas and initial bag conditions are shown in Table 4. The example on which we perform the study is a typical passenger side airbag. The geometric de- tails have been measured from a commercially avail- able airbag. The initial state of the airbag is a closed rectangular whose sides are to be finished to 482 × 635 mm2 and is shown in Fig. 4.

Table 3 Material properties of airbag and rigid plate used in FE

–

Table 4 Initial values used for FE simulation of the swelling of

3.33 × 10−4

Fig. 4 The initial airbag geometry in the form of a rectangular Bendjaballah et al. International Journal of Mechanical and Materials Engineering (2017) 12:12

Related Journal Articles & DOI Links

Selected peer-reviewed publications relevant to 12 Lead ECG Acquisition. Click the DOI to access the full paper (may require institutional access).

Why Choose Us?

Bangalore guidance for robotics, Spectre and autonomous systems projects.