Baptiste Percier

Thibaut Divoux

Nicolas Taberlet

Insights on the local dynamics induced by thermal cycling in granular matter

Keywords: 45, 70, Cc - 47, 57, Gc - 83, 80, Fg -

In this letter, we report results on the effect of temperature variations on a granular assembly through Molecular Dynamic simulations of a 2D granular column. Periodic dilations of the disks are shown to perfectly mimic the behavior of grains submitted to thermal cycling, and to provide a simultaneous description of the compaction process and the local grains dynamics while accounting for finite size effects. Through extensive simulations, we show quantitatively that the individual disks properties, namely their roughness and elastic modulus define a critical cycling amplitude of temperature ∆Tc. Below ∆Tc, the dynamics is intermittent and spatially heterogeneous while confined into localized regions, whereas above ∆Tc, the whole column flows and the disks dynamics ranges continuously from cage-like at the bottom of the column to purely diffusive at the top. Our results rationalize a series of experimental results reported in literature on dry or immersed granular materials and provide a simple simulation method for the study of thermal cycling as an alternate driving method for soft glassy materials.

Flows of soft glassy materials (SGM) such as foams [1,2], emulsions [3], colloidal glasses [4] or granular media [5,6] exhibit complex dynamics whose common features include aging, spatially inhomogeneous flows (shearbanding) [7,8], non-local effects in confined geometries [9], etc. Recent breakthroughs in the field have been twofold. First, dense flows of density matched suspensions and dry granular materials have been described successfully as a dynamical critical phenomenon unifying suspension and granular rheology [10]. Second, the Herschel-Bulkley scaling law, which successfuly accounts for the steady state flow of a large number of SGMs has been linked on the one hand to the material microstructure [11,12], and on the other hand to the transient yielding process [13]. However, because of subbtle non-local effects and cooperative behaviors [9,14,15], the inertial approach mentionned above fails to describe flows in the limit of vanishingly small shear rates [6]. Furthermore, for such quasistatic flows transient regimes may become extremely long-lived [13], making irrelevant the use of any constitutive, steady state equations (e.g. the Herschel-Bulkley model) to describe such creep-ing flows. Consequently, quasitatic flows have triggered a large amount of studies either by continuous shear [1,16] or by periodic shear [17][18][19], and still urge for more experimental results.

Among the latter means, periodic variations of temperature, coined thermal cycling, are used to induce slow motion in both thermal [20] and athermal systems with granularity [21][22][23]. On top of being a delicate alternative to other driving methods such as cyclic shear [17,24,25], tapping [19,26] or flow-induced fluidization [27], this method has successfully brought key information on quasistatic flows of granular matter : grains may exhibit stringlike motions [28] hinting towards a glassy-like dynamics, while force chains and thus finite size effects are expected to play a key role [29]. However, to date the individual grain dynamics has received very little attention and its link to the macroscopic control parameters remains an open issue.

In this letter, we report Molecular Dynamic (MD) simulations of a 2D granular column sumitted to periodic variations of temperature. This thermal cycling is modeled through the quasi-static sinusoidal dilation of grains in a fixed container. This purely geometric approach accounts for the experimental features of previous macroscopic studies and allows to quantify the key role of the grains properties. Here we show that there exists a critical cycling amplitude ∆T c fully determined by the grain surface roughness and elastic modulus: for cycling amplitudes lower than ∆T c , the granular pile compacts by jumps associated with localized plastic events taking place over the whole column. Whereas for cycling amplitudes larger than ∆T c , the column flows continuously. In this last case, the finite size of the container sets the limit between two different regions: in the upper part of the column, the grains motion is purely diffusive, whereas in the lower part of the column the grains exhibit a glassy-like dynamics including cage rearrangements at short timescales and a diffusive motion at long timescales. Finally, we show that the diffusion coefficient is inversely proportionnal to the pressure inside the column and scales as a non trivial power law with the cycling amplitude and the grain diameter.

Set-up.-MD simulations of N deformable disks (N = 7500 unless otherwise mentionned) of mean diameter 2R = 1 (20% polydispersity) and densely packed (initial packing fraction: φ 0 = 0.81) were performed in a rectangular container (width L = 30, unless otherwise mentionned). Periodic variations of temperature are mimicked through the periodic dilation of the grains by a relative radius ∆T , ranging from 10 -5 to 10 -2 . Therefore, the driving is here only geometrical and as such our study strongly differs from previous approaches modelling the temperature field and the heat transfer across the grains [29]. Furthermore, the container is rigid and kept fixed during the cycles of temperature. Thus the grains rearrangements only result from their own dilation and not from the dilation of the boundaries [START_REF] Clément | Proceedings of the IIIrd Intern. Conf. on Powders and Grains[END_REF]. Each cycle of temperature is devided into 7.5 × 10 4 timesteps during which the disks trajectories are obtained by integrating the equations of motion using Verlet's method. Grains are submitted to gravity and contact forces. A spring-dashpot scheme is used to model the normal force F n ij between two contacting disks i and j:

F n ij = 2k n (r eff /R)δ + η(R/2r eff )∂ t δ
, where δ denotes the disks overlap; k n is the normal stiffness and r -1

ef f ≡ r -1 i + r -1 j
with r i and r j the disks radii [START_REF] Schäfer | [END_REF]. Unless otherwise mentionned, we choose k n = 10 5 and η = 100 leading to a restitution coefficient e = 0.48, but the results reported in this letter are robust and do not depend on the specific values of these parameters. A standard tangential force model with elasticity was implemented [32]:

F t ij = -sgn(ζ t ) min(k t ζ t , µF n ij ) where ζ t = t t0 v ij (t ′ )dt ′ is the net tangential displacement after contact is first established at time t = t 0 ; v ij
is the relative velocity; µ = 0.3 is the dynamic friction coefficient, and k t = 3.3 × 10 4 the spring stiffness. Consequently, the ratio µF n /k t is much smaller than a disk diameter [START_REF] Schäfer | [END_REF]. Finally, the disks are introduced at the top of the column and left to settle with an additional fluid drag, which is switched off prior to the start of the periodic dilations. Such process produces an initial packing fraction varying from φ 0 = 0.79 to φ 0 = 0.83. We carrefully checked that the following results do not depend on φ 0 within this range.

Global dynamics.-Under periodic dilations the granular column flows in a very similar fashion to what has been previously reported in the literature (fig. 1): the packing fraction slowly increases with the number of applied cycles [21,22]. Depending on the cycling amplitude ∆T , the compaction process exhibits two different regimes: at low cycling amplitudes, the granular column flows by successive jumps [Fig. 1(a)], while the packing fraction remains almost constant. For larger cycling amplitudes, the successive dilations of the disks lead to a slow and continuous compaction of the column which does not reach any steady state, even after n = 1000 cycles [Fig. 1(c)] 1 .

In order to quantify these two regimes, the total increase in the packing fraction ∆φ = φ -φ 0 after n = 1000 cycles was computed for various cycling amplitudes ∆T , and disks properties: different elastic modulii (k n ) and surface properties (k t , and µ). Surprinsingly, the increase in ∆φ with ∆T can be scaled onto a single curve while introducing a critical cycling amplitude ∆T c ≡ mg/R • µ/(k n k t ) which corresponds to the minimum amplitude to fully mobilize the contact between two disks during a cycle of temperature [Fig. 1(b)]. Indeed, at low cycling amplitudes the disks do not slide on each other which reads k t R∆T < µmg following the Amontons-Coulomb law, and the dilationinduced forces do not overcome gravity: k n R∆T < mg. A transition is expected when both disks may slide and the dilation becomes comparable to gravity, which reads: (k n R∆T ).(k t R∆T )= µ(mg) 2 . As a key result, the am- plitude ∆T c , which quantifies the limit between the two regimes described above, is fully determined by the individual grain properties.

Confinement effects.-In order to assess the impact due to the finite size of the system, we compute the average displacement δr achieved by the disks between two successive cycles of temperature along the column height h, for different cycling amplitudes below and above ∆T c [Fig. 2(a)]. Interestingly, the grains displacement is not homogeneous along the vertical, but rather displays two different regimes separated by a characteristic altitude h

(∆T ) c . Below h (∆T ) c
, the grains displacement scales as √ h, and only the prefactor depends on the cycling amplitude ∆T within the explored range, whereas above h

(∆T ) c
, δr increases faster with h. The existence of these two different regimes is fully accounted for by the friction of the disks on the sidewalls. Indeed, the evolution along the column height of the vertical component of the stress tensor σ yy averaged over the column width also exhibits two regimes [Fig. 2

(b)]. It defines another characteristic height h (σ) c

below which the presure is constant, and above which the presure increases strongly with h. The two characteristic heights h are in fact equal [Fig. 2(c)], which illustrates that the screening effect due to the redirection of the stress, known as the Janssen effect [33], is responsible for the different grain displacement behavior along the column height. Furthermore, it is worth emphasizing that the size of the upper and unscreened part of the granular column is proportionnal to the column width L [Fig. 2(d)], and thus controlled by the finite size of the column. Local dynamics.-Let us turn to the individual grain dynamics and start with the case of large cycling amplitudes (∆T > ∆T c ). To correlate the macroscopic compaction process to the local disks dynamics, we focus on the horizontal displacement (the vertical motion being biased by gravity) of the grains during τ time steps: ∆x(τ ) = x(t + τ ) -x(t), and compute its mean square ∆x 2 at different altitudes of the column. ∆x 2 is averaged over horizontal slices of one grain diameter and over the duration of the experiment, except for the first twenty cycles which strongly depend on the preparation protocol of the column [Fig. 3(a)]. At the bottom of the container, the grains exhibit a subdiffusive behavior at short lag times [ ∆x(τ ) 2 ∝ τ α1 , with α 1 < 1] and a diffusive behavior at larger lag times [ ∆x(τ ) 2 ∝ τ α2 , with α 2 ≃ 1]. The plateau behavior of the MSD at the bottom of the container is the signature of cage-like rearrangements [Fig. 3(c)], as first evidenced in [34] and also reported for other cycling methods in granular matter [17,24]. From the bottom to the top of the column, one observes a continuous change of the MSD which tends toward a purely diffusive behavior: the exponent α 1 increases continuously up to 1 [Fig. 3(b)]. We can thus define a diffusion coefficient D ≡ ∆x 2 (τ ) /2τ of the motion induced by thermal cycling, including in the lower part of the column for τ ≥ 100 and which a priori depends on the altitude h inside the column, and the cycling amplitude. We found that D increases along the column height, especially above h c and at a given altitude that D increases with the cyling amplitude [Fig. 4 (inset)]. As a central result, both the influence of ∆T and h on D are captured by the following mastercurve:

D(h, ∆T ) ∼ h L 4 ∆T 4 k n d σ yy (1)
The continuous change of the grain dynamics along the column height is also confirmed by computing the probability distribution function of the disks horizontal displacement P (∆x), between two successive cycles of temperature [Fig. 3(d)]. At the bottom of the column, where the cagelike dynamics takes place, P (∆x) is well described by a gaussian: the particles diffuse inside their cage, whereas in the upper part of the column the grains are more mobile and P (∆x) follows a decreasing exponential scaling [35].

Finally, focusing on the population of disks which exhibits the largest displacements allows to capture the heterogeneous nature of the dynamics that the diffusion coefficient misses. Indeed, for large cycling amplitudes, all the disks move and those presenting the highest displacement are all located in the upper part of the column [Fig. 5 (a) and movie 1 in the supplemental material], whereas decreasing the amplitude leads to a very different scenario [Fig. 5 (b), and movie 2 in the supplemental material]. For ∆T < ∆T c , the grains with the largest displacement are heterogeneously spread inside the column into localized regions which temporal dynamics is intermittent [Fig. 5 (c) and movie 3 in the supplemental material].

Discussion.-The master curve reported in figure 1(c) rationalizes the two types of quasistatic flows reported for dry grains in [22] while varying the cycling amplitude: discontinuous compaction process at low amplitude vs. continuous compaction at large amplitude. Our results also explain why lubricated grains for which the inter-grain friction is reduced, experience a faster compaction dynamics than dry grains [28]. The master curve allows to get one step forward, since the grain elastic modulus had not been varied on previous experimental reports. Indeed, fig. 1(c) demonstrates that the grains surface roughness (µ) and their elastic modulus (k t ) play a key and equivalent role on the flow dynamics through a single control parameter ∆T c . A practical consequence is, for a given cycling amplitude, that smooth and soft grains will be less sensitive to temperature variations than rough and hard grains. One can thus expect temperature variations to trigger mainly localized rearrangements in jammed assemblies of soft objects. Nontheless, these rearrangements are non negligible and contribute to the system dynamics as reported in [20]. A systematic investigation of the effect of temperature variations on dense static packings of emulsions, microgels or other soft colloids is still lacking, and figure 1(c) provides an experimentaly testable prediction.

Compared to previous experimental studies mostly limited to macroscopic observations, our results show that the local grain dynamics is different on each side of the critical cycling amplitude. Above ∆T c , finite size effects split the column into two: the upper part where the grains exhibit a diffusive behavior, and the lower part where the grains exhibit cage-like rearrangements at short lagtimes. The diffusion coefficient displays an unusual linear scaling with the disks diameter, which contrasts with the scaling reported for a 2D Couette shear-cell (d 2 ) [36], and the one observed in the split bottom cell geometry (d 3 ) [37]. The key difference lays in the scale at which the driving takes place: here, each grain is driven individually by periodic dilation as opposed to any other macroscopic shear where the energy is injected at the macroscopic scale and cascades down towards the grain scale. Nonetheless, the diffusion coefficient is, in a robust fashion, inversly proportionnal to the pressure inside the column as observed in other geometries [37]. Below ∆T c , the compaction process is spatially heterogeneous and takes place in localized regions which are strikingly similar to the "hot spots" recently reported in [38]. A more detailed analysis is needed to properly make the link between these localized rearrangements and any non-local effect.

Conclusion.-We have introduced a simple and efficient method to mimic the effect of temperature variations on a dense packing of disks. Focusing on the impact of the first thousand cycles of temperature has allowed to establish the link between the individual grain properties and the compaction dynamics. Of course, a systematic study over a larger amount of cycles is needed. This 2D work also paves the way for future studies on the effect of thermal variations on 3D assemblies of hard grains and soft colloids. Future work in our group will focus on the dynamics of single and multiple intruders of different mechanical properties (density, size, shape, etc) in dense packings submitted to temperature variations.

Fig. 1 :

 1 Fig. 1: Evolution of the packing fraction φ -φ0, where φ0 = 0.81 is the initial packing fraction, with the number n of cycles of temperature for a cycling amplitude ∆T = 10 -5 (a) and ∆T = 2.5 × 10 -3 (b). (c) Packing fraction increase ∆φ after n = 1000 cycles normalized by ∆φp, the value of ∆φ in the limit of low cycling amplitudes vs. the dimensionless cycling amplitude ∆T /∆Tc, where ∆Tc is the critical cycling amplitude (see text). The symbols correspond to disks of different mechanical properties: [bullet: kn = 10 5 , kt = 3300, µ = 0.05, 0.25, 0.5, 0.75, 1]; [+: kn = 10 5 , kt = 100; µ = 0.5];[△: kn = 10 4 , kt = 3.3 × 10 4 , µ = 0.5]. The vertical dashed line separates the two compaction regimes: by jumps at low ∆T and continuous for large ∆T .

Fig. 2 :

 2 Fig. 2: (a) Average displacement δr of the grains along the column height h, for 4 different cyling amplitudes (∆T = 2.5×10 -4 , 5×10 -4 , 10 -3 and 2.5×10 -3 from left to right). (b) Vertical component of the stress tensor σyy averaged over the column width vs. the column height h for ∆T = 2.5 × 10 -3 . In both (a) and (b), the horizontal dashed line defines a characteristic height h (∆T ) c and h (σ) c respectively (see text) which are ploted versus one another in (c) for different parameters (column width: 25 ≤ L ≤ 60, and disks number 5.10 3 ≤ N ≤ 2.10 4 ). The black line has a slope equal to 1. (d) Extent of the upper part of the column H -hc vs. the column width L. The black line is the best linear fit: H -hc = 2.5L.

Fig. 3 :

 3 Fig. 3: (Color online). (a) Mean square value of the disks horizontal displacement ∆x 2 vs. the lagtime τ in cycles unit at different heights of the column [from z = (bottom, red) to z = (top, pink)]. Each MSD is averaged over an horizontal slice of one disk diameter and fitted by a single or two powerlaws of exponents α1 and α2 at short and long lagtimes respectively. (b) MSD exponents α1 (black) and α2 (grey) along the column height h. (c) Exemple of a disk cagelike trajectory located at the center of the column (d) Normalized probability distribution σP of the horizontal grain displacement ∆x between two consecutive cycles of temperature, where σ stands for the standard deviation of the displacements, computed at different heights of the column: z = 20, 40, 105, 192 and 240. Each curve is shifted from the previous one by a factor of 10 for the sake of clarity. The lower (upper resp.) fit is a gaussian (exponential resp.) presented as a guide for the eye. ∆T = 5.10 -3 > ∆Tc. Same color code as in (a).

Fig. 4 :

 4 Fig. 4: (Color online). Diffusion coefficient rescaled by the pressure D(h, ∆T ) × σyy vs. the column height rescaled by the cycling amplitude (h∆T /L) 4 knd. Colors from cyan to red, stand for different cycling amplitudes: from ∆T = 2.5 × 10 -4 , 5 × 10 -4 , 10 -3 , 2.5 × 10 -3 to 5 × 10 -3 . Inset: log-linear plot of the diffusion coefficient D vs. the column height h/d.

Fig. 5 :

 5 Fig. 5: (Color online). Colormap for the whole column of the disks displaying the largest displacement, i.e. larger than 5 times the mean value computed every 5 consecutive cycles of temperature, reported for three different cycling amplitudes: ∆T = 3∆Tc (a), ∆Tc (b) and 0.1 ∆Tc (c). For ∆T > ∆Tc, all the disks move and the ones presenting the largest displacement are located at the top of the column where the pressure is minimum. Whereas, for ∆T < ∆Tc, the disks presenting the largest displacement are heterogeneously distributed in localized regions inside the column.

Note that for even larger cycling amplitudes, convection rolls develop inside the granular pile, while the packing fraction quickly reaches a steady state. This regime, very similar to what has been reported for tapping experiments[39], corresponds to extremely large dilation amplitudes which are out of the scope of the present letter.

We appreciate useful conversations with B. Blanc, J.-C. Géminard, S. Joubaud, M. Leocmach, S. Manneville & L. Ponson.


Marc Bonnet
email: bonnet@athena.polytechnique.fr

Differentiability of strongly singular and hypersingular boundary integral formulations with respect to boundary perturbations 1

Differentiability of strongly singular and hypersingular

boundary integral formulations with respect to boundary perturbations.

Introduction

In e.g. shape design analysis, inverse problems or fracture mechanics, one is often faced with the need of computing sensitivities of integral functional or physical variables with respect to perturbations of the shape of the geometrical domain Ω under study. This goal is achievable by resorting to either finite-difference methods, considering small but finite domain perturbations, or analytical differentiation followed by discretization.

The analytical approach is a priori clearly superior in terms of both accuracy and efficiency. It relies on either the adjoint variable approach or a direct differentiation of the field equations formulated in weak or BIE fashion. A substantial research effort has been devoted in the last decade or so to various formulations and applications of sensitivity analyses based on analytical differentiation with respect to shape parameters, or on the related mathematical concept of domain derivative. As a result, these concepts are successfully applied to more and more engineering problems. Among a fairly abundant literature, the reader is referred to [START_REF] Haug | Design Sensitivity Analysis of Structural Systems[END_REF], [START_REF] Simon | Second variations for domain optimization problems[END_REF], [START_REF] Sokolowski | Introduction to shape optimization[END_REF], [START_REF] Dems | On a class of conservation rules associated with sensitivity analysis in linear elasticity[END_REF].

Further, since problems with variable or unknown domains put a great emphasis on the shape of the boundary ∂Ω of Ω, it is often found convenient, or even essential, to resort to the boundary element method (BEM). Some of these investigations concern the adjoint approach (see e.g. [START_REF] Burczyński | Shape design sensitivity analysis via material derivative -adjoint variable approach[END_REF]), which is not discussed further in this paper. Others formulate the direct differentiation approach, which consists of taking the material derivative of the relevant governing boundary integral equation (BIE), so that a governing BIE for the sensitivities of field variables on the boundary is available. The sensitivity of any objective function of interest is then evaluated using (boundary) field variable sensitivities. Both approaches are the main subject of a recent journal special issue, [START_REF] Bui | Unknown or variable domains, inverse problems[END_REF].

The material differentiation formula applies a priori to nonsingular or at most weakly singular integrals on moving surfaces. On the other hand, usual BIE formulations are either strongly singular or hypersingular. When regularized BIE formulations are considered, the boundary integrals are at most weakly singular and the material differentiation formula can be applied in a straightforward way. Such formulations are studied, up to second-order material derivatives, in [START_REF] Bonnet | Regularized BIE formulations for first-and second-order shape sensitivity of elastic fields[END_REF] for three-dimensional problems and in (Zhang and Mukherjee, 1991) for two-dimensional problems; see also [START_REF] Matsumoto | Optimum design of cooling lines in injection moulds by using boundary element design sensitivity analysis[END_REF] and [START_REF] Nishimura | A boundary integral equation method for an inverse problem related to crack detection[END_REF]. Other authors, e.g. [START_REF] Barone | A boundary element approach for recovery of shape sensitivities in threedimensional elastic solids[END_REF], [START_REF] Erman | Design sensitivity analysis of three-dimensional elastostatic solids using the boundary integral equation method[END_REF], [START_REF] Mellings | Flaw identification using the boundary element method[END_REF], apply the material differentiation formula in an equally straightforward way but to strongly singular or hypersingular BIE formulations, without prior regularization. The mathematical validity of this operation, although unquestioned, is not obvious at first. Indeed, the singular BIE formulations are defined as limiting cases of representation formulas when a vanishing neighbourhood of the collocation (singular) point is removed from the domain under consideration. Then, the vanishing neighbourhood is generally affected by the domain perturbations considered for sensitivity analysis, with a priori possible consequences on the result of the limiting process. There lies the source of possible mathematical difficulties when, as in the above references, the effect of domain perturbation on the vanishing neighbourhood is ignored. In the past, failure to take into account a similar effect led to an erroneous value of the free-term associated with the singular representation of strain in the presence of initial strain or stress, until corrected by [START_REF] Bui | Some remarks about the formulation of three-dimensional thermoelastoplastic problems by integral equations[END_REF].

The goal of this paper is to put material differentiation of usual singular BIE formulation on safe ground. More precisely, we prove that the usual material differentiation formula for nonsingular surface integrals still yields the correct result when applied to BIE-type strongly singular or hypersingular integrals. In order to do so, the singular integrals are expressed in a fixed (i.e. independent of the domain perturbation) parameter space. The direct approach of [START_REF] Guiggiani | A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method[END_REF] and [START_REF] Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF], although initially devised in con-nexion with the mapping of a boundary element onto its parent element, is in fact applicable to any regular surface parametrization and provides the key tool for the present analysis. As a result, the usual, straightforward material derivative technique applies without modification to singular BIE formulations. The material derivative of a singular integral takes the form of another singular integral, with the same level of singularity as the initial one. The present analysis encompasses a wide range of singular kernels, including most usual fundamental solutions used in static or dynamic BIE formulations.

Material differentiation

Let us consider, in the three-dimensional Euclidean space R 3 equipped with a Cartesian orthonormal basis (e 1 , e 2 , e 3 ), a body Ω p whose shape depends on a finite number of shape parameters p = (p 1 , p 2 , . . .). The latter are treated as time-like parameters using a continuum kinematics-type lagrangian description and introducing an "initial" (nonperturbed) configuration Ω = Ω 0 conventionally associated with p = 0:

y ∈ Ω → y p = Φ(y; p) ∈ Ω p where (∀y ∈ Ω) Φ(y; 0) = y (1)
Similarly to continuum kinematics, it is assumed that any such mapping y → Φ(y; p), also termed geometrical transformation, is a diffeomorphism between Ω and Ω p . A given domain perturbation considered as a whole, as is the case e.g. in shape optimization, admits many different representations (1).

As first-order derivatives with respect to p are considered here, attention is focused without loss of generality to the effect of infinitesimal variations δp about p = 0 of a generic shape parameter p while the others are kept fixed and equal to 0. The initial transformation velocity θ(y) is defined by

θ(y) = ∂ ∂p Φ(y; 0) (2)
In other words, the geometrical transformation (1) has the form

y → y δp = y + θ(y)δp + o(δp)
The material derivative f (y) (also termed total or lagrangian derivative) of a generic field variable f (y p , p) in the domain transformation, taken at p = 0, is defined as:

f (y) = lim δp→0 1 δp [f (y δp , δp) -f (y, 0)] = ∂f ∂p + f ,m θ m (y, 0) (3)
using the notation f ,m ≡ ∂f /∂y m . Also, the material derivative of a two-point function K(y px p ) (e.g. the kernel functions that appear in boundary integral equations), assuming both points x, y follow the geometrical transformation (1), is given by

K (x, y) = K ,m (y -x)[θ m (y) -θ m (x)] (4)
The material derivatives of the unit normal n and the differential area dS on any moving generic smooth surface S p changing according to (1) are given [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF], [START_REF] Bonnet | Regularized BIE formulations for first-and second-order shape sensitivity of elastic fields[END_REF] by:

dS= D m θ m dS n i = -n m D i θ m (5)
where

D m f ≡ f ,m -n m n p f ,p = f ,m -n m ∂f ∂n
denotes the projection of the cartesian partial derivative ∂f /∂y m onto the tangent plane at y ∈ S. The numbers D i θ m define in fact the tangential gradient ∇ S θ while the scalar D m θ m is the surface divergence of the vector field θ.

The material derivative of a generic regular integral over S p is then given, from e.g. [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF], by the formula:

d dp Sp f dS = S f +f D m θ m dS (6)
Note that in formulas (3), ( 5), ( 6) and the sequel, all p-derivatives are understood for p = 0. The argument p is omitted for brevity when p = 0.

Weakly singular integral on a changing surface

Let us consider weakly singular integrals of the form

I W (p) = Sp K w (y p -x p )u(y p , p) dS y ( 7
)
where u is regular and possibly depends on p, explicitly and/or implicitly (e.g. u is the field variable which solves a boundary-value problem over the changing domain Ω p ) and K w is a weakly singular fundamental solution. For three-dimensional potential theory or elastostatic, such kernel functions are usually linear combinations of

K w (z) = z i z j |z| 3 (8)
where 1 ≤ i, j ≤ 3 are fixed indexes; besides, they are homogeneous of degree -1 and symmetric:

K w (az) = 1 a K w (z) K w (-z) = K w (z) (9)
In the BIE context, the surface S p is any portion of the (changing) boundary ∂Ω p containing the singular point x p , which is also assumed to follow the domain transformation (1) (this is emphasized by the chosen notation). It follows that the material derivative of the kernel function K w is given by (4). Then, application of formula (6) to the weakly singular integral (7) and of (4) leads to:

d dp I W = S K w ,m (y -x)[θ m (y) -θ m (x)]u(y) dS y + S K w (y -x)[ u +uD m θ m ](y) dS y (10)
It is important to note that, under the smoothness assumption made on the domain transformation, one has:

|θ(y) -θ(x)| ≤ C |y -x| (11)
for some C > 0. As a consequence, the integrand in eq. ( 10) is weakly singular, and the straightforward application of the material derivative formula (6) to I W yields the correct result.

Strongly singular integral on a changing surface

Let us now consider CPV singular integrals of the form

I S (p) = lim ε→0 sε(x p ) K s (y p -x p ) dS y + Sp-eε(x p ) K h (y p -x p )u(y p , p) dS y ≡ - Sp K s (y -x)u(y p , p) dS y (12)
and K s is a strongly singular fundamental solution. For three-dimensional potential theory or elastostatics, such kernel functions are usually linear combinations of

K s (z) = z i z j z k |z| 5 (13)
where 1 ≤ i, j, k ≤ 3 are fixed indexes; besides, they are homogeneous of degree -2 and antisymmetric:

K s (az) = 1 a 2 K s (z) K s (-z) = -K s (z) (14)
When a fundamental solution for flux or traction vector is used in the definition of the integral (12), the unit normal components are included in the nonsingular factor u, in order to be still able to invoke the homogeneity property (14a).

In the BIE context, the surface S p is any portion of the (changing) boundary ∂Ω p containing the singular point x p , which is also assumed to follow the domain transformation (1).

CPV integrals are defined as the result of a limiting process where a spherical exclusion neighbourhood v ε (x) of vanishing radius ε centered at x is removed around x. The difficulty that arises when a changing domain is considered is that the geometrical transformation (1) is likely to alter the shape of v ε (x p ), and this is not allowed by the definition of CPV convergence. Therefore, it is not a priori obvious that a mere application of (6) to (12) yields the correct value for the p-derivative of I S .

To investigate this issue on a firm ground, we introduce a parametrization of the initial surface S on a parameter space ∆:

ξ = (ξ 1 , ξ 2 ) ∈ ∆ → y(ξ) ∈ S (15)
This allows the changing surface S p to be mapped onto the fixed parameter space domain ∆:

(ξ ∈ ∆, p ≥ 0) → y p (ξ) = Φ(y(ξ), p) ∈ S p (16)
so that one can perform the change of variable y → ξ in both I S (0) and I S (δp). Of course, v ε (x p ), or rather e ε (x p ), is also distorted when mapped onto ∆. However, [START_REF] Guiggiani | A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method[END_REF]) solved this difficulty by giving the correct value of the limiting case ε → 0 expressed in the parameter space ∆; their analysis assumes that ∆ is the parent element but in fact remains equally true for any smooth mapping ∆ → S. As a consequence, a rigorous and computable expression of I S using the parametrization ( 16) is available for any p > 0, and indeed provides a sound basis for the present discussion. The latter result, which we now recall, is based upon the introduction of polar coordinates (ρ, α), centered at the image η of the singular point x p in the parameter space ∆

ξ 1 = η 1 + ρ cos α ξ 2 = η 2 + ρ sin α so that ξ ∈ ∆ ⇔ 0 ≤ ρ ≤ ρ(α)
(no loss of generality occurs in assuming that ∆ is star-shaped around η). Since the mapping ( 16) is smooth, the position vector r p = y px p takes the form:

r p = ρr(ρ, α, p) = ρ[a(α, p) + o(1)] (17)
where a(α, p) = r(0, α, p) = a 1 (p) cos α + a 2 (p) sin α (18) and (a 1 , a 2 ) is the natural basis, calculated at ξ = η, associated with the mapping ( 16):

a β (p) = ∂y p ∂ξ β ξ=η
An immediate consequence of definition ( 18) is:

a(α + π, p) = -a(α, p) (19)
The kernel function K s (r p ) then takes the form

K s (r p ) = 1 ρ 2 K s (r) = 1 ρ 2 {K s (a) + O(ρ)} (20)
The CPV integral ( 12) is then given, following [START_REF] Guiggiani | A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method[END_REF], by:

I S (p) = 2π 0 ρ(α) 0 F (ρ, α, p) - f (α, p) ρ dρdα + 2π 0 f (α, p) ln[ρ(α)a(α, p)]dα (21)
using the notations a(α, p) = |a(α, p)| and

F (ρ, α, p) = 1 ρ K s (r)u(y p , p)J(ξ, p) f (α, p) = K s (a)u(x p , p)J(η, p)    (x p = Φ(y(η), p), y p = Φ(y(ξ), p)) (22)
J(ξ, p) being the Jacobian of the parametrization ( 16); a = a(α, p) and r are defined by eqs. ( 17) and ( 18). Note that the last term in (21) accounts for the preservation of the shape of v ε (x). Besides, the antisymmetry property (14b) of K s plays a key role in the validity of the result (21) (I S would otherwise be unbounded when ε → 0) and implies, according to definition (22), that:

f (α + π, p) = -f (α, p) (23)
5 Material derivative of the strongly singular integral

The formula (21) holds for any fixed p ≥ 0, the parameter space ∆ being independent of p; moreover the double integral is regular because of the singularity removal achieved through the introduction of f (α, p). Thus, the derivative of I S (p) at p = 0 is obtained using ordinary differentiation of (21) under the integral sign:

I S = d dp I S = 2π 0 ρ(α) 0 F (ρ, α) - 1 ρ f (α) dρdα + 2π 0 f (α) ln[ρ(α)a(α)]dα - 2π 0 f (α) a (α) a(α) dα (24)
with (using formula (4)):

F (ρ, α) = 1 ρ J(ξ) K s (r)[ u +uD m θ m ](y) + K s ,m (r)[θ m (y) -θ m (x)]u(y) (25) f (α) = J(η) K s (a) a i u(x) + [ u +uD m θ m ](x) ( 26
)
and where all p-dependent quantities, like r, a, are evaluated at p = 0 (note that ∆, i.e the function ρ(α), does not depend on p). Moreover, the derivatives a and a are given by

a (α) = ∂ ∂p {a 1 cos α + a 2 sin α} = ∂ 2 y ∂ξ 1 ∂p ξ=η cos α + ∂ 2 y ∂ξ 2 ∂p ξ=η sin α = ∂ ∂ξ 1 ∂y ∂p ξ=η cos α + ∂ ∂ξ 2 ∂y ∂p ξ=η sin α = ∂θ ∂ξ 1 | ξ=η cos α + ∂θ ∂ξ 2 | ξ=η sin α = ∇θ(x).a(α) (27) a (α) = a(α).∇θ(x).a(α) a(α) (28)
Using the latter formulas and the expansion

θ(y) -θ(x) = ρ∇θ(x).a + O(ρ 2 ) = ρ a (α) + O(ρ 2 ) ( 29
)
about ρ = 0, it is easy to show that f /ρ is indeed the singular part of F , and also that K s has the same degree of singularity than K s (y -x). As a result, { F -f /ρ} is regular. From ( 25),

(26), it is then readily seen that the material derivative I S as given by ( 24) has the form

I S = - S {[K s u] + K s uD m θ m } dS + R (30)
where the residual term, given by

R = - 2π 0 f (α) a (α) a(α) dα ( 31
)
remains to be examinated. Indeed, from ( 28), one has

f (α) a (α) a(α) = f (α) a(α).∇θ(x).a(α) a 2 (α)
and then, from property (14b) and the definition (18) of A(α):

f (α + π) a (α + π) a(α + π) = -f (α) a (α) a(α) (32)
As a result, the residual integral R (31) vanishes:

R = - 2π 0 f (α) a (α) a(α) dα = - π 0 + 2π π f (α) a (α) a(α) dα = 0 (33)
Thus we have proved that the material differentiation formula ( 6) is generalizable to CPV singular integrals with the kernel (13):

d dp - Sp K s (y p -x p )u(y p , p) dS y = - S {[K s (y -x)u(y)] + K s (y -x)u(y)D m θ m (y)} dS y (34)
Comments. The result (34) obviously extends for any linear combination of kernels (13); also, the unit normal n, if present, can be incorporated without difficulty in the regular factor u and is therefore not excluded from the analysis. All usual potential or isotropic elastostatic strongly singular BIE formulations are encompassed as a direct consequence of the result (34).

The analysis presented here assume that S p is smooth at x p . The kernel ( 13) can be generalized to any combination of kernels of the form |r| -2 (e r ⊗ . . . ⊗ e r ) where the unit vector e r = r/ |r| appears an odd number of times in the tensor product, since properties ( 14) still hold.

The approach of [START_REF] Barone | A boundary element approach for recovery of shape sensitivities in threedimensional elastic solids[END_REF] and others thus receives an a posteriori validation, at least when the free term present in singular BIE formulations remains constant under the geometrical transformation (1). This is normally the case when dealing with initially smooth singular points x. Although the present analysis does not answer this issue, it is conjectured that the present result remains true for a corner or edge point x, for which the free-term is likely to vary under the geometrical transformation. The fact that the material derivative of a regularized displacement BIE formulation is valid and has the same mathematical form for smooth and corner singular points x ( [START_REF] Bonnet | Regularized BIE formulations for first-and second-order shape sensitivity of elastic fields[END_REF]) provides in our view an indirect proof of this conjecture.

Finally, the fundamental solutions for dynamic problems (wave propagation, elastodynamics, parabolic heat equation) do not in general verify the homogeneity properties (9a), (14a), (37a). However, their singularity is identical to that of its static counterpart. As a result, upon splitting a dynamic fundamental into the sum of its (singular) static counterpart and a (nonsingular) complement, it is apparent that the result (34) validates the material differentiation of dynamic strongly singular BIE formulations as well.

Hypersingular integral on a changing surface

Let us now consider the case of hypersingular kernels. Such kernel functions are usually derivatives with respect to source point coordinates of strongly singular kernels:

K h (y -x) = ∂ ∂x K s (y -x) = -K s , (y -x) (35)
In particular, K h (y -x) is assumed to be divergence-free:

∂ ∂x K h (y -x) = -K h , (y -x) = 0 (y = x)
For three-dimensional potential theory or elastostatics, K h (y -x) usually appears as a linear combination of

K h (z) = K s , = 4z i z j z k z |z| 7 - δi z j z k |z| 5 - δj z i z k |z| 5 - δk z i z j |z| 5 (36)
where 1 ≤ i, j, k, ≤ 3 are fixed indexes; besides, they are homogeneous of degree -3 and symmetric:

K h (az) = 1 a 3 K h (z) K h (-z) = K h (z) (37)
If need be, the unit normal components are included in the nonsingular factor u, in order to be still able to invoke the homogeneity property (37a). Following [START_REF] Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF], hypersingular boundary integral formulations involve terms defined by a limiting process ('finite part') of the form

I H (p) = lim ε→0 u(x p , p) sε(x p ) K h (y p -x p ) dS y + Sp-eε(x p ) K h (y p -x p )u(y p , p) dS y ≡ = S K h (y -x)u(y) dS y (38)
(where again a spherical exclusion neighbourhood v ε (x p ) of vanishing radius ε centered at x p is removed around x p ). We follow again the line of investigation used in sections 4, 5. From [START_REF] Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF], the limit ( 38) is known to exist and a rigorous and computable expression of I H using the parametrization ( 16) is available for any p > 0. We now recall the latter result, using again the notations introduced in section 4, together with some additional notation related to the second-order expansion of the position vector r:

r p = ρr(ρ, α, p) = ρ[a(α, p) + ρb(α, p) + o(ρ)] (39
The kernel function K h (r p ) then takes the form

K h (r p ) = 1 ρ 3 K h (r) = 1 ρ 3 K h (a) + ρK h ,m b m (α, p) + O(ρ 2 ) (42)
The integral ( 38) is given by:

I H (p) = 2π 0 ρ(α) 0 F (ρ, α, p) - f (α, p) ρ 2 - g (α, p) ρ dρdα + 2π 0 g (α, p) ln[ρ(α)a(α, p)] + f (α, p) (a.b) a 2 - 1 ρ (α, p) dα (43)
where

F (ρ, α, p) = 1 ρ 2 K h (r)u(y p , p)J(ξ, p) (44) f (α, p) = K h (a)u(x p , p)J(η, p) (45) g (α, p) = K h ,m (a)b m (α, p)u(x p , p)J(η, p) + K h (a) ∂ ∂ρ [uJ] ρ=0 (46)
and a and r are still defined by eqs. ( 17) and ( 18). In fact, f (α, p) and g (α, p) (denoted F -2 and F -1 in [START_REF] Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF]) are obtained as the coefficients of ρ 0 = 1 and ρ in the Taylor expansion of K h (r)u(y p , p)J(ξ, p) with respect to ρ about ρ = 0. The symmetry property (37b) of K h , applied to eqs. ( 45) and ( 46), imply that:

f (α + π, p) = f (α, p) g (α + π, p) = -g (α, p) (47)
7 Material derivative of the hypersingular integral

The formula (43) holds for any fixed p ≥ 0, the parameter space ∆ being independent of p; moreover the double integral is regular because of the singularity removal achieved through the introduction of f (α, p) and g (α, p). Thus eq. ( 43) can be differentiated under the integral sign at p = 0:

d dp I H = 2π 0 ρ(α) 0 F (ρ, α) - 1 ρ 2 f (α) - 1 ρ g (α) dρdα + 2π 0 g (α) ln[ρ(α)a(α)] -g (α) a (α) a(α) dα + 2π 0 f (α) (a.b) a 2 - 1 ρ (α) + f (α) (a.b) a 2 - 1 ρ (α)dα (48)
with

F (ρ, α) = 1 ρ 2 J(ξ) K h ,m (r)[θ m (y) -θ m (x)]u(y) + K h (r)[ u +uD m θ m ](y) (49) f (α) = J(η) K h ,m (a) a m u(x) + K h (a)[ u +uD m θ m ](x) (50) g (α) = K h ,mn (a) a n (α)b m (α)u(x)J(η) + K h ,m (a) b m (α)u(x) + b m (α)[ u +uD m θ m ](x) J(η) + K h ,m (a) a m (α) ∂ ∂ρ [uJ] ρ=0 + K h (a) ∂ ∂ρ [ u +uD m θ m ]J ρ=0 (51)
where a and r are still defined by eqs. ( 17) and ( 18). We note that a is given by ( 27) and that, similarly

b (α) = ∂ ∂p b 11 cos 2 α + 2b 12 cos α sin α + b 22 sin α = ∂ 2 θ ∂ξ 1 ∂ξ 1 | ξ=η cos 2 α + 2 ∂ 2 θ ∂ξ 1 ∂ξ 2 | ξ=η cos α sin α + ∂ 2 θ ∂ξ 2 ∂ξ 2 | ξ=η sin 2 α (52)
Moreover, differentiation of eq. ( 40) with respect to p gives, accounting for the fact that the variables p and ξ β are independent, the following expansion:

θ(y) -θ(x) = ρ[ a (α) + ρ b (α) + o(ρ)] (53)
Using the previous formulas, it is then easy to show that

F (ρ, α) - 1 ρ 2 f (α) - 1 ρ g (α) = O(1)
is indeed nonsingular. Moreover, using the same argument as in section 5, the antisymmetry property (47b) implies that Then, using eqs. ( 19), ( 41), which in turn imply that

a (α + π) = -a (α) b (α + π) = b (α)
it is easy to conclude that a.b a 2 -

1 ρ (α + π) = a.b a 2 - 1 ρ (α)
and consequently, using (47a), that Comments. The result (54) obviously extends for any linear combination of kernels (36); also, the unit normal n, if present, can be incorporated without difficulty in the regular factor u and is therefore not excluded from the analysis. All usual potential or isotropic elastostatic hypersingular BIE formulations are encompassed as a direct consequence of the result (34). All other comments made at the end of section 5 still apply when transposed in the obvious way to the result (34) for hypersingular boundary integrals.

Zhang, Q.; Mukherjee, S. 1991: Second-order design sensitivity analysis for linear elastic problems by the derivative boundary element method., Comp. Meth. in Appl. Mech. Engng.,

  ) where b(α, p) = b 11 (p) cos 2 α + 2b 12 (p) cos α sin α + b 22 (p) sin 2 α (40) and b αβ (p) = ∂ 2 y p ∂ξ α ∂ξ β | ξ=η Note that formula (40) implies that: b(α + π, p) = b(α, p)

  , we have proved that the material differentiation formula (6) is generalizable to hypersingular boundary integrals with kernels (36):d dp = Sp K h (y px p )u(y p , p) dS y = = S [K h (y -x)u(y)] + K h (y -x)u(y)div S θ(y) dS y (54)

Acknowledgements

This manuscript has benefited from various discussions with Prof. M. Guiggiani of Siena University, Italy.


Marc Bonnet
email: bonnet@athena.polytechnique.fr

Massimo Guiggiani
email: guiggiani@unisi.it

Prof M Guiggiani

Dip Di Costruzioni

Tangential derivative of singular boundary integrals with respect to the position of collocation points 1

Keywords: sensitivity analysis, singular boundary integral equations, hypersingular integrals, boundary element method

This paper investigates the evaluation of the sensitivity, with respect to tangential perturbations of the singular point, of boundary integrals having either weak or strong singularity. Both scalar potential and elastic problems are considered. A proper definition of the derivative of a strongly singular integral with respect to singular point perturbations should accommodate the concomitant perturbation of the vanishing exclusion neighbourhood involved in the limiting process used in the definition of the integral itself. This is done here by resorting to a shape sensitivity approach, considering a particular class of infinitesimal domain perturbations that "move" individual points, and especially the singular point, but leave the initial domain globally unchanged. This somewhat indirect strategy provides a proper mathematical setting for the analysis. Moreover, the resulting sensitivity expressions apply to arbitrary potential-type integrals with densities only subjected to some regularity requirements at the singular point, and thus are applicable to approximate as well as exact BEM solutions. Quite remarkable is the fact that the analysis is applicable when the singular point is located on an edge and simply continuous elements are used. The hypersingular BIE residual function is found to be equal to the derivative of the strongly singular BIE residual when the same values of the boundary variables are substituted in both SBIE and HBIE formulations, with interesting consequences for some error indicator computation strategies.

Introduction

Boundary integral equation (BIE) formulations lead to numerical techniques, namely boundary element methods (BEMs), that are now well established and widely applied. After a boundary value problem has been numerically solved using a standard boundary element method, it is sometimes necessary, or merely worthwhile, to be able to formulate and compute directional derivatives of field variables, or more generally of BIE residuals, at an observation point on the boundary and along a tangential direction. This of course leads to investigation of such derivatives for singular integrals of scalar or elastic potential theory. A relevant practical issue is the evaluation of in-plane strain, hoop stresses, etc. from a given (usually approximate) BEM solution. Another one concerns the evaluation of derivatives of BIE residuals under perturbations of the collocation points; a related idea is used to define local a posteriori error indicators in [START_REF] Guiggiani | Sensitivity analysis for boundary element error estimation and mesh refinement[END_REF].

Concerning derivatives of the most singular integral operators encountered in standard BIE formulations (i.e. of the form of double-layer potentials), the Lyapunov-Tauber theorems of classical potential theory give results for the limit on the boundary of the normal derivative of scalar doublelayer potentials or of the traction operator applied to vector elastic potentials, see e.g. the books by Kellogg [11] and Kupradze et al. [START_REF]Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity[END_REF]. These results allow one to formulate Neumann problems using either single-or double-layer potentials. On the other hand, general results for tangent derivatives of strongly singular potential integrals are not given in the classical literature on potential theory.

More recently, many studies have been devoted to so-called hypersingular boundary integral equation (HBIE) formulations, which basically result from differentiation of strongly singular (SBIE) formulations with respect to components of the collocation point. Of course, both SBIE and HBIE give an identically zero residual (considered as a function of the collocation point) when the exact solution of a given problem is substituted into them. In practice, BIEs are discretized and collocated at a finite number of boundary points, and an approximate solution is obtained. For instance, the discrete SBIE is solved, and the SBIE residual may take nonzero values at boundary points that were not used for collocation. Then, if the previous approximate solution is substituted into the HBIE, is the HBIE residual function equal to the derivative of the SBIE function? This question is related to the study of e.g. local error estimation [START_REF] Guiggiani | Sensitivity analysis for boundary element error estimation and mesh refinement[END_REF]. Its answer is not a priori obvious because both BIE formulations are defined in terms of a limiting process involving a vanishing exclusion neighbourhood of the singular point y, so that in a proper definition of the derivative of a strongly singular integral with respect to y the vanishing neighbourhood should "move" with y.

We address this particular difficulty in a somewhat indirect way, by resorting to a domain sensitivity approach. Specifically, the main idea used in this paper is the introduction of domain perturbations, depending on a time-like parameter p, that "move" individual points but leave the initial domain (p = 0) globally unchanged. This way, the concomitant change of singular point and exclusion neighbourhood is properly taken into account. The analysis then uses adequately chosen domain differentiation formulas. The sensitivity with respect to tangential perturbations of the singular point of boundary integrals having either weak or strong singularity, of the type used for both scalar and elastic BIE formulations, is formulated. From the results obtained, it is apparent that the HBIE residual function is indeed equal to the derivative of the SBIE residual function when the same values of the boundary variables are substituted in both SBIE and HBIE formulations. It is worth pointing out that this relationship between residuals issue is not raised in the formulation of hypersingular representations of potential gradients or stress tensors on the boundary, like in [START_REF] Guiggiani | Hypersingular formulation for boundary stress evaluation[END_REF], although they allow one to evaluate tangential derivatives of the potential or displacement field. In particular the latter approach uses the fact that the density variables solve the relevant governing equations, whereas the present analysis addresses arbitrary potential-type integrals with densities only subjected to some (fairly weak) regularity requirements at the singular point. An advantage of this viewpoint is that the resulting formulas apply to approximate as well as exact BEM solutions, which is convenient in devising strategies for error indicator computation.

Basic definitions and relations

In this section we summarize some basic results that will be used extensively throughout the paper.

More details can be found in previous papers dealing with shape sensitivity of BIE formulations [START_REF] Bonnet | Regularized BIE formulations for first-and second-order shape sensitivity of elastic fields[END_REF], [START_REF] Bonnet | Differentiability of strongly singular and hypersingular boundary integral formulations with respect to boundary perturbations[END_REF], [START_REF] Zhang | Second-order design sensitivity analysis for linear elastic problems by the derivative boundary element method[END_REF] and in Appendix A of the recent book [START_REF] Bonnet | Equations intégrales et éléments de frontière[END_REF].

Let us consider, in a three-dimensional Euclidean space equipped with a Cartesian orthonormal basis (e 1 , e 2 , e 3 ), a body Ω p whose shape depends on a time-like parameter p. The body Ω = Ω 0 corresponding to p = 0 is traditionally assumed as the "initial" configuration. Points of Ω p are denoted by x p and y p , whereas points of Ω are simply denoted by x and y.

Let the initial domain Ω be bounded by a Kellogg's regular surface ∂Ω = Γ with outer unit normal n [START_REF] Kellogg | Fundations of potential theory[END_REF]. On the boundary we can define some tangential differential operators. First, we define the tangential gradient ∇ S f of a scalar field f defined on Ω

∇ S f (x) := ∇f -n (∇f • n) = ∇f -n ∂f ∂n (1)
which is the projection of the gradient ∇f onto the tangent plane at x ∈ ∂Ω.

The r-th Cartesian component of ∇ S f can be denoted by D r f if we adopt the following definition for the tangential partial derivative D r

D r f := f, r -n r ∂f ∂n = f, r -n r (f, m n m ), (2)
where the notation f, m = ∂f /∂x m has been introduced.

The surface divergence div S v of a vector field v can now be defined as

div S v := ∇ S • v = D m v m . (3)
It is convenient to define another tangential differential operator D rs

D rs f := n r D s f -n s D r f = n r f, s -n s f, r . (4)
Among the many versions of Stokes' theorem, we are mainly interested in the following ones

S D rs f dS = e jrs ∂S f τ j ds, (5)
S D r g dS = ∂S gν r ds - S gn r 2κ dS ( 6
)
where ∂S is the contour of the (possibly) open surface S and

κ(x) = - 1 2 div S n(x) = - 1 2 D a n a (x)
is the mean curvature of S at x. The vector τ is defined by

τ = n × ν, (7)
where ν is the unit vector normal to ∂S and tangent to S and pointing towards the outside of S. e jrs denote the permutation symbols (e 123 = e 231 = e 312 = 1, e 132 = e 213 = e 321 = -1, and equal to zero otherwise). Formula ( 6) is obtained by putting f = gn s in [START_REF] Guiggiani | Error indicators for adaptive mesh refinement in the Boundary Element Method-A new approach[END_REF].

A geometrical transformation from Ω onto Ω p is given by a mapping

x ∈ Ω → x p = Φ(x; p) ∈ Ω p , (8)
where ∀x ∈ Ω, Φ(x; 0) = x, as already stated. Similarly to continuum mechanics, it is assumed that any such mapping is a diffeomorphism between Ω and Ω p .

The corresponding initial transformation velocity θ(x) is a vector field on Ω defined by

θ(x) := ∂Φ ∂p p=0 . (9)
In other words, the geometrical transformation (8) can be expanded in the following form

x p = x + θ(x) p + o(p).
Now, let us consider a generic scalar field f (p, x p ) defined on the changing domain Ω p and also depending explicitly on p. Its total derivative * f with respect to p, evaluated at p = 0, is given by the following well-known expression

* f (x) := df dp p=0 = ∂f ∂p + f, m θ m (x) p=0 = ḟ + ∇f • θ(x) p=0 , (10)
where ḟ = ∂f /∂p. Depending on the context, the total derivative * f can also be termed material or Lagrangian derivative. The function ḟ can be consistently called Eulerian derivative.

The kernel functions that appear in boundary integral equations are two-point functions like K(x py p ) which do not depend explicitly on p. Adopting the notation just defined, and assuming that both points follow the geometrical transformation (8), we have from [START_REF] Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF] the following expression for the total derivative of K with respect to p * K(y, x)

:= dK dp p=0 = K, m (x -y)[θ m (x) -θ m (y)], (11)
where K, m = ∂K/∂x m = -∂K/∂y m .

The total derivative of the integral of a nonsingular function f (p, x p ) over a changing (possibly open) surface S p is given, from e.g. [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF], by either of the formulas

d dp Sp f (p, x p ) dS p=0 = S * f + f div S θ dS (12) = S ḟ + ∂f ∂n -2κf (θ • n) dS + ∂S (θ • ν)f ds, (13)
where S = S 0 . Formulas ( 12) and ( 13)) are equivalent; they relate the derivative of the integral to Lagrangian or Eulerian derivatives of the function f , respectively.

According to Petryk and Mroz [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF], we can express the total (material) derivative of each Cartesian component of the unit normal n(x p ) as

* n i (x) := dn i dp p=0 = -n m (D i θ m ). ( 14
)

Singular boundary integrals

Singular boundary integrals arise naturally in the boundary element method (BEM) whenever the source (collocation) point y is taken to the boundary. Due to the singularity of the kernel functions K involved, a limiting process is necessary for the definition of these integrals.

Let y be a fixed point on the boundary ∂Ω of a three-dimensional domain Ω. We consider an exclusion neighbourhood v ε (y) of y, of radius ≤ ε (Figure 1). For any ε > 0, y is always an external point for the domain Ω ε (y) = Ω -v ε (y) whose boundary ∂Ω ε is given by

∂Ω ε = (∂Ω -e ε ) + s ε = Γ ε + s ε ,
where

s ε = Ω ∩ ∂v ε , e ε = ∂Ω ∩ vε , and Γ ε = ∂Ω -e ε .
The typical singular boundary integral can be defined by means of a limiting process

I(y) = lim ε→0 Γε(y)+sε(y) K(x -y)v(x) dS x , (15)
where the density v(x) is a sufficiently regular function.

There are instances (e.g., [START_REF] Guiggiani | Sensitivity analysis for boundary element error estimation and mesh refinement[END_REF]) in which we may assume that the position of the source point y on ∂Ω depends on a time-like parameter p. Since the vanishing neighbourhood v ε (y) is attached to y, the parameter p also affects the surface of integration in (15) (although p does not affect ∂Ω). Somewhat similar considerations occur for the correct evaluation of the Cartesian derivatives of the integral representation for displacement at interior points y in the presence of initial strain [START_REF] Bui | Some remarks about the formulation of three-dimensional thermoelastoplastic problems by integral equations[END_REF]. Obviously, a necessary condition for the function y p = y(p) to be a diffeomorphism is that the point y lies either on a smooth part of ∂Ω or, at most, on an edge (two tangent planes). In the last case, the point y p has to move along the edge itself.

Still following ref. [START_REF] Guiggiani | Sensitivity analysis for boundary element error estimation and mesh refinement[END_REF] and for greater generality, the density function v can be taken to explicitly depend on p (this is the case, e.g., of approximate numerical solutions that clearly depend on the specific set of collocation points).

Hence, we can define the function I(p)

I(p) = lim ε→0 Γε(y p )+sε(y p ) K(x -y p )v(p, x) dS x , (16)
which can be seen as I(y p ) with a p-dependent density function.

In this paper we address the problem of the evaluation of the derivative

dI dp p=0 . ( 17
)
Derivatives of this type are basically what are required for the evaluation of the sensitivity of an approximate BEM solution with respect to the position of the collocation points [START_REF] Guiggiani | Sensitivity analysis for boundary element error estimation and mesh refinement[END_REF].

Tangential derivatives of singular boundary integrals

To evaluate the derivative dI/dp we rely on a domain sensitivity approach. Although in the present context neither the shape of the body nor the discretization are affected by the parameter p, the sensitivity viewpoint provides a convenient mathematical means to describe exclusion neighbourhoods v ε (y) moving according to perturbations of the singular point y.

To this end, let us consider a geometrical transformation x → Φ(x, p) as in ( 8) but with initial transformation velocity θ (defined in ( 9)) such that

θ(x) • n(x) = 0 (∀x ∈ ∂Ω),
and θ(y) = 0 (18)

In other words, the transformation velocity is almost everywhere orthogonal to the boundary normal of the body. Under any such transformation, Ω is globally invariant (up to order O(p)), while

y → y + a p + o(p),
where

a = θ(y) (19)
is the transformation velocity of the source point. Therefore, the initial (i.e. of first-order with respect to p about p = 0) effect of the perturbation of the position of the singular point y on I(p) is adequately taken into account. At the same time, v ε (y) is "moved" (and possibly distorted, but this is irrelevant)

along with the singular point, so that the singular integral ( 16) is properly defined for small p ≥ 0.

Notice that if ∂Ω has more than two tangent planes at y, then there are no transformation velocity fields such that a • n(y) = 0.

Essentially, with a transformation satisfying (18) we can use all the results already available from sensitivity analysis, yet we are only considering the effect of the perturbation of the singular point y and of the attached neighbourhood v ε . Formulas for the derivative of an integral over a changing surface were given in [START_REF]Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity[END_REF] and [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF]. The derivative of (16) can be conveniently evaluated applying formula [START_REF]Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity[END_REF] to the integral on s ε and formula [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF] to the integral on Γ ε (all derivatives are understood at p = 0)

dI dp = lim ε→0 sε [(Kv) * + Kv D m θ m ] dS + Γε (Kv) q dS + ∂Γε (θ m ν m )Kv ds . (20)
In this way, full advantage is taken from the orthogonality condition (18).

Although quite obvious, we may also observe that

(Kv) * = * Kv + K * v = K, m [θ m (x) -a m ]v + K * v, (21) and
(Kv)

q = Kv + K v = ∂K ∂y m a m v + K v = -K, m a m v + K v. ( 22
)
In this paper we consider the tangential derivative of singular integrals with respect to the source point y for three-dimensional potential and elastic problems. It is worth recalling that the analysis is not restricted to smooth boundaries, but also covers the case of y on an edge. In Appendix A it is shown that the final results are not affected by the shape of v ε . Therefore, it is always assumed to employ a symmetric neighbourhood v ε (y) of uniform radius ε.

Weakly singular kernels

The formulation of cartesian derivatives of single-layer scalar or elastic potentials has been studied

in classical potential theory. The results of this section are in effect nothing more than extensions to p-dependent density functions of classical results on cartesian derivatives of single-layer scalar or elastic potentials, see e.g. [START_REF] Kellogg | Fundations of potential theory[END_REF] for scalar potentials and [START_REF]Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity[END_REF] for elasticity. However, it is interesting to see the present sensitivity approach at work in this context before dealing with the more difficult case of double-layer potentials.

Potential problems

Let us consider the following boundary integral with weakly singular kernel

I(p) = lim ε→0 Γε(y p )+sε(y p ) G(x -y p ) q(p, x) dS x , (23)
where the kernel K = G = 1/(4πr), with r = |x -y|, is the fundamental solution for 3D potential problems, and the density function q is a sufficiently regular function (usually an approximation to the normal derivative of the potential).

A first expression for dI/dp is given by (20), where (Gq) q can be expanded according to ( 22)

dI dp = lim ε→0 sε [(Gq) * + Gq D m θ m ] dS + Γε G q + a m ∂G ∂y m q dS + ∂Γε (θ m ν m )Gq ds . ( 24
)
For 3D problems we have on

s ε that G = O(ε -1 ), G, m = O(ε -2 ), and θ m (x) -a m = O(ε α ), with
α > 0, because of the regularity assumption on the geometrical transformation Φ. If we also assume that | * q(x) - * q(y)| ≤ Cr α-1 (which is quite reasonable), we can conclude that lim ε→0 sε

[(Gq) * + Gq D m θ m ] dS = 0. ( 25
)
Moreover, when y is on an edge of tangent a, the limit of the line integral on ∂Γ ε is zero

lim ε→0 ∂Γε (θ m -a m )ν m Gq ds + a m ∂Γε ν m Gq ds = 0 ( 26
)
since it is zero on each side (provided q is piecewise C 0,α at y). Essentially, owing to the scalar product a m ν m , the line integral on each side amounts in the limit to the integral from 0 to π of cos ϕ dϕ.

Therefore, we have that

dI dp p=0 = lim ε→0 Γε G q + a m ∂G ∂y m q dS = ∂Ω G q dS + a m - ∂Ω W m q dS = ∂Ω G q dS + - ∂Ω Ġq dS (27)
where W m = ∂G/∂y m are the kernel functions of the hypersingular boundary integral equations [START_REF] Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF] and Ġ = a m W m . The weak singularity of G allows the first integral in (25) to be defined as an ordinary integral on ∂Ω. On the other hand, the second integral with strongly singular kernel Ġ is equivalent (owing to the assumed circular shape of ∂Γ ε ) to a Cauchy principal value [START_REF] Guiggiani | Error indicators for adaptive mesh refinement in the Boundary Element Method-A new approach[END_REF].

Elastic problems

In this case we have to consider the tangential derivative of the weakly singular integral

I(p) = lim ε→0 Γε(y p )+sε(y p ) U ij (x -y p )t j (p, x) dS x , (28)
where the kernels U ij , are the fundamental solutions for 3D elastic problems, and the density functions t j are just sufficiently regular functions. Again, a possible expression for dI/dp is given by (20).

Using the same arguments as above, we can show that lim ε→0 sε

[(U ij t j ) * + U ij t j D m θ m ] dS = 0. (29)
Similarly, when y is on an edge of tangent a the line integral on ∂Γ ε is zero since it is zero on each side (provided, on each element, t j ∈ C 0,α at y).

Therefore, we have that

dI dp p=0 = lim ε→0 Γε U ij ṫj + a m ∂U ij ∂y m t j dS = ∂Ω U ij ṫj dS + a m - ∂Ω W ijm t j dS = ∂Ω U ij ṫj dS + - ∂Ω Uij t j dS (30)
where

W ijm = ∂U ij /∂y m and Uij = a m W ijm .
Formulas ( 27) and (30) have the same structure. They show that the derivative dI/dp is given by a term that accounts for the dependence on p of the density function and another term that takes into account the effect of the changing singular point on the kernel function. This last term is the one predicted by classical potential theory.

6 Strongly singular kernels

Potential problems

Let us consider the following boundary integral with strongly singular kernel T = G, j n j

I(p) = lim ε→0 Γε(y p )+sε(y p ) T (y p , x)u(p, x) dS x = Γε(y p )+sε(y p ) G, j (x -y p )n j (x)u(p, x) dS x , (31)
where, in general, the density function u is only an approximation to a potential.

First we consider the derivative of the integral on s ε , for fixed ε > 0. According to (20), we apply formula ( 12) and obtain

d dp sε(y p ) T u dS p=0 = sε(y) [(T u) * + T uD i θ i ]dS, (32)
This expression can be given a more explicit form using [START_REF] Kellogg | Fundations of potential theory[END_REF] and ( 14)

(T u) * + T uD i θ i = (G, j n j u) * + G, j n j uD i θ i (33) = G, j n j * u + * G, j n j u + G, j * n j u + G, j n j uD i θ i (34) = G, j n j * u + [G, ij n j (θ i -a i ) + G, j (n j D i θ i -n i D j θ i )] u (35) = G, j n j * u + [(D ji G, j )(θ i -a i ) + G, j D ji θ i ] u (36) = G, j n j * u + D ji (G, j [θ i -a i ]) u ( 37
)
where a i = θ i (y). Expression (36) comes from (4) once with f = θ i , and once with f = G, j with the additional fact that G, ii = 0. The final step is possible since D ji a i = 0 (a i do not depend on x).

We can now add and subtract the same terms to obtain

(T u) * + T uD i θ i = G, j n j [ * u(x) - * u(y)] + D ji (G, j [θ i -a i ])[u(x) -u(y)] + * u(y)G, j n j + u(y)D ji (G, j [θ i -a i ]). ( 38
)
At this stage we have to make some assumptions on the regularity of the density function u(x) at y. We assume u ∈ C 0,α (y) and * u ∈ C 0,α (y), along with θ ∈ C 0,1 (y) which comes from the geometrical transformation Φ being a diffeomorphism. It is worth noting that we do not require u ∈ C 1,α (y) (i.e., ∇u ∈ C 0,α (y)), but only the total derivative * u(y) = u(y) + ∇u • a (39) to be Hölder continuous at y. This is a condition that is satisfied by standard continuous boundary elements provided a is tangent to the interelement border between two elements.

Under these assumptions, the integrals on s ε of the first two terms on the r.h.s. of eq. (38) vanish in the limit, as ε → 0.

Therefore, we have that

sε(y) [(T u) * + T uD i θ i ]dS = * u(y) sε(y) G, j n j dS + u(y) sε(y) D ji (G, j [θ i -a i ]) dS + o(1). ( 40
)
The limit of the first integral in the above expression is readily obtained

lim ε→0 * u(y) sε(y) G, j n j dS = * u(y) lim ε→0 sε(y) T dS = * u(y)c(y) = u(y) + ∂u ∂a c(y), (41)
where, by definition, c(y) is the free-term coefficient of the usual boundary integral equation.

We already found that D ji G, j = G, ij n j . Hence from the final term in (40) we get -u(y)a i sε(y)

D ji G, j dS = -u(y)a i sε(y) G, ij n j dS = u(y)a i sε(y) V i dS = u(y)a i b i ε + d i + O(ε) , (42)
that is the same free-term coefficients of the hypersingular boundary integral equations as obtained in [START_REF] Guiggiani | Hypersingular boundary integral equations have an additional free term[END_REF]. In fact the hypersingular kernel is

V i = -G, ij n j .
Stokes' theorem in the form (5) can be applied to the function D ji (G, j θ i ) in ( 40)

u(y) sε(y) D ji (G, j θ i )dS = u(y)e jip ∂sε G, j θ i τp ds = -u(y) ∂Γε T θ i ν i ds, (43)
where we used the fact that ∂s ε and ∂Γ ε are identical curves with opposite tangents τ = -τ . The last step in (43) comes from the following general result

e jip w i τp = -(w × τ ) j = -(w • ν)n j + (w • n)ν j ,
where w is a general vector, τ = -τ , and (n, ν, τ ) is a Cartesian orthonormal basis defined in [START_REF] Guiggiani | Hypersingular boundary integral equations have an additional free term[END_REF]. If the vector w is replaced by θ we have that

e jip θ i τp = -(θ × τ ) j = -(θ • ν)n j ,
because of the assumed orthogonality between θ and n (eq. ( 18)).

Collecting the results of (41), (42), and (43), we can write the derivative of the integral on s ε as

d dp sε(y p ) T u dS p=0 = sε(y) [(T u) * + T uD i θ i ]dS, = u(y) + ∂u ∂a c(y) + u(y)a i b i ε + d i + O(ε) -u(y) ∂Γε T θ i ν i ds. ( 44
)
We can now turn our attention to the integral on Γ ε in (31). To compute the derivative with respect to p, for fixed ε, we can employ the general relation [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF] with the obvious simplification due

to θ • n = 0 on Γ ε d dp Γε T udS p=0 = Γε T udS + Γε Ṫ udS + ∂Γε T uθ i ν i ds = Γε T udS -a i Γε (G, ji n j )udS + u(y) ∂Γε T θ i ν i ds + ∂Γε T [u -u(y)]θ i ν i ds = Γε T udS + a i Γε V i udS + u(y) ∂Γε T θ i ν i ds + o(1), (45)
where Ṫ = a i V i .

When contributions from s ε and Γ ε are taken together, we observe that the line integrals in (44) and (45) cancel each other. Therefore, we are left with the following final expression

dI dp p=0 = lim ε→0 d dp Γε(y p )+sε(y p ) T udS p=0 = lim ε→0 * u(y) sε T dS + Γε T udS + u(y) sε Ṫ dS + Γε Ṫ udS + a i ∂Γε T [u -u(y)]ν i ds = lim ε→0 u(y) + ∂u ∂a sε T dS + Γε T udS + u(y)a i sε V i dS + a i Γε V i udS = u(y)c(y) + ∂u ∂a c(y) + - ∂Ω T u dS + lim ε→0 a i Γε V i u dS + u(y) a i b i ε + d i . ( 46
)
We recall that the above expression has been obtained by assuming u and only its directional derivative ∂u/∂a to be C 0,α at the point y lying on an edge of tangent a. This is typically satisfied between two standard continuous boundary elements.

Elastic problems

Let us consider the following boundary integral with strongly singular kernel T ij = Σ ijk n k , where Σ ijk is the stress tensor associated to the displacement field U ij (fixed index i)

I(p) = lim ε→0 Γε(y p )+sε(y p ) T ij (y p , x)u j (p, x) dS x = lim ε→0 Γε(y p )+sε(y p ) Σ ijk (x -y p )n k (x)u j (p, x) dS x , (47)
where, in general, the density functions u j are only an approximation to displacement field components.

As for the potential case, we consider first the derivative of the integral on s ε , for fixed ε > 0.

According to (20), we apply formula ( 12) and obtain

d dp sε(y p ) T ij u j dS p=0 = sε(y) [(T ij u j ) * + T ij u j D θ ]dS, ( 48
)
This expression can be given a more explicit form using ( 11) and ( 14)

(T ij u j ) * + T ij u j D θ = (Σ ijk n k u j ) * + Σ ijk n k u j D θ (49) = Σ ijk n k * u j + * Σ ijk n k u j + Σ ijk * n k u j + Σ ijk n k u j D θ (50) = Σ ijk n k * u j + [Σ ijk, n k (θ -a ) + Σ ijk (n k D θ -n D k θ )]u j (51) = Σ ijk n k * u j + [(D k Σ ijk )(θ -a ) + Σ ijk D k θ ]u j (52) = Σ ijk n k * u j (x) + D k (Σ ijk [θ -a ])u j (x) (53)
Expression (52) comes from (4) once with f = θ , and once with f = Σ ijk with the additional fact that Σ ijk,k = 0.

Following the development of (38) in the potential case, we can now add and subtract the same terms in (53) which, along with the assumptions u j ∈ C 0,α (y) and * u j ∈ C 0,α (y), allows us to write (cfr. eq. ( 40))

sε(y) [(T ij u j ) * +T ij u j D θ ]dS = * u j (y) sε(y) Σ ijk n k dS +u j (y) sε(y) D k (Σ ijk [θ -a ]) dS +o(1). ( 54
)
The limit of the first integral on the r.h.s. of the above expression is readily obtained

lim ε→0 * u j (y) sε(y) Σ ijk n k dS = * u j (y) lim ε→0 sε(y) T ij dS = * u j (y)c ij (y) = uj (y) + ∂u j ∂a c ij (y), (55)
where, by definition, c ij (y) are the free-term coefficients of the usual boundary integral equation.

The last term in (54) can be converted in a more convenient form if we observe that

D k Σ ijk = Σ ijk, n k -u j (y)a sε(y) D k Σ ijk dS = -u j (y)a sε(y) Σ ijk, n k dS = u j (y)a sε(y) V ij dS = u j (y)a b ij ε + d ij + O(ε) , (56)
that is the same free-term coefficients of the hypersingular boundary integral equations as obtained in [START_REF] Guiggiani | Hypersingular boundary integral equations have an additional free term[END_REF]. In fact the hypersingular kernel is

V ij = -Σ ijk, n k .
Finally, Stokes' theorem in the form ( 5) can be applied to the function D k (Σ ijk θ ) in (54)

u j (y) sε(y) D k (Σ ijk θ )dS = u j (y)e k p ∂sε Σ ijk θ τ ds = -u j (y) ∂Γε T ij θ ν ds, (57)
where we took exactly the same steps as in (43).

Summing up, the derivative of the integral on s ε can be written as

d dp sε(y p ) T ij u j dS p=0 = sε(y) [(T ij u j ) * + T ij u j D θ ]dS = uj (y) + ∂u j ∂a c ij (y) + u j (y)a b ij ε + d ij + O(ε) -u j (y) ∂Γε T ij θ ν ds. (58)
We can now turn our attention to the integral on Γ ε in (47). To compute the derivative with respect to p, for fixed ε, we can employ again the general relation [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF], simplified by the fact that

θ • n = 0 on Γ ε d dp Γε T ij u j dS p=0 = Γε T ij uj dS + Γε Ṫij u j dS + ∂Γε T ij u j θ ν ds = Γε T ij uj dS -a i Γε (Σ ijk, n k )u j dS + u j (y) ∂Γε T ij θ ν ds + ∂Γε T ij [u j -u j (y)]θ ν ds = Γε T ij uj dS + a Γε V ij u j dS + u j (y) ∂Γε T ij θ ν ds + a ∂Γε T ij [u j -u j (y)]ν ds + o(1), (59)
where Ṫij = a V ij .

It is worth noting that when i = j, and contrary to the potential case (cfr. (45)), the final line integral in the above expression may not be zero in the limit. In Appendix B it is shown that, when y is between two boundary elements with border of tangent a, we obtain the following result (i = j)

lim ε→0 ∂Γε T ij [u j -u j (y)](a ν )ds = (1 -2ν) 12π(1 -ν) ( ∂u ∂b 1 × b 1 ) i -( ∂u ∂b 2 × b 2 ) i , (60)
where, on each side 1 and 2, b = n × a is a unit vector normal to the edge and tangent to the surface.

It is clear that the result in (60) is zero when at y the boundary is smooth and the gradient ∇u is continuous. Smoothness of the boundary is, in general, not sufficient since only ∂u/∂a is continuous across the two elements.

When contributions from s ε and Γ ε are taken together, we observe that in (58) and (59) there are two line integrals that cancel each other. Therefore, we obtain the following final expression, if

u ∈ C 1,α (y) dI dp p=0 = lim ε→0 d dp Γε(y p )+sε(y p ) T ij u j dS p=0 = lim ε→0 * u j (y) sε T ij dS + Γε T ij uj dS + u j (y) sε Ṫij dS + Γε Ṫij u j dS = lim ε→0 uj (y) + ∂u j ∂a sε T ij dS + Γε T ij uj dS + u j (y)a sε V ij dS + a Γε V ij u j dS = uj (y)c ij (y) + ∂u j ∂a c ij (y) + - ∂Ω T ij uj dS + lim ε→0 a Γε V ij u j dS + u j (y) a b ij ε + d ij . (61)
If u ∈ C 1,α (y), the term in (60) must be added on the r.h.s.

Some general remarks

In the analysis presented in this paper the density functions q, u, t j and u j are just sufficiently regular functions that do not necessarily solve a boundary value problem as in the case, e.g., of approximate or numerical solutions. In particular, Green's or Betti's reciprocity identities have never been used.

Moreover, the resulting formulae for the tangential derivatives of the integral I(p) allow for a change of density ( q, u, etc.) induced by the change of the singular point.

Among other things, these results support the approach presented in [START_REF] Guiggiani | Sensitivity analysis for boundary element error estimation and mesh refinement[END_REF] to BEM error estimation through sensitivity analysis as shown in Section 8.

Quite remarkable is the fairly weak requirement posed on the smoothness of u and u j at the singular point y. The presented analysis only needs the density u and its total derivative * u (eq. ( 39)) be C 0,α at the singular points. In other words, we do not require u ∈ C 1,α (y), but only u and its directional derivative ∂u/∂a to be Hölder continuous at y. As already stated, directional derivatives are continuous between two standard boundary elemements provided a is taken tangent to the interelement border. Therefore, the analysis is applicable when the point y is located on an edge and simply continuous elements are used.

Numerical examples

Some numerical computations are carried out to test the theoretical findings presented in this paper.

In particular, we consider formula (61) for the tangential derivative of a double-layer elastic potential, which is probably the most severe test.

Formula (61) covers the general case. However, if we make some simplifying assumptions, the same computation can also be performed by alternative classical means that can be used for comparison.

In particular, if we consider a piece-wise flat surface Γ, the kernel function T ij becomes

T ij (y, x) = - (1 -2ν) 8π(1 -ν) (n i r, j -n j r, i ) r 2 = (1 -2ν) 8π(1 -ν) D ij 1 r ,
and hence, according to (5), the following Cauchy principal value integral can be converted into a line integral

- Γ T ij dS = lim ε→0 Γε(y p ) T ij (y p , x) dS x = (1 -2ν) 8π(1 -ν) lim ε→0 ∂Γε e ijk τ k (x) ds r ,
where ∂Γ ε = ∂Γ ∪ ∂e ε (Figure 1). If the singular point lies on the edge between two panels Γ 1 and Γ 2 we obtain that

I(y) = - Γ T ij dS = 1 -2ν 8π(1 -ν) ∂Γ e ijk τ k (x) ds r + 2e ijk (b 1 -b 2 ) • e k ,
where, on each panel Γ , b = n × a like in (60). If the source point y lies within one panel, the last term vanishes.

In either case, the line integrals on the r.h.s. involve only regular functions since y ∈ ∂Γ. Thus, differentiation with respect to y can be expressed by differentiation under the integral sign. If, as usual in this paper, y p = y(p) = y + a p + o(p), we have that

dI dp p=0 = 1 -2ν 8π(1 -ν) ∂Γ e ijk τ k a (x -y ) r 3 ds. ( 62
)
The regular integral on the r.h.s. can be evaluated by standard quadrature rules.

The last result is precisely what we need to evaluate the derivative in (61) if we assume a constant density u j = 1 (and hence u ∈ C 1,α (y)), with uj = 0. Therefore, we can compare these results to check the general formula (61) which employs the hypersingular integration algorithm described in [START_REF] Guiggiani | Hypersingular boundary integral equations: a new approach to their numerical treatment[END_REF] and [START_REF] Guiggiani | A general algorithm for the numerical solution of hypersingular boundary integral equations[END_REF]. In particular, we can test the exactness of our results regarding the differentiation along the edge. Two flat panels of distorted shape and general orientation were considered (Figure 2). The coordinates of the four nodes of each element are as follows: (1.03, 0.65, 1.03), (0.3, 0.4, 1), (0, 0, 0),

(1, 0.3, -0.1) for one element, and (0.3, 0.4, 1), (0.24, 1.52, 0.8), (0, 1, 0), (0, 0, 0), for the other element. The singular point y = (0.189, 0.252, 0.63) was on the common edge of direction a = (0.268328, 0.357771, 0.894427). A Poisson ratio ν = 0.3 was assumed.

Both methods, i.e., equations (61) and (62), provided the same results up to the seventh significant digit, thus confirming the correctness of the theoretical results presented even when the singular point is located on an edge. The numerical results are summarized in Table 1 for all relevant indices of the kernel function T ij (kernels with equal indices are not singular on flat elements).

For completeness, the same comparison was carried out also for the case of the singular point y located inside one flat element. Again, for more generality, the element had a distorted shape. The coordinates of the four nodes are as in the first element of the former example. The singular point y was set to have intrinsic c oordinates (0.3, -0.2), and the derivative was taken along the direction a = (0.645733, 0.407501, 0.645733). The results obtained with both methods, i.e., equations ( 61) and

(62), agreed within six significant digits and are reported in Table 2.

Sensitivity of approximate boundary element solutions

In the standard boundary element method (BEM), approximate solutions are obtained by collocation at a suitable number of boundary points. Details can be found in any BEM book.

Let, as in [START_REF] Guiggiani | Sensitivity analysis for boundary element error estimation and mesh refinement[END_REF], the whole set of collocation points {y (k) } be dependent upon a single scalar parameter p. Therefore, given the value of p we know the location of each collocation point y (k) (p) on ∂Ω. The value p = 0 can be associated to the standard nodal collocation. Clearly, the approximate BEM solution u j ({y (k) (p)}, x) = u j (p, x) and t j ({y (k) (p)}, x) = t j (p, x) do depend on p through the full set {y (k) (p)} of p-dependent collocation points.

After the system of equations has been solved, the approximate solution is known all over the boundary ∂Ω and we have the following identities

lim ε→0 ∂Ωε(y (k) (p)) T ij (y (k) (p), x)u j (p, x) -U ij (y (k) (p), x)t j (p, x) dS = R k (p) ≡ 0, (63)
where y (k) (p) is any collocation point. Some comments are in order here. First, we must observe that R k is not the residual of the boundary integral equation. Obviously, since R k (p) ≡ 0, also its derivative is identically zero

dR k dp ≡ 0.
The results obtained in this paper allow for a more explicit expression of the above derivative. It is just a matter of collecting the results in (30) and ( 61)

lim ε→0 d dp ∂Ωε(y (k) (p)) T ij (y (k) (p), x)u j (p, x) -U ij (y (k) (p), x)t j (p, x) dS = c ij (y) uj (y) + - ∂Ω T ij uj dS - ∂Ω U ij ṫj dS + c ij (y) ∂u j ∂a + a d ij u j (y) + lim ε→0 Γε a [V ij u j -W ij t j ]dS + a b ij ε u j (y) ≡ 0, ( 64
)
where y = y (k) (p). If u j ∈ C 1,α (y (k) ), the additional term (60) must be added on the r.h.s.

Expression (64) supports the approach presented in [START_REF] Guiggiani | Sensitivity analysis for boundary element error estimation and mesh refinement[END_REF] where the sensitivity of two-dimensional approximate boundary element solutions with respect to the positions of collocation points was obtained. Each row in the system of equations ( 14) in [START_REF] Guiggiani | Sensitivity analysis for boundary element error estimation and mesh refinement[END_REF] is identical to eq. ( 64), provided the BEM solution is used for the density functions u j and t j . (The notation is different: dotted quantities in this paper correspond to starred quantities in the other paper).

If we rewrite the identity (64) in this fashion The functions u and ṫ show that the approximate solution is, in fact, sensitive to the selection of collocation points whenever the same approximate solution fails to satisfy the HBIE for tangential derivatives. The exact solution would have all dotted functions identically equal to zero since the r.h.s.

would be identically zero as well.

In [START_REF] Guiggiani | Sensitivity analysis for boundary element error estimation and mesh refinement[END_REF] it was heuristically suggested to employ suitable norms of these sensitivities to steer a self adaptive mesh refinement process.

Still by integration by parts we already knew that sε G, ij n j dS = ∂Γε G, j (ν i n j -ν j n i )ds.

Therefore, an alternative expression for (46) is lim Since all integrals are weakly singular the limit exists and it does not depend on the selected shape of v ε .

Similar proofs can be provided for the elastic case. Since u, 1 is continuous across the edge, the contribution for the two half-circles add up to zero.

On the other hand, the contributions of terms involving (u, 2 ×b) i will give rise to nonzero free terms in eq. ( 60) which depend on the two derivatives of u along the directions b .

Figure 1 :

 1 Figure 1: Exclusion of the singular point y by a vanishing neighbourhood v ε (y).

Figure 2 :

 2 Figure 2: Flat distorted elements employed in numerical examples.

c

  ij (y) uj (y) + -∂Ω T ij uj dS -∂Ω U ij ṫj dS = -c ij (y) ∂u j ∂a + a d ij u j (y) + lim ε→0 Γε a [V ij u j -W ij t j ]dS + a b ij ε u j (y) . (65)we have on the l.h.s. a standard boundary integral equation for the "sensitivities" u and ṫ, and on the r.h.s. the residual (with minus sign) of the hypersingular boundary integral equation (HBIE) for the tangential derivatives in the a-direction.

  n j dS + Γε (G, j n j ) udS + u(y)a i ∂Γε G, j (ν i n j -ν j n i )ds + ∂Γε [(a i n i )ν j -(a i ν i )n j ]G, j uds + a i Γε G, j (D ji u)dS + a i ∂Γε G, j n j [u -u(y)]ν i ds = lim ε→0 * u(y) sε G, j n j dS + Γε (G, j n j ) udS + a i Γε G, j (D ji u)dS + ∂Γε G, j ν j (a i n i )[u -u(y)]ds , (A3)where the final line integral clearly tends to zero if θ • n = 0 and u ∈ C 0,α (y).Further manipulations are possible since sε G, j n j dS = -Γε G, j n j dS anda i Γε G, j (D ji u)dS = Γε [(a i u, i )G, j n j -(a i n i )G, j u, j ]dS,where a i n i = O(ε α ). Moreover, we can seta i u, i = a i [u, i -u, i (y)] + a i u, i (y) = a i [u, i -u, i (y)] + u(y) + a i u, i (y) -u(y) = a i [u, i -u, i (y)] + * u(y) -u(y).Summing up, the following regularized expression for (A3) or (46) is obtainedlim u(y)] + a i [u, i -u, i (y)] G, j n j -(a i n i )G, j u, j dS = ∂Ω[ u -u(y)] + a i [u, i -u, i (y)] G, j n j -(a i n i )G, j u, j dS. (A4)

T

  ij [u j -u j (y)](a ν )ds = A[(δ j1 δ m2 + δ j2 δ m1 )δ i3 -(δ i1 δ m2 + δ i2 δ m1 )δ j3 ]u j,m (y) = A[(δ i3 δ j1 -δ i1 δ j3 )δ m2 + (δ i3 δ j2 -δ i2 δ j3 )δ m1 ]u j,m (y) = A[e kij e k31 u j,2 (y) + e kij e k32 u j,1 (y)] = A[(u, 2 ×b) i -(u, 1 ×a) i ],(B1)where A = (1 -2ν)/(12π(1 -ν)) and the derivatives of u are intended in the local frame (a, b, n).

Table 1 :

 1 Derivatives of strongly singular integrals: pole along the edge between two elements.

	indices (i, j)	derivative
	(1,2)	-0.00277012
	(2,3)	-0.0719483
	(3,1)	0.0590433

Table 2 :

 2 Directional derivatives of strongly singular integrals: pole within one element.

	indices (i, j) derivative
	(1,2)	0.113809
	(2,3)	-0.0711122
	(3,1)	0.0107791

  Changing p means selecting a new set of collocation points and obtaining, by collocation and through the solution of a new system of equations, a different numerical solution which has again zero residual at each new point y (k) (p).

The residual would be obtained with a fixed approximate solution and a changing singular point, whereas both depend on p in (63). In fact, R k (p) is, by construction, identically zero for any p.

Acknowledgments

This work was carried out while M. Bonnet was visiting the University of Pisa (Italy) sponsored by the French Academy of Sciences. Support to M. Guiggiani was provided by MURST.

Appendix A: Existence of the limits

The integrals in (24) whose value might be affected by the shape of the exclusion neighbourhood are those with more singular kernels. For simplicity, the analysis is presented for smooth boundary at y.

However, the same path can be followed when y lies on an edge and q is only piecewise C 0,α . Basically, using (2) and then integrating by parts using eq. ( 6) we have that where κ is the mean curvature. In the last equation, the original singular integrals have been converted

into a weakly singular integral on Γ ε . Thus the limit investigated exists and its value does not depend on the shape of v ε (y). Now, let us consider eq. ( 46). We recall that the third line in (46) follows from the first or second line only when a circular shape of uniform radius ε is taken for the vanishing neighbourhood v ε . We have therefore to show that this choice for v ε is not relevant.

Integrating by parts as in [START_REF] Guiggiani | Error indicators for adaptive mesh refinement in the Boundary Element Method-A new approach[END_REF] we have that

Appendix B: Proof of eq. ( 60)

Let the singular point y be on an edge of tangent a. On each side we define a unit vector b , = 1, 2 b = n × a normal to the edge and tangent to the surface. Clearly, (a, b , n ), = 1, 2, are two Cartesian orthonormal bases.

Let (a, b, n) = (a, b 1 , n 1 ) be the basis with b pointing out of the edge. We have on the half part ∂Γ ε of ∂Γ ε lying on the corresponding element

where polar coordinated centered at y have been introduced.

Assume that on each side of the edge (i.e., for 0 ≤ θ ≤ π here), u has a C 0,α gradient. Then

Then, one has