Asymptotic expansion of steady-state potential in a high contrast medium with a thin resistive layer

Applied Mathematics and Computation 221 (2013) 48–65

Contents lists available at SciVerse ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Asymptotic expansion of steady-state potential in a high contrast medium with a thin resistive layer Ronan Perrussel a, Clair Poignard b,⇑ a

Laboratoire Plasma et Conversion d’Énergie, CNRS UMR5213, INPT & UPS, Université de Toulouse, F-31071 Toulouse, France Team-Projet MC2, INRIA Bordeaux-Sud-Ouest, Institut de Mathématiques de Bordeaux, CNRS UMR 5251, Université de Bordeaux1, 351 cours de la Libération, 33405 Talence Cedex, France b

a r t i c l e

i n f o

Keywords: Laplace equation Asymptotic expansion Thin layer High contrast media Cell electroporation

a b s t r a c t We study the steady-state potential in a high contrast medium with a resistive thin layer. We provide asymptotic expansions of the potential at any order, whatever small the layer conductivity is. Approximate transmission conditions at any order and the corresponding variational formulation are given. We prove uniform estimates with respect to the thickness of the layer and with respect to its resistivity. The main insight consists in this expansion at any order and in the corresponding non-standard variational formulation whatever small the layer conductivity is. In particular the limit potential is not continuous across the limit surface (when the membrane thickness shrinks to zero) unlike the soft contrast case and its jump is proportional to the electric flux. Numerical simulations illustrate the theoretical results. Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction In the present paper the electro-quasistatic approximation of the Maxwell equations in a high contrast medium with an insulating thin layer is considered. We aim at providing asymptotic expansions at any order with respect to the membrane thickness, whatever small the conductivity of the layer is. The asymptotic for the thickness tending to zero are very different compared with the soft contrast medium [23,24], since the effect of the thin layer appears at the zeroth order in the high contrast case. 1.1. Motivation The distribution of the steady-state potential in a biological cell is important for bio-electromagnetic investigations. A sufficiently large amplitude of the transmembrane potential (TMP), which is the difference of the electric potentials between both sides of the cell membrane, leads to an increase of the membrane permeability [25,30]. Molecules such as bleomycin can then diffuse across the plasma membrane. This phenomenon, called electropermeabilization, has been already used in oncology and holds promises in gene therapy [21,28], justifying precise assessments of the TMP. Since the experimental measurements of the TMP on living cells are limited — mainly due to the membrane thinness, which is a few nanometers thick — a numerical approach is often chosen [25,27]. However, these computations are confronted with the heterogeneous parameters of the biological cells. Therefore in this paper we derive a rigorous asymptotic analysis to tackle these numerical difficulties. ⇑ Corresponding author. E-mail addresses: [email protected] (R. Perrussel), [email protected] (C. Poignard). 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.06.047

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

(a) The domain with a thin layer.

49

(b) The “background” domain.

Fig. 1. Geometry of the problem.

We consider the three-dimensional model of biological cell given by Schwan [14,15] for different frequency ranges. This model considers the cell as a highly heterogeneous medium composed with a thin resistive membrane surrounding a conductive cytoplasm. Electro-quasistatic approximation1 of the Maxwell equations in the time-harmonic regime is studied here. This approximation is usually considered to describe the behavior of a cell submitted to an electric field of frequency smaller than a few giga Hertz. Depending on the frequency, the modulus of the complex conductivity of the thin layer is either very small (for the frequency range under 10 kHz), or of the same order (for the frequency range between 100 kHz and 100 MHz) compared with the membrane thickness d, which is a small parameter. Such singular perturbation problem can be tackled thanks to appropriate numerical method, which is often time and resource consuming [1–29], which prevents the use of such techniques for time-dependent or non linear problems. For such reason, we propose here to approximate the partial differential equation (P.D.E) problem set in the domain with thin layer by a sequence of P.D.E problems set in a domain without thin layer: the effect of the thin membrane is then translated in the transmission conditions. At any order k 2 N of accuracy dk , we aim at giving an asymptotic expansion of the potential, that is valid whatever the frequency is, and that avoids meshing the thin layer. More precisely, we provide uniform approximate transmission condition, that describes the effect of the layer without meshing it. This uniform estimate over the proposed range of frequency (up to 100 MHz) enables us to apply safely the transmission conditions also for time-transient and non-linear problem, which are more relevant for modeling the electropermeabilization process. The reader may refer to the works of Neu and Krassowska [22] and DeBruin and Krassowska [11] for detailed description of theirnon-linear model of electroporation. Note also that Kavian et al. have recently proposed a simplified non linear model of electroporation [19]. The results of the present paper can also be seen as rigorous justification, and extension of the so-called Kedem–Katchalsky equation obtained in 1958 to describe the transport of molecules across a membrane [20]. We provide some clues to tackle numerically these problems in the concluding remarks at Section 7. 1.2. The studied problem The geometry of the problem is given in Fig. 1. We denote by X a bounded domain with smooth boundary @ X. Let Odc be a subdomain of X surrounded by a thin layer Odm with thickness d. We assume that the domain Odm [ Odc is independent on d and that its distance to @ X is strictly positive. We denote by Oe ¼ X n Odm [ Odc . Observe that Oe is also d-independent. Notation 1.1. Present now the notations used throughout the paper.  We generically denote by n the normal to a closed smooth surface of R3 outwardly directed from the domain enclosed by the surface (see Fig. 1) to the outer domain.  Let C be a surface of R3 , and let u be a function defined in a tubular neighborhood of C. We define ujC by

8x 2 C;

ujC ðxÞ ¼ limþ uðx  tnðxÞÞ; t!0

1 The electro-quasistatic approximation consists in neglecting the curl part of the electric field, which therefore derives from a so-called electric potential. This amounts to considering the steady-state potential equations with complex coefficients.

50

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

moreover if u is differentiable, we define @ n ujC by

8x 2 C;

@ n ujC ðxÞ ¼ limþ ruðx  tnðxÞÞ  nðxÞ; t!0

where  denotes the Euclidean scalar product of R3 . In addition we define the jump ½uC by

½uC ¼ ujCþ  ujC : Let rc be the inner complex conductivity2 of Odc and we denote similarly re and rm the respective conductivities of Oe and Odm . We suppose that both imaginary and real parts of re ; rc , and rm are positive. Let Oc be the smooth bounded domain defined by Oc ¼ X n Oe and denote by C its boundary. We define the domain conductivity by

8 d > < rc ; in Oc ; r ¼ rm ; in Odm ; > : re ; in Oe ;

~¼ and r



rc ; in Oc ; re ; in Oe :

Throughout the paper the characteristic length of X as well as the characteristic conductivity of the domain are assumed equal to 1 so that we only deal with dimensionless quantities, however for the sake of simplicity we omit the term ‘‘dimensionless’’ that should be in front of each physical quantity. The electro-quasistatic formulation is given by

(

r  ðrrud Þ ¼ 0; ud j@ X ¼ g;

in X

ð1Þ

on @ X;

g is the electric potential imposed on the boundary of X. We suppose that g is as regular as we need. We aim at providing the rigorous asymptotic expansion of the potential u for d tending to zero and jrm j 6 d. Remark 1.2. Multiply (1) by ud and integrate by parts. There exists a constant C such that

kud kH1 ðOe Þ 6 CjgjH1=2 ð@ XÞ ;

krud kL2 ðOdc Þ 6 CjgjH1=2 ð@ XÞ

and

C krud kL2 ðOdm Þ 6 pffiffiffiffiffiffiffi jgjH1=2 ð@ XÞ :

rm

ð2Þ

According to the last inequality, the more resistive the thin layer is, the worse the estimate of the electric field rud in the membrane is since rm is small for a resistive medium. Estimate (2) makes us think that ‘‘the electric field is trapped into the membrane’’. Therefore, unlike the case of a soft contrast medium, the resistive thin layer has an influence on the zeroth approximation of the potential. 1.3. State of the art and plan of the paper The main insight of the paper is to consider that the thin layer is highly insulating in the sense that

rm ¼ dnSm ;

ð3Þ

where jSm j ¼ 1 and 0 < jnj 6 1. Depending on the frequency of the studied phenomenon, the parameter jnj, which is a nondimension frequency, can be of the order 1, d, or even smaller. Observe that Sm can be seen as a dimensionless surface conductivity. For each frequency n ¼ dq , with q 2 N, it is possible to provide a precise asymptotic expansion by solving specific partial differential equations. However, in this paper we aim at giving a rigorous uniform variational formulation, that allows to approach the steady-state potential accurately and for all the different cases in a same manner. Several papers in the mathematical research area are devoted to domains with thin layer. In [9] (see also the numerous works of Ammari et al. [3,7,2,4–6]), Capdeboscq and Vogelius provide a general representation of the conductivity in domains with soft contrast small inclusions. They prove that the potential can be approached by the ‘‘background’’ potential — the potential in the domain without inclusion — with the help of the so-called polarization tensor. Their result proves that for a soft contrast small inclusion the influence of the inhomogeneity is the same as a dipole characterized by the polarization tensor (see Theorem 1, page 161 of [9]). More recently, Vogelius and Nguyen [31] show that this general representation formula holds for high contrast small diameter inclusion. They also mention that the representation formula does not hold in 2

The complex conductivity r of a given material is defined by

r ¼ ixe þ s; where s and e denote respectively the (real) conductivity and permittivity of the material, and x is the frequency of the time-harmonic field.

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

51

the case of high contrast thin inclusion. The aim of the present paper is to provide a precise description of the potential for such thin high resistive domains. The method used to derive the effective transmission conditions has been extensively described previously in the case of soft contrast media. We refer to [13] for the original paper and to [24,23] for a more general description. It consists in a suitable change of variables in the thin layer in order to make the small parameter appear in the equations. Recently Schmidt and Tordeux provided the asymptotic expansion of the solution to a bidimensional eddy-current like problem for thin conductive layers [26]. Note also the work of Chun et al. [10] that consider transmission conditions for Helmholtz equation in a soft contrast medium. Our problem is different since the high contrast appears in the transmission conditions satisfied by the potential (that is not the case of Schmidt and Tordeux since the high contrast holds on the zeroth order operator of their second order partial differential equation). In addition, unlike the work of Chun et al. [10], we write the approximate transmission conditions in the limit background domain, with only one surface C (see Fig. 1(b)) whereas Chun et al. [10] consider the domain with the two curves Cd and C (see Fig. 1(a)): this leads to a transmission condition that involves the curvature at the first order, whereas the curvature terms disappear at the first order when dealing with only one curve (but of course they appear at the higher orders). Moreover we consider the three-dimensional resistive case without any assumption on the magnitude of the layer conductivity. More precisely, we only suppose that jrm j is smaller than the layer thickness d whereas Schmidt and Tordeux choose their conductivity proportional to 1=d. In this sense, our result is more general for the resistive case. Another difference lies in the fact that the solution to our limit problem does not belong to the Sobolev space H1 but it is only piecewise H1 (denoted by PH1 ), and therefore an appropriate analysis have to be performed. Despite these differences we let the reader observe that in our result as in the result of Schmidt and Tordeux, the layer influence appears at the zeroth order term, which highlights the main insight of high contrast thin domains. The outline of the paper is the following. In Section 2, we perform the change of variables, that will be useful to derive the asymptotic expansion. We then perform the formal identification of the coefficients in Section 3, and we provide preliminary estimates in Section 4. Section 5 is devoted to the proof of our main result, and in the last section, numerical simulations with the finite element method illustrate the theoretical results. We end by concluding remarks that present two applications of the results provided in this paper: the dynamic electric potential in biological cells, and the magnetostatic field in carbon nanotubes trapped in the wall of foam bubbles. In the following subsections we present the model obtained by electrical engineers by physical considerations, and then we present our main result, that justify and extend the electrical engineering results. 1.4. Resistive thin layer in electrical engineering Several papers in the electrical engineering research area have studied the electric potential in domains with thin resistive layers. We refer the reader to the survey of Bossavit [8] — and references therein — for the modelling principle. More recently, Pucihar et al. have proposed in [25] an electrical model of biological cells where the electric potential uapprox satisfies

8 approx ¼ 0; in Oe [ Oc ; uapprox j@ X ¼ g; > < Du approx Sm ½ u C  re @ n uapprox jCþ ¼ 0; > : e approx ½r@nu C ¼ 0;

ð4Þ

where Sm is the surface membrane conductivity with n ¼ 1 in (3). Observe that in this model the membrane influence appears in the transmission condition. We aim at justifying the model of Pucihar et al. and at providing more accurate approximations of ud by performing a rigorous asymptotic expansion of the potential at any order, that are valid from the resistive case to the very resistive case. An interesting result of the paper consists in the variational formulation of problem (4). More precisely, in their paper Pucihar et al. propose to solve problem (4) using an iterative method, that consists in making explicit the Neumann trace of the potential to compute the jump of the potential at the following iteration. For such a method several computations are needed before reaching an accurate value of the potential. This can be time – and resource – consuming, especially when studying a large frequency range, or when dealing with non-linear problems (in the electroporation modelling, the membrane conductivity depends on the TMP [18]). To avoid this drawback we propose a variational formulation that allows to solve the problem in only one computation. In addition, this formulation is useful in the justification of the asymptotics. 1.4.1. The functional space PH1 ðXÞ Denote by PH1 ðXÞ the space of functions in X, which have the H1 – regularity in Oe and in Oc (but not in X) and let PH10 ðXÞ be the subspace of functions belonging to PH1 ðXÞ that vanish on @ X:

PH1 ðXÞ ¼ fv :

v jO

PH10 ðXÞ ¼ fv :

v 2 PH1 ðXÞ; v j@X ¼ 0g:

e

2 H1 ðOe Þ;

v jO

c

2 H1 ðOc Þg;

ð5Þ ð6Þ

52

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

For any function g 2 H1=2 ð@ XÞ we define the space g þ PH10 ðXÞ by

g þ PH10 ðXÞ ¼

n

v 2 PH1 ðXÞ : v j@X ¼ g

o

ð7Þ

:

1.4.2. Variational formulation of problem (4) The variational formulation of problem (4) satisfied by uapprox writes

(

Find uapprox 2 g þ PH10 ðXÞ such that for all / 2 PH10 ðXÞ R R approx e  r/ d x þ n C Sm ½uapprox C ½/C ds ¼ 0: X r ru

1.5. Main result Before stating our main result we present the change of variables used in order to derive the asymptotic expansion. 1.5.1. Geometry We suppose that the inner boundary C of the domain Oe is a smooth compact 2D – manifold of R3 without boundary. Let xT ¼ ðx1 ; x2 Þ be a system of local coordinates on C ¼ fwðxT Þg. Define the map U by

8ðxT ; x3 Þ 2 C  R;

UðxT ; x3 Þ ¼ wðxT Þ þ x3 nðxT Þ;

where n is the normal vector of C defined in Notation 1.1. The thin layer Odm is then parameterized by

Odm ¼ fUðxT ; x3 Þ;

ðxT ; x3 Þ 2 C  ðd; 0Þg:

The Euclidean metric in ðxT ; x3 Þ is given by the 3  3 – matrix ðg ij Þi;j¼1;2;3 where g ij ¼ h@ i U; @ j Ui:

8a 2 f1; 2g;

g 33 ¼ 1; 2

g a3 ¼ g 3a ¼ 0; 0

8ða; bÞ 2 f1; 2g ;

g ab ðxT ; x3 Þ ¼ g ab ðxT Þ þ 2x3 bab ðxT Þ þ

g 0ab ¼ h@ a w; @ b wi;

bab ¼ h@ a n; @ b wi;

ð8aÞ x23 cab ðxT Þ;

ð8bÞ

where

cab ¼ h@ a n; @ b ni:

ð8cÞ

ij

We denote by ðg Þ the inverse matrix of ðg ij Þ, and by g the determinant of ðg ij Þ. For all l 2 N define

 pffiffi  8 @ i ð g g ij Þ  l > l  > pffiffi a ¼ @ ; for ði; jÞ 2 f1; 2; 3g2 ; < ij 3  g x3 ¼0  > > : Alab ¼ @ l3 ðg ab Þ ; for ða; bÞ 2 f1; 2g2

ð9Þ

x3 ¼0

l

and let S C be the differential operator on C of order 2 defined by

S lC ¼

X

alab @ b þ Alab @ a @ b :

ð10Þ

a;b¼1;2

For all q P 0 the differential operators RqN and RqD are defined by the following recurrence formulae:

R0N ¼ 0;

R0D ¼ Id; R1N ¼ Id; R1D ¼ 0

and for all p P 0

Rpþ2 ¼ N ¼ Rpþ2 D

p   X l C lp al33 Rpþ1l þ Rpl N N SC ;

ð11bÞ

l¼0 p X

  l C lp al33 Rpþ1l þ Rpl D D SC ;

ð11cÞ

l¼0

where C lp ¼ p!=ðl!ðp  lÞ!Þ. 1.5.2. Asymptotic at any order The main result of the present paper holds in the following Theorem, that provides asymptotic expansion of the steadystate potential at any order. Theorem 1.3 (Main result). Let g belong to H5=2þn ð@ XÞ, for n 2 N. The integer n denotes the order of the approximation. Let the four sequences ðg pD Þ06p6n ; ðg pN Þ06p6n ; ðup Þ26p6n and ðupm Þ26p6n be defined by induction:

u2 ¼ u1 ¼ 0; in X;

1 u2 m ¼ um ¼ 0; in ð1; 0Þ  T;

g 0D ¼ g 0N ¼ 0; on C:

ð12aÞ

53

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

For p P 0 the function up is defined in X by

Dup ¼ 0;

in Oe [ Oc ;

up j@X ¼ dp0 g; where dp0 is the classical Kronecker symbol;

ð12bÞ

nSm ½up C  re @ n up jCþ ¼ g pD ;

ð12cÞ

e @ n up C ¼ g pN ; ½r

ð12dÞ

the profile function

upm

upm ¼ up jCþ þ g

is defined in the cylinder C  ð1; 0Þ by

re nSm

  @ n up 

þ Cþ

Z g

0

ðg  sÞ a033 @ g up1 m þ

! p2 l  X s s l p2l p2l alþ1 @ u þ S u ds g m C m l! l þ 1 33 l¼0

ð12eÞ

and for p P 1 the function g pD and g pN are defined on C by p  X 1 l RN ð@ n upl jC Þ þ RlD ðupl jC Þ þ ð1Þp nSm l! l¼1 ! Z 0 p2 l  X g g lþ1 p2l l p2l 0 p1  ð1 þ gÞ a33 @ g um þ a @ g um þ S C um dg: l! l þ 1 33 1 l¼0

g pD ¼ nSm

g pN ¼ rc

p X 1 l¼1



Z

l!

 pl pl Rlþ1 jC Þ þ Rlþ1 jC Þ  ð1Þp nSm N ð@ n u D ðu

0

1

ð13aÞ

a033 @ g up1 m þ

! p2 l  X g g lþ1 p2l p2l dg: a33 @ g um þ S lC um l! l þ 1 l¼0

ð13bÞ

The above functions are uniquely determined and satisfy the following regularity for 0 6 p 6 n

g pD 2 H5=2þnp ðCÞ;

g pN 2 H3=2þnp ðCÞ;

ð14Þ

3þnp

up jOe 2 H ðOe Þ; up jOc 2 H4þnp ðOc Þ;   upm 2 C1 ð1; 0Þ; H5=2þnp ðCÞ :

ð15Þ ð16Þ

Denote by

unapp ¼

n X ðdÞp up ; p¼0

then there exists a constant C n independent on d and n such that

ku  unapp kH1 ðOe Þ 6 C n ndnþ1 jgjH5=2þn ð@ XÞ

ð17Þ

and for all x compactly embedded in Oc

ku  unapp kH1 ðxÞ 6 C n dnþ1 jgjH5=2þn ð@ XÞ ;

ð18Þ

krðu  unapp ÞkL2 ðxÞ 6 C n ndnþ1 jgjH5=2þn ð@ XÞ : Observe that in the elementary problem at the order p, the source terms the elementary problems at the order q for q ¼ 0; . . . ; p  1.

ð19Þ g pD

and

g pN

depend on the solutions u and uqm of q

For n ¼ 0, the approximate potential is the same as the potential of Pucihar et al., therefore our result is a rigorous generalization of the electrical engineering models of resistive thin layers. We provide now the two first order terms of the expansion. 1.5.3. Approximate transmission conditions at the zeroth andat the first orders To obtain the two first order terms, we first observe that according to definition (10), S 0C is the Laplace–Beltrami operator on C, therefore we denote it by S 0C ¼ DC . Moreover the mean curvature H of the surface C equals

H¼

a033 : 2

Hence, applying recurrence formula (12), we get the two first orders of the asymptotics. 1.5.4. The order 0

8 0 in Oe [ Oc ; > < Du ¼ 0; e @ n u0 C ¼ 0; nSm ½u0 C  re @ n u0 jCþ ¼ 0: ½r > : 0 u j@ X ¼ g:

ð20Þ

54

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

In the rescaled membrane ð1; 0Þ  C

u0m ¼ g½u0 C þ u0 jCþ : 1.5.5. The order 1

8 Du1 ¼ 0; in Oe [ Oc ; > > > > > e @ n u1 C ¼ rc DC u0 jC ; < ½r

    > nSm ½u1 C  re @ n u1 jCþ ¼ nSm 1  H nSrcm @ n u0   : > > C > > : 1 u j@ X ¼ 0:

ð21Þ

In the rescaled membrane ð1; 0Þ  C 1 um 1 ¼ u j Cþ þ g

re nSm

  @ n u1 

Cþ

 g2 H½u0 C :

Observe that the effect of the geometry appears at the first order term, unlike the case of a soft contrast medium with a thin layer where it appears at the second order term [23,24]. If we were considering an infinite cylinder, and still denoting by C the curve of the cross section, a033 would be the curvature of C. The transmission condition of (20) is nothing but the so-called Kedem–Katchalsky equation [20], while the transmission condition (21) is the first order correction of this condition. Therefore, Theorem 1.3 provides a generalization at any order of the Kedem–Katchalsky equation. We recall that the zeroth order approximation is u0 and the first order approximation is u0  du1 . In Section 6 we illustrate numerically the theoretical results by studying the rate of convergence of the errors u  u0 and u  ðu0  du1 Þ in H1 – norm in the outer medium Oe . We prove now the asymptotic expansion. 2. Preliminary results 2.1. Laplace operator for functions The Laplace–Beltrami operator acting on functions for the metric ðg ij Þi;j¼1;2;3 is given by3 [12]

1 X pffiffiffi ij

Dg ¼ pffiffiffi @i gg @j : g i;j¼1;2;3

ð22Þ

According to (8b) and using the smoothness and the compacity of C, the functions g ij are analytic on the cylinder C  ½d; 0, for d small enough. Hence the functions alij and Alab defined by (9) are smooth functions defined on C. We infer the expansion of the Laplacian in local coordinates:

Dg ¼ @ 23 þ

 X xl  3 al33 @ 3 þ S lC ; l! lP0

8ðxT ; x3 Þ 2 C  ðd; 0Þ;

ð23Þ

where the differential operators S lC are given by (10). Performing the change of variable g ¼ x3 =d we derive

Dg ¼

 X l gl  g 1 0 2 l lþ1 @ þ a @ þ d a @ þ S g g g C ; d 33 l! l þ 1 33 d2 lP0 1

8ðxT ; gÞ 2 C  ð1; 0Þ:

ð24Þ

The following proposition will be useful in order to derive the asymptotics. Proposition 2.1. Let u be a smooth harmonic function in Oc . We still denote by u the function u  U in the tubular neighborhood of C. For all p 2 N there exists two surfacic operators on C denoted by RpN and RpD respectively of order p  1 and p, and defined by (11c) such that

@ p3 ujC ¼ RpN ð@ n ujC Þ þ RpD ðujC Þ:

ð25Þ

Proof. The proof is an application of the Leibniz rule for the differentiation applied to (23). h The above proposition leads to the following corollary by applying the Taylor formula with integral remainder term in the variable g.

3

@ i denotes the partial derivative with respect to the variable xi , for i ¼ 1; 2; 3.

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

55

Corollary 2.2. Let u be a smooth harmonic function in Oc , and denote by Cd the boundary of Odc . There exists a sequence of operators ðT p ÞpP0 of order p such that for all p 2 N

ujC ¼ d

p X ðdÞk  k¼0

@ n ujC ¼ d

k!

 ðdÞpþ1 RkN ð@ n ujC Þ þ RkD ðujC Þ þ T pþ1 ðuÞ: ðp þ 1Þ!

p  ðdÞpþ1 X ðdÞk  kþ1 RN ð@ n ujC Þ þ Rkþ1 T pþ2 ðuÞ; D ðujC Þ þ k! ðp þ 1Þ! k¼0

ð26Þ ð27Þ

where there exists a constant C such that for all p P 0

jT pþi ðuÞjH1=2i ðCÞ 6 CkukHpþ2 ðOc Þ ;

i ¼ 0; 1:

We are now ready to derive formally the asymptotic expansion.

3. Formal asymptotics Define the so-called profile function um in C  ð1; 0Þ by

8ðxT ; gÞ 2 C  ð1; 0Þ; Since

um ðxT ; gÞ ¼ u  UðxT ; dgÞ:

rm =d ¼ nSm , the transmission conditions inherent to (1) write now

rc @ n ujCd  Uð; dÞ ¼ nSm @ g um jg¼1 ; ujCd  Uð; dÞ ¼ um jg¼1 ; re @ n ujCþ  w ¼ nSm @ g um jg¼0 ; ujCþ  w ¼ um jg¼0 :

ð28Þ ð29Þ

To obtain our transmission conditions, we set

(

u ¼ u0  du1 þ d2 u2     ; in Odc [ Oe ; um ¼ u0m þ du1m þ d2 u2m þ    ; in C  ð1; 0Þ

ð30Þ

and we perform the identification of the terms with the same power in d. Observe that the formal series for u is in ðdÞk whereas the formal series of the profile terms involves dk . We choose this convention for the sake of simplicity since the dk appears naturally in the equality (24), while ðdÞk appears in 26,27. P We emphasize that the coefficients uj will be defined in the whole domain X even if pj¼0 ðdÞj uj approximates u only in d Oe [ Oc . Observe that Eqs. 26,27are written on Cd , while we want to write transmission conditions on C. Insert the ansatz (30) into equalities 26,27 of Corollary 2.2 to derive

ujCd ¼ @ n ujC

d

! p  X X 1 l ðdÞp up jC þ RN ð@ n upl jC Þ þ RlD ðupl jC Þ ; l! pP0 l¼1

! p  X X 1  lþ1 p lþ1 p pl pl ¼ ðdÞ @ n u jC þ R ð@ n u jC Þ þ RD ðu jC Þ : l! N pP0 l¼1

ð31Þ ð32Þ

3.1. Recurrence formulae Let us now identify the terms with the same power in d. We necessarily have

Dup ¼ 0;

u0 j@X ¼ dp0 g:

in Oe [ Oc ;

Remind that dp0 is the Kronecker symbol equal to 1 if p ¼ 0 and 0 otherwise. Using (24), write now the equations satisfied by upm in C  ð1; 0Þ:

 p2 l  X g g lþ1 p2l p2l : a33 @ g um þ S lC um l! l þ 1 l¼0

@ 2g upm ¼ a033 @ g up1 m 

According to transmission conditions (28) and (29), by integrating relations (33) and (32), we infer

! p  X 1  lþ1 pl Þ RN ð@ n upl jC Þ þ Rlþ1 ðu j C D l! l¼1 ! p2 l  X g g lþ1 p2l l p2l 0 p1 dg þ ð1Þp re @ n up jCþ : a33 @ g um þ a @ g um þ S C um l! l þ 1 33 l¼0

nSm @ g upm jg¼1 ¼ ð1Þp rc @ n upC þ ¼ nSm

Z

0

1

ð33Þ

56

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

This leads to the following Neumann transmission condition for up : p  X 1  lþ1 pl RN ð@ n upl jC Þ þ Rlþ1 jC Þ  ð1Þp nSm D ðu l! l¼1 ! Z 0 p2 l  X g g lþ1 p2l l p2l 0 p1 dg:  a33 @ g um þ a @ g um þ S C um l! l þ 1 33 1 l¼0

e @ n up C ¼ rc ½r

To obtain the Dirichlet transmission condition, we use the Taylor equality

8g 2 ð1; 0Þ;

upm ð; gÞ ¼ upm jg¼0 þ g@ g upm jg¼0 þ

Z

g

0

ðg  sÞ@ 2g upm ð; sÞ ds:

According to (31)–(33)–(28)–(29) and applying the above equality with g ¼ 1, we infer p  X 1 l RN ð@ n upl jC Þ þ RlD ðupl jC Þ þ ð1Þp nSm l! l¼1 ! Z 0 p2 l  X g g lþ1 p2l l p2l  ð1 þ gÞ a033 @ g up1 þ a @ u þ S u dg: g m m C m l! l þ 1 33 1 l¼0

nSm ½up C  re @ n up jCþ ¼ nSm

Therefore, we have formally obtain the asymptotic coefficients of the electric potential given by Theorem 1.3. 4. Preliminary estimates Observe that the elementary problems of the expansion are non standard since they involve particular transmission conditions: according to Eq. (12c) the potential is not continuous but it is given by the flux itself and the parameter n. Since we aim at deriving uniform estimates with respect to n ! 0, in this section we present estimates for a generic problem similar to the elementary problems. Let s P 0. Let g; g N , and g D belong respectively to H1=2þs ð@ XÞ; H1=2þs ðCÞ, and to H1=2þs ðCÞ. Let v satisfy

Dv ¼ 0; in Oe [ Oc ; v j@ X ¼ g; e @ n v C ¼ ng N ; ½r   re n½v C  @ n v  ¼ ng D : Sm Cþ

ð34aÞ ð34bÞ ð34cÞ

4.1. Variational formulation The variational formulation of Problem (34c) is

n

Find v 2 g þ PH10 ðXÞ such that for all / 2 PH10 ðXÞ

R X

re rv r/dx þ n

R

C Sm ½

v C ½/C ds ¼ n

R

R

C Sm g D ½/C ds þ n C g N /jC

ds:

ð35Þ Observing that the bilinear form að; Þ defined by

h

i2

8ð/; wÞ 2 PH1 ðXÞ ;

að/; wÞ ¼

Z

r/rw dx þ

Z C

X

½/C ½wC ds

1

is an hermitian product on PH ðXÞ we infer the following proposition. Proposition 4.1. Let s P 0. Let g; g N , and g D belong respectively to H1=2þs ð@ XÞ; H1=2þs ðCÞ, and to H1=2þs ðCÞ. There exists a unique solution v to problem (34c), which belongs to PH1 ðXÞ. More precisely, using elliptic regularity

v jO

e

2 H1þs ðOe Þ;

v jO

and

c

2 H2þs ðOc Þ:

4.2. Uniform estimates in n for the generic problem In this section, we aim at providing uniform estimates for n tending to zero, for the solution v to (34c). Let I be defined by



DI ¼ 0; in Oe ; I j@X ¼ 0; IjC ¼ 1:

Multiply (34c) by I in Oe and by 1 in Oc , and integrate by parts twice in Oe , and once in Oc . We obtain the following equality:

Z

C

re v jCþ @ n I jCþ ds ¼

Z

@X

re g@ n I ds þ n

Z

C

g N ds:

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

57

Using the above equality, and multiplying (34c) by I in Oe and by 0 in Oc integration by parts leads to

Z C

½v C ds ¼

 Z  1 g D þ g N ds: Sm C

We therefore define the sequence ðv le ; v lc ÞlP0 by

(

Dv 0e ¼ 0; in Oe ; v 0e j@X ¼ g; re @ n v 0e jC ¼ 0

and

v 0c ¼

1 mesðCÞ

Z C

v 0e ds 

Z 

gD þ

C

  1 g N ds Sm

and for l P 1

(

Dv le ¼ 0;

in Oe ;

ð36aÞ

l v le j@X ¼ 0; re @ n v le jC ¼ Sm ðv l1  v l1 e c ÞjC  Sm d1 g D

and

8 l > < Dv c ¼ 0; in Oc ; l rc @ n v lc jC ¼ Sm ðv l1  v l1 c ÞjC  d1 ðSm g D  g N Þ; > R l e :R l C v c ds ¼ C v e ds:

ð36bÞ

Proposition 4.2. For all k P 0; ðv ke ; v kc Þ 2 H1þsþk ðOe Þ  H2þsþk ðOc Þ, and which only depends on the domain Oe and Oc such that for all k P 0

v 0c

is constant. Moreover there exists a constant C Oe ;Oc ,

  kv ke kH1þs ðOe Þ þ kv kc kH2þs ðOc Þ 6 C kOe ;Oc jgjH1=2þs ð@ XÞ þ jg D jH1=2þs ðCÞ þ jg N jH1=2þs ðCÞ :

ð37Þ

Proof. According to (36b) there exists a constant C Oe ;Oc such that for all k P 1

kv ekþ1 kH1 ðOe Þ þ kv kþ1 kH1 ðOc Þ 6 C Oe ;Oc ðkv ke kH1 ðOe Þ þ kv kc kH1 ðOc Þ Þ; c hence

  kv ekþ1 kH1 ðOe Þ þ kv kþ1 kH1 ðOc Þ 6 C kOe ;Oc ðkv 1e kH1 ðOe Þ þ kv 1c kH1 ðOc Þ Þ 6 C kþ1 c Oe ;Oc jgjH1=2 ð@ XÞ þ jg D jH1=2 ðCÞ þ jg N jH1=2 ðCÞ : Then we use the elliptic regularity satisfied by

v ke and v kc

to conclude. h

P P Proposition 4.3. For n small enough, the series lP0 nl v le and lP0 nl v lc normally converge in H1þs ðOe Þ and in H1þs ðOc Þ respectively. Moreover, the sums satisfy Problem (34c), hence by uniqueness

v jO

e

¼

X l n v le ;

v jO

lP0

Particularly, since

c

¼

X l n v lc : lP0

v

0 c

is constant, we infer that all the derivatives of

v jO

c

are of order n.

Proof. Choose n 2 ð0; 1=ð2C Oe ;Oc ÞÞ and Proposition 4.2 leads to Proposition 4.3.

h

4.3. A priori estimates for the asymptotic coefficients The following proposition is an application of both Propositions 4.2 and 4.3. The proof performed using the induction process is left to the reader. Proposition 4.4. Let n 2 N and s P 0 and suppose that g 2 H1=2þsþn ð@ XÞ. The functions defined by (12) and (13) are uniquely determined and satisfy the following regularity for 0 6 p 6 n

g pD 2 H1=2þsþnp ðCÞ;

g pN 2 H1=2þsþnp ðCÞ;

ð38Þ

58

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

up jOe 2 H1þsþnp ðOe Þ;

up jOc 2 H2þsþnp ðOc Þ;

ð39Þ

  upm 2 C1 ð1; 0Þ; H1=2þsþnp ðCÞ :

ð40Þ

In addition there exists a constant C X;n such that for all 0 6 p 6 n the following estimates hold

kupe kH1þsþnp ðOe Þ 6 C X;n jgjH1=2þsþn ð@XÞ ; kupc kH2þsþnp ðOc Þ 6 C X;n jgjH1=2þsþn ð@ XÞ ; max jupm ð; gÞjH1=2þsþnp ðCÞ 6 C X;n jgjH1=2þsþn ð@ XÞ :

g2½1;0

Moreover, we have the following estimate uniformly in n:

krulc kH1þsþnp ðOc Þ 6 C X;n njgjH1=2þsþn ð@XÞ : Remark 4.5. More precisely, according to (12) for all p P 0 there exists a sequence of functions ðf q Þ06q6p defined on C such that in the cylinder C  ð1; 0Þ, the coefficients upm writes

8ðxT ; gÞ 2 C  ð1; 0Þ;

upm ðxT ; gÞ ¼

p X

gq f q ðxT Þ:

q¼0

The functions ðf q Þ06q6p are defined by induction using (12) and they satisfy the following regularity:

f q 2 H1=2þsþpq ðCÞ;

for q 2 0; . . . ; p:

5. Error estimates For p P 1, we denote by Lp the operator on C  ð1; 0Þ of order 1 in the g-variable and 2 in the xT -variables defined as

Lp ¼ Dg 

1 d2

2

@g þ

a033 d@ g

þ

p X

l

d

l¼2

gl2  g ðl  2Þ! l  1

al1 33 @ g

l2

þ SC

! 

:

For any function u 2 C 1 ð½1; 0; Hsþ3=2 ðCÞÞ

  max jLp ðuÞð; gÞjH1=2þs ðCÞ 6 C p dp1 max juð; gÞjHsþ3=2 ðCÞ þ j@ g uð; gÞjHsþ3=2 ðCÞ :

g2ð1;0Þ

g2½0;1

ð41Þ

According to Proposition 4.3, for all k P 0, the derivatives of uk jOc are of order n. Define U k as follows:

( k

U ¼

u0  du1 þ    þ ðdÞk uk ; in Oe [ Odc ; k m m um 0 þ du1 þ    þ d uk ; in C  ð1; 0Þ:

Our main result provides uniform estimates of the error performed by the approximate expansion at any order, for jnj tending to zero. The following theorem straightforwardly leads to Theorem 1.3. Theorem 5.1. Let k P 0. Suppose that g belongs to H5=2þk ðCÞ. Let wk ¼ u  U k . There exists a constant C k such that

kwk kH1 ðOe Þ 6 C k ndkþ1 jgjH5=2þk ð@ XÞ

ð42Þ

and for all x compactly embedded in Oc

kwk kH1 ðxÞ 6 C k dkþ1 jgjH5=2þk ð@XÞ ; k

kþ1

krw kL2 ðxÞ 6 C k nd

jgjH5=2þk ð@ XÞ :

Proof. By definition of the functions ðul ÞlP0 and ðum l ÞlP0 , we infer

rDwk ¼ 0; in Oe [ Odc ; rm Dg wk ¼ rm Lk ðU k Þ; in C  ð1; 0Þ; wk j@ X ¼ 0;

ð43Þ ð44Þ

59

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

with the following transmission conditions

rc @ n wk jCd  rm @ n wk jCþd ¼ rm @ n U k jCþd  rc @ n U k jCd ; rm @ n wk jC  re @ n wk jCþ ¼ 0; wk jC  wk jCþ ¼ U k jCþ  U kC ; d

d

d

d

wk jCþ  wk jC ¼ 0: Moreover, by definition we have

rm @ n U k jCþd  rc @ n U k jCd ¼ rc dkþ1 T kþ2 ðU k Þ; U k jCþ  U k jC ¼ dkþ1 T kþ1 ðU k Þ: d

d

According to Proposition 4.4 and to Corollary 2.2 we infer

   kþ1þi k  ðU Þ T

H1=2i ðCÞ

6 Cn;

i ¼ 0; 1:

Multiply by / in H1 X and integrate by parts with the help of Green’s formula. Using equality (41) and the metric ðg ij Þ in local coordinates in the layer, by taking / ¼ wk and using the trace theorem on C we infer that

krrwk kL2 ðXÞ ¼ Oðndkþ1=2 Þ: d k ~ k ¼ wk  dkþ1 um ~k Then defining w kþ1 in C  ð1; 0Þ and w ¼ w in Oe [ Oc we infer

~ k kL2 ðXÞ ¼ Oðndkþ1 Þ; krrw hence the theorem. h

6. Numerical simulations for a 3D – resistive thin layer For the numerical purpose, we compare the exact solution to problem (1) to the approximate solution u0  du1 , where u0 and u1 are given by 20,21, in order to illustrate the theoretical convergence rate by numerical simulations. Transmission conditions of 20,21 show that the potential is discontinuous across the membranes, which is consistent with resistive thin layer properties. However these discontinuities have to be taken into account properly to obtain accurate numerical simulations. With the help of physical arguments Pucihar et al. present conditions [25] similar to (20), and Theorem 1.3 generalizes such conditions at any order of accuracy. In order to solve such a problem, Pucihar et al. propose an iterative method that consists in expliciting the flux rc Þ@ n ujC , however this method is numerically instable, time consuming, and can lead to inaccurate results, especially for the time-transient case or when including cytoplasmic constituents. Actually, observe that the jump of the potential is given in terms of its normal derivative, and therefore this cannot be rigorously considered as a source term imposed on the surface C. We propose here the variational method coming from (35) that leads to a direct computation of the involved transmission conditions. We suppose that the space-discretization of X is already performed. Since the numerical difficulties lie in the transmission conditions – the jumps of the potentials and the fluxes across the membrane – we present how to compute the rigidity matrix in the finite element method (F.E.M). 6.1. Rigidity matrix of the static problem (20) Ne

Nc

Denote by ð/ke Þk¼1 , and ð/kc Þk¼1 the bases of finite dimensional subspaces respectively of H1 ðOe Þ and H1 ðOc Þ, and extend these functions by 0 in the whole domain X. Let N ¼ N e þ N c . Therefore, the sequence of functions N

ðwk Þk¼1 ¼ ð/1e ; . . . ; /Ne e ; /1c ; . . . /Nc Þ is a basis of a finite dimensional subspace of the space PH1 ðXÞ, which is the space of functions of X that are H1 in Oe , and Oc . In each subdomain Oe , and Oc u can be written as:

u¼

N X

ak wk

k¼0

Table 1 Electric parameters of problem (1).

re

rc

Sm

2

1

1

60

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

Fig. 2. Isovalues of u0  du1 computed by the finite element method with d ¼ 0:05. The Dirichlet condition imposed on the outer circle is gðhÞ ¼ cosðhÞ. The membrane is resistive but not perfectly: the potential is very small in the domain but it is not completely equal to zero.

Fig. 3. Log–log diagram of the error versus the membrane thickness for n ¼ 1.

and according to problem (20) (and equivalently to problem (21)) the rigidity matrix ðMkl Þ is then

Mkl ¼

Z X

re rwk :rwl dx þ

Z

Sm ½wk ½wl  dr:

ð45Þ

C

6.2. 2D simulations by Fourier expansion We first start with bidimensional simulation. Domains X and Oi are two disks centered in 0 of respective radius equal to 2 and 1. Observe that in such a case, the mean curvature of the infinite cylinder equals 1=2, hence H of Eq. (21) simply equals 1/ 2. Parameters re ; rc and Sm are given by Table 1 The Dirichlet boundary condition of (1) is given by

gðhÞ ¼ cosð3hÞ: Analytic solution to the problems (1) is obtained through a calculus in Fourier series. The solution to problems (20) and (21) are computed by a finite elements method. In Fig. 2, we plot the isovalues of the approximate solution u0  du1 , which approaches the solution to (1), with d ¼ 0:05. Observe that the membrane is resistive but not perfectly, thus the potential in the inner domain is small but it is uniformly equal to zero. Fig. 2.

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

61

Fig. 4. Log–log diagram of the error in Oe of the zeroth order for different small conductivities.

Fig. 5. Log–log diagram of the error in Oe of the zeroth order for different small conductivities. Computations are performed with the finite element method.

Fig. 6. Computational domain and the solution for one particular d.

62

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

Fig. 7. Log–log diagram of the error versus the membrane thickness.

Fig. 8. Description of polymer foams with CNT at different scales.

From the analytic solutions the expected rates of convergence with d are verified since the numerical estimated exponents are respectively 1 for the order 0 and 2 for the order 1; see Fig. 3. The green line corresponds to the zeroth order for the soft contrast media. In the soft contrast case, the effect of the layer appears only at the first order term, meaning that the influence of the layer is of order d. We observe that this potential does not provide an accurate approximation of the potential, since the error estimate does not converge, when d tends to zero. As predicted by Theorem 5.1 and Theorem 1.3, the rate of convergence of the zeroth order term u0 defined by (20) increases when the membrane conductivity decreases, justifying the use of such approximate conditions whatever small the membrane conductivity is; see Fig. 4. This behavior is also recovered when using the finite element method based on the variational formulation (35) for solving our model problem; see Fig. 5. We conclude by observing that the convergence rates are asymptotic meaning that if the membrane is not thin enough, it is not sure that our approximate transmission conditions are relevant. For instance, in the configurations we chose, the accuracy of the transmission conditions seems to increase for a thickness around 101 but if we were choosing other configurations, the accuracy could have decreased as well. Therefore, the shape of the error curves are less significant for a thickness larger than 102 . 6.3. 3D simulations by the finite element method For this simple 3D example, the geometry and the source term are supposed to be axisymmetric (independent of the angular variables in a cylindrical coordinate system). It enables to perform the computations on an equivalent 2D problem. The reduced 2D computational domain is given on Fig. 6(a) and the solution ud satisfies the following reduced problem:

   8  @ @ud @ @ud > > < @r rr @r þ @z rr @z ¼ 0; d d > u ¼ 1; on Cu and u ¼ 0; on Cd ; > : @ud ¼ 0; on CN : @n Moreover, for sake of simplicity, n equals d and

ð46Þ

rc , re and Sm are equal to 1.

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

63

As no analytic solution is known for this example, the reference solution is computed by the finite element method. On Fig. 6(b), this solution is given for d ¼ 0:05. The asymptotic convergence in d is obtained for the order 0 and in d2 for the order 1 as shown on Fig. 7. In this particular situation, the pre-asymptotic regime is shorter than in the 2D example. Note the great interest of such approximate condition, compared with the full problem with a thin membrane: Instead of meshing the thin membrane, we just have to mesh the surface of the cell, which is very significant for 3D-simulationst. Moreover the variational formulation (35) makes possible to compute the linear potential in one iteration, unlike the iterative methods of Pucihar et al. or Weaver’s group [16]. 7. Concluding remarks In this paper, we have provided asymptotic expansions of the electric potential for a resistive thin layer, whatever small the membrane conductivity is. The expansion leads to particular elementary problems in the sense that the jump of the potential across the membrane is given by the flux. We have provided a variational formulation that allows to compute the potential without any lifting of the flux, allowing to compute the approximate potential in only one computation whatever small the membrane conductivity. We conclude by pointing out some possible applications of our results. 7.1. Transient electric potential in a cell. Such results are ‘‘important’’ when dealing with time dependent problem. More precisely, the time dependent electric potential in a biological cell satisfies (see part 5 of [17])

DV ¼ 0; in Oe [ Oc ; Vðt; xÞjx2@ X ¼ g; C m @ t ½V C þ Sm ½V C  re @ n VjCþ ¼ 0 e @ n V C ¼ 0: ½r

Vð0; Þ ¼ 0

ð47aÞ ð47bÞ ð47cÞ

Guyomarc’h and Lee propose a discontinuous Galerkin method to compute V by making explicit the flux. The main drawback of their time explicit scheme is that a Courant–Friedrich–Lewy type condition, that links the time-step and the mesh, has to be satisfied leading to time – and resource – consuming computations. Our variational formulation that can be easily computed by discontinuous Galerkin method should lead to unconditionally stable scheme since the flux is implicited. More precisely, denoting by Dt the time-step of the semi-discretization in times, (47b) writes

    Cm Cm t þ Sm ½V tþ1 C  rc @ n V tþ1  ¼ ½V C Dt Dt  C and thus the rigidity matrix is very similar to the static case (45), since Sm has to be replaced by C m =Dt þ Sm :

Mkl ¼

Z X

re rwk :rwl dx þ

Z  C

 Cm þ Sm ½wk ½wl  dr: Dt

7.2. Diffusion of molecules in domains separated by a thin membrane Note that the variational formulation is also interesting when dealing with the diffusion of a molecule across a thin membrane of thickness d. For such problem, let M be the molecule that diffuses in the above domain X, and let De ; Dm and Dc be the coefficients of diffusion of M in the three domain. Assume also that in the thin membrane Om the diffusion coefficient is e be the diffusion coefficient defined as very low, equal to Dm ¼ dP m . Let D

e ¼ D



De ;

in Oe ;

Dc ;

in Oc ¼ X n Oe :

Then it is well-known that the Kedem–Katachalsky equation describes the diffusion of M accross the layer:



@ t M  r:ðDe rM Þ ¼ 0;

in Oe ;

@ t M  r:ðDc rM Þ ¼ 0;

in Oc ;

ð48Þ

with the transmission conditions on C

De @ n MjCþ ¼ Dc @ n MjC ;

Pm ½M  Dc @ n MjC ¼ 0:

Actually, discretizing (48) in time as

  Mtþ1  Dt r: De;c rMtþ1 ¼ Mt ; in Oe;c ; then the rigidity matrix writes

Mkl ¼

Z X

wk :wl dx þ Dt

Z X

e rwk :rwl dx þ D D

Z C

Pm ½wk ½wl  dr:

ð49Þ

64

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

7.3. Magnetic potential in composite materials. Another important application of our present paper holds for magnetic properties of composite materials composed of foams with carbon nanotubes (CNTs) trapped into the wall of bubbles (see Fig. 8). The magnetic permeability of the CNTs is very high (depending on the concentration of CNT) compared with the air permeability l0 , but the thickness of the bubble wall is thin compared with the bubble diameter. The domain X is composed of a domain O0 with a network Od containing the CNTs. In a 2D approximation of the problem, the vector magnetic potential Ad (reduced to a scalar in 2D) satisfies:

 r:



1

ld



rAd ¼ Jsource ;

in X

Adj@ X ¼ 0: Due to the high contrast in the permeabilities and since the thickness of the CNT network is small, the computation of the magnetostatic field is time- and resource- consuming. Our result shows that, denoting the CNT permeability by

ld jXd ¼ d1 lm ; the magnetic potential Ad can be approached at the order d by A0 , the variational solution in PH10 ðXÞ to:

1

l0

Z

rA0 r/ dx þ

X

l0 lm

Z C

½A0 ½/ ds ¼

Z X

Jsource / dx;

8/ 2 PH10 ðXÞ;

where C is the bidimensional surface corresponding to the network of CNT. Acknowledgments The authors thank warmly P. Dular and R. Sabariego from the Electrical Engineering Department of the University of Liège, and L. Krähenbühl from the Laboratoire Ampère of Lyon for helpful discussions for the motivations of the work presented here. The authors are also very grateful to K. Schmidt from DFG Research Center for well-considered advices and suggestions. C.P. is granted by the projects of the French National Agency ‘‘projet blanc’’ ANR INTCELL and ANR MEMOVE. References [1] G.M. Amiraliyev, The convergence of a finite difference method on layer-adapted mesh for a singularly perturbed system, Appl. Math. Comput. 162 (3) (2005) 1023–1034. [2] H. Ammari, E. Beretta, E. Francini, Reconstruction of thin conductivity imperfections, II. The case of multiple segments, Appl. Anal. 85 (1–3) (2006) 87– 105. [3] H. Ammari, S. He, Effective impedance boundary conditions for an inhomogeneous thin layer on a curved metallic surface, IEEE Trans. Antennas and Propagation 46 (5) (1998) 710–715. [4] H. Ammari, H. Kang, Reconstruction of conductivity inhomogeneities of small diameter via boundary measurements, in: Inverse Problems and Spectral Theory, Contemp. Math., vol. 348, Amer. Math. Soc., Providence, RI, 2004, pp. 23–32. [5] H. Ammari, H. Kang, Reconstruction of Small Inhomogeneities from Boundary Measurements, Lecture Notes in Mathematics, vol. 1846, SpringerVerlag, Berlin, 2004. [6] H. Ammari, H. Kang, Polarization and Moment Tensors, Applied Mathematical Sciences, vol. 162, Springer, New York, 2007. [7] H. Ammari, S. Moskow, Asymptotic expansions for eigenvalues in the presence of small inhomogeneities, Math. Methods Appl. Sci. 26 (1) (2003) 67– 75. [8] A. Bossavit, Modèles et modélisation en électrotechnique, Techniques de l’Ingénieur A1207 (1991) 1–25. [9] Y. Capdeboscq, M.S. Vogelius, A general representation formula for boundary voltage perturbations caused by internal conductivity inhomogeneities of low volume fraction, M2AN Math. Model. Numer. Anal. 37 (1) (2003) 159–173. [10] S. Chun, H. Haddar, J.S. Hesthaven, High-order accurate thin layer approximations for time-domain electromagnetics, Part II: transmission layers, J. Comput. Appl. Math. 234 (8) (2010) 2587–2608. [11] K. DeBruin, W. Krassowska, Modelling electroporation in a single cell I Effects of field strength and rest potential, Biophys. J. 77 (1999) 1213–1224. [12] B.A. Dubrovin, A.T. Fomenko, S.P. Novikov, Modern geometry—methods and applications. Part I, Graduate Texts in Mathematics, second ed., vol. 93, Springer-Verlag, New York, 1992. The geometry of surfaces, transformation groups, and fields, Translated from the Russian by Robert G. Burns.. [13] B. Engquist, J.C. Nédélec, Effective boundary condition for acoustic and electromagnetic scattering thin layer. Technical Report of CMAP, 278, 1993. [14] E.C. Fear, M.A. Stuchly, Modelling assemblies of biological cells exposed to electric fields, IEEE Trans. Biomed. Eng. 45 (10) (1998) 1259–1271. [15] K.R. Foster, H.P. Schwan, Dielectric properties of tissues and biological materials: a critical review, CRC Biomed. Eng. 17 (1) (1989) 25–104. [16] Thiruvallur R. Gowrishankar, James C. Weaver, An approach to electrical modeling of single and multiple cells, Proc. Natl. Acad. Sci. 100 (6) (2003) 3203–3208. [17] Grégory Guyomarc’h, Chang-Ock Lee, Kiwan Jeon, A discontinuous Galerkin method for elliptic interface problems with application to electroporation, Comm. Numer. Methods Eng. 25 (10) (2009) 991–1008. [18] A. Ivorra, J. Villemejane, L.M. Mir, Electrical modeling of the influence of medium conductivity on electroporation, Phys. Chem. Chem. Phys. 12 (34) (2010) 10055–10064. [19] O. Kavian, M. Leguèbe, C. Poignard, L. Weynans, Classical electropermeabilization modeling at the cell scale, J. Math. Biol. (2012). . [20] O. Kedem, A. Katchalsky, Thermodynamic analysis of the permeability of biological membranes to non-electrolytes, Biochim. Biophys. Acta 27 (1958) 229–246. [21] M. Marty, G. Sersa, J.-R. Garbay, et al, Electrochemotherapy – an easy, highly effective and safe treatment of cutaneous and subcutaneous metastases: results of esope (european standard operating procedures of electrochemotherapy) study, E.J.C Supplements 4 (2006) 3–13. [22] J. Neu, W. Krassowska, Asymptotic model of electroporation, Phys. Rev. E 53 (3) (1999) 3471–3482. [23] C. Poignard, Approximate transmission conditions through a weakly oscillating thin layer, Math. Methods Appl. Sci. 32 (4) (2009) 435–453.

R. Perrussel, C. Poignard / Applied Mathematics and Computation 221 (2013) 48–65

65

[24] C. Poignard, P. Dular, R. Perrussel, L. Krähenbühl, L. Nicolas, M. Schatzman, Approximate conditions replacing thin layer, IEEE Trans. Mag. 44 (6) (2008) 1154–1157. [25] G. Pucihar, T. Kotnik, B. Valicˇ, D. Miklavcˇicˇ, Numerical determination of transmembrane voltage induced on irregularly shaped cells, Ann. Biomed. Eng. 34 (4) (2006) 642–652. [26] K. Schmidt, S. Tordeux, Asymptotic modelling of conductive thin sheets, Zeitschrift fr Angewandte Mathematik und Physik (ZAMP) 61 (4) (2010) 603– 626. [27] J.L. Sebastián, S. Muñoz, M. Sancho, J.M. Miranda, Analysis of the influence of the cell geometry and cell proximity effects on the electric field distribution from direct rf exposure, Phys. Med. Biol. 46 (2001) 213–225 (electronic). [28] G. Serša, Application of electroporation in electrochemotherapy of tumors, in: Electroporation based technologies and treatment: Proceedings of the International Scientific Workshop and Postgraduate Course, pages 42–45, 14-20 November 2005. Ljubljana, SLOVENIA. [29] Clara Sussmann, Dan Givoli, Yakov Benveniste, Combined asymptotic finite-element modeling of thin layers for scalar elliptic problems, Comput. Methods Appl. Mech. Eng. 200 (47–48) (2011) 3255–3269. [30] J. Teissié, M. Golzio, M.P. Rols, Mechanisms of cell membrane electropermeabilization: a minireview of our present (lack of?) knowledge, Biochimica et Biophysica Acta 1724 (2005) 270–280. [31] M.S. Vogelius, H-M. Nguyen, A representation formula for the voltage perturbations caused by diametrically small conductivity inhomogeneities. Proof of uniform validity, Ann. I. H. Poincaré 26 (2009) 2283–2315.