Modeling Electronic Skin Response to Normal Distributed Force

The reference electronic skin is a sensor array based on PVDF (Polyvinylidene fluoride) piezoelectric polymers, coupled to a rigid substrate and covered by an elastomer layer. It is first evaluated how a distributed normal force (Hertzian distribution) is transmitted to an extended PVDF sensor through the elastomer layer. A simplified approach based on Boussinesq’s half-space assumption is used to get a qualitative picture and extensive FEM simulations allow determination of the quantitative response for the actual finite elastomer layer. The ultimate use of the present model is to estimate the electrical sensor output from a measure of a basic mechanical action at the skin surface. However this requires that the PVDF piezoelectric coefficient be known a-priori. This was not the case in the present investigation. However, the numerical model has been used to fit experimental data from a real skin prototype and to estimate the sensor piezoelectric coefficient. It turned out that this value depends on the preload and decreases as a result of PVDF aging and fatigue. This framework contains all the fundamental ingredients of a fully predictive model, suggesting a number of future developments potentially useful for skin design and validation of the fabrication technology.


Introduction
Touch-sensitive electronic skin (e-skin) provides information on contact events occurring on its surface and can be used in a variety of contexts involving tactile interactions [1][2][3][4][5]. In recent decades, different tactile sensing technologies have been extensively developed especially for autonomous robotics [6][7][8]. However, there are a lack of models of the overall skin behavior, which would allow for better skin design and for the validation of the fabrication technology. Some interesting attempts in this direction have been made [9][10][11][12] e.g., to show how grasping information can be extracted from strain sensors beneath a compliant skin using simplified solid mechanics models and basic contact theory [13]. However, consideration of skin mechanics has been infrequent in the design and use of tactile sensors.
From a system perspective, e-skin includes stack-wise arrangements of functional and structural materials together with adequate interface electronics to read sensor signals [14]. In this paper, the reference architecture is a basic multilayer involving a rigid substrate, a Polyvinylidene fluoride (PVDF) piezoelectric polymer sensor array and an elastic layer on top for stress transmission and sensor protection. It is important to note that the whole theoretical analysis considers a single sensor skin. The design of a sensor array can be optimized based on the single sensor analysis, which studies how a distributed (Hertzian) force is transmitted to an extended sensor through an elastic layer. Indeed, needless to say, the presence of an array of sensors rather than a single sensor does not modify the analysis.
where, D 3 is the charge density on the sensor surface, T 3 is the normal stress component (i.e., the pressure acting on the bottom of the protective layer), and d 33 is the piezoelectric coefficient. As PVDF sensors directly convert mechanical stimuli into charge, electronic circuits for data acquisition are based on charge amplifiers [18]. Equation (1) does not include the electric field across the PVDF sensor as it is assumed to be negligible due to the virtual ground at the operational amplifier inverting input. Hence, measuring the charge density on the sensor surface provides a direct measure of the normal stress acting on the sensor surface. Assuming circular sensors, the total sensor charge measured by the charge amplifier is given by: where r T is the sensor radius andT 3 is the normal stress component T 3 averaged over a single extended sensor. small film thickness (e.g., 100 µ m for commercial films), when integrated on a rigid substrate PVDF sensors work in thickness mode, such that [17]: where, D3 is the charge density on the sensor surface, T3 is the normal stress component (i.e., the pressure acting on the bottom of the protective layer), and d33 is the piezoelectric coefficient. As PVDF sensors directly convert mechanical stimuli into charge, electronic circuits for data acquisition are based on charge amplifiers [18]. Equation (1) does not include the electric field across the PVDF sensor as it is assumed to be negligible due to the virtual ground at the operational amplifier inverting input. Hence, measuring the charge density on the sensor surface provides a direct measure of the normal stress acting on the sensor surface. Assuming circular sensors, the total sensor charge measured by the charge amplifier is given by: where T r is the sensor radius and 3 T is the normal stress component T3 averaged over a single extended sensor. The protective layer, which is usually polymer-based (e.g., PDMS), implements a mechanical filtering of the tactile stimulus applied at the skin surface and distributes the mechanical stimulus onto the sensor array below. The thickness of the soft layer is a very important parameter in the design of a flexible skin. In the case of the present study-where PVDF is placed at the device bottom and works in thickness mode-thicker layer means weaker sensor response signals. However, a thicker elastomer layer exalts its surface deflection under contact forces and warrants more accuracy in the measurement of changes of curvature if sensors are reduced to tiny compliant strips and placed near the outer elastomer surface.

Electronic Skin Model
Equation (2) relates PVDF sensor output (i.e., the measured charge Qmeas) to the normal stress component T3 averaged over the sensor (i.e., 3 T ). Scope of the following two sections is to retrieve the average normal stress component 3 T transmitted to a single extended tactile sensor as a function of the normal force F3 applied at the surface of the skin protective layer. This would allow to estimate the sensor charge Qmeas (output) as a function of the normal force F3 (input), which leads to an estimation of the system transfer function FRF, defined as the ratio between the Fourier transforms of the output charge and of the input force. In this paper, we only focus on normal contact forces.
Basic assumption of the proposed model (which can be relaxed in future steps) is that the e-skin protective layer can be treated as an incompressible elastic medium (Poisson ratio sufficiently close to 0.5). This is appropriate for an elastomer (ν = 0.48), but not for other materials (e.g., a foam).
In the present section we make an additional major assumption: we assume that our system, namely an elastic medium consisting of a layer with finite thickness, is part of an elastic half-medium bounded by the free surface on which the external force is applied. As shown below, this assumption is quite useful, allowing for an analytical solution of the mathematical problem in closed form, based on Boussinesq's classical solution [19].
This approach provides a qualitative picture of the response of the e-skin to external forces that lets the relevant physical parameters emerge naturally from the analysis. However, the price paid is The protective layer, which is usually polymer-based (e.g., PDMS), implements a mechanical filtering of the tactile stimulus applied at the skin surface and distributes the mechanical stimulus onto the sensor array below. The thickness of the soft layer is a very important parameter in the design of a flexible skin. In the case of the present study-where PVDF is placed at the device bottom and works in thickness mode-thicker layer means weaker sensor response signals. However, a thicker elastomer layer exalts its surface deflection under contact forces and warrants more accuracy in the measurement of changes of curvature if sensors are reduced to tiny compliant strips and placed near the outer elastomer surface.

Electronic Skin Model
Equation (2) relates PVDF sensor output (i.e., the measured charge Q meas ) to the normal stress component T 3 averaged over the sensor (i.e.,T 3 ). Scope of the following two sections is to retrieve the average normal stress componentT 3 transmitted to a single extended tactile sensor as a function of the normal force F 3 applied at the surface of the skin protective layer. This would allow to estimate the sensor charge Q meas (output) as a function of the normal force F 3 (input), which leads to an estimation of the system transfer function FRF, defined as the ratio between the Fourier transforms of the output charge and of the input force. In this paper, we only focus on normal contact forces.
Basic assumption of the proposed model (which can be relaxed in future steps) is that the e-skin protective layer can be treated as an incompressible elastic medium (Poisson ratio sufficiently close to 0.5). This is appropriate for an elastomer (ν = 0.48), but not for other materials (e.g., a foam).
In the present section we make an additional major assumption: we assume that our system, namely an elastic medium consisting of a layer with finite thickness, is part of an elastic half-medium bounded by the free surface on which the external force is applied. As shown below, this assumption is quite useful, allowing for an analytical solution of the mathematical problem in closed form, based on Boussinesq's classical solution [19].
This approach provides a qualitative picture of the response of the e-skin to external forces that lets the relevant physical parameters emerge naturally from the analysis. However, the price paid is the inability of the approach to account for the actually finite thickness of the elastomer as well as to impose the boundary condition at the rigid substrate. In the next section we then relax the latter assumption and solve the problem numerically for the actual physical configuration. This will allow us to obtain a quantitatively more reliable picture suitable for experimental validation.

Effect of Sensor Size on Sensor Response to a Single Normal Point Force
Consider a single sensor of radius r T at a fixed position on a rigid surface (not necessarily plane) with unit normal vector n. Assume a point force F is applied at a given position on the outer surface of the layer. The surface is coated with an elastomer layer of constant thickness h (Figure 2a). A stress field is then generated in the elastomer and transmitted to the sensor. Let T denote the stress tensor. The stress vector acting on the sensor reads n·T. This stress response can be conveyed into an appropriate circuit and retrieved by an electronic device. For example, PVDF sensors directly convert the T 3 stress component into charge, as indicated by Equation (1). the inability of the approach to account for the actually finite thickness of the elastomer as well as to impose the boundary condition at the rigid substrate. In the next section we then relax the latter assumption and solve the problem numerically for the actual physical configuration. This will allow us to obtain a quantitatively more reliable picture suitable for experimental validation.

Effect of Sensor Size on Sensor Response to a Single Normal Point Force
Consider a single sensor of radius T r at a fixed position on a rigid surface (not necessarily plane) with unit normal vector n. Assume a point force F is applied at a given position on the outer surface of the layer. The surface is coated with an elastomer layer of constant thickness h (Figure 2a). A stress field is then generated in the elastomer and transmitted to the sensor. Let T denote the stress tensor. The stress vector acting on the sensor reads n·T. This stress response can be conveyed into an appropriate circuit and retrieved by an electronic device. For example, PVDF sensors directly convert the T3 stress component into charge, as indicated by Equation (1). In order to determine the relation between a point load F applied on the outer surface and the stress at a given point inside the cover layer, we may get advantage of the solution of the so-called Boussinesq's problem [19]. This problem considers an elastic half-medium bounded by a surface on which a point force is applied. For such a configuration the stress field determined by Boussinesq reads: where r is the radial distance of the generic point of the medium from the application point of the load F, all bold-faced symbols represent tensors or vectors, ek is the unit vector in the k-direction and  is the symbol of tensor product. Let us next assume that Equation (3) may be applied to our configuration, approximating our finite thickness layer with Boussinesq's semi-infinite medium. Hence, the stress tensor at the center of the sensor will be evaluated by Equation (3), with r separation vector of the force application point from the center of the sensor area. The advantage of Equation (3) is its simplicity, as the stress is uniaxial in the radial direction, and its independence of elastic parameters.
In this paper, we develop Equation (  In order to determine the relation between a point load F applied on the outer surface and the stress at a given point inside the cover layer, we may get advantage of the solution of the so-called Boussinesq's problem [19]. This problem considers an elastic half-medium bounded by a surface on which a point force is applied. For such a configuration the stress field determined by Boussinesq reads: where r is the radial distance of the generic point of the medium from the application point of the load F, all bold-faced symbols represent tensors or vectors, e k is the unit vector in the k-direction and ⊗ is the symbol of tensor product. Let us next assume that Equation (3) may be applied to our configuration, approximating our finite thickness layer with Boussinesq's semi-infinite medium. Hence, the stress tensor at the center of the sensor will be evaluated by Equation (3), with r separation vector of the force application point from the center of the sensor area. The advantage of Equation (3) is its simplicity, as the stress is uniaxial in the radial direction, and its independence of elastic parameters.
In this paper, we develop Equation (3) considering the sole T 3 stress component on the bottom of the elastic cover of thickness h, as received by a PVDF sensor working in thickness mode Equation (1). Lettingr be the radial distance of the point where the force is applied from the sensor center projected on the outer surface (see Figure 2a), we have r 2 =r 2 + h 2 and e r = sinα (e 1 cosβ + e 2 sinβ) − e 3 cosα (see Figure 2b), where cosα = h/r, sinα =r/r.
For a vertical contact force F3 (F1 = F2 = 0), Equation (4) reduces to where r h  .  As the sensor is not point-like, it is convenient to average the stress over the surface of the sensor. With the notations of Figure 3  , i.e., the dimensionless distance of the sensor center from the projection of the force application point on the sensor plane. Note that this is the first important result of the analysis: indeed, γ is a measure of how the normal point force F3 is transmitted to the sensor through the elastic layer leading to an average normal stress 3 T acting on the sensor.
In Figure 4, the dependence of  As the sensor is not point-like, it is convenient to average the stress over the surface of the sensor. With the notations of Figure 3 we find: Asr 2 = r 2 +r 2 − 2 rrcosθ, λ =r h and r h is defined as λ. In Equation (6), the γ (ˆr h , r T h ) coefficient includes the double integral, which can be solved numerically as a function ofr h , i.e., the dimensionless distance of the sensor center from the projection of the force application point on the sensor plane. Note that this is the first important result of the analysis: indeed, γ is a measure of how the normal point force F 3 is transmitted to the sensor through the elastic layer leading to an average normal stressT 3 acting on the sensor.
In Figure 4, the dependence of γ (ˆr h , r T h ) onr h is plotted for different skin designs, measured by the ratio of sensor radius to the thickness of the layer, r T h . The reference black curve corresponds to the case of "point-like" sensor (r T = 0). The color legend is included in the figure caption.   (6)) is plotted versus the distance r from sensor center (see Figure 2 for notations). Different curves are associated with different sensor sizes as measured by the From an analysis of Figure 4, it follows that, for sufficiently small values of T r h (say 0.2, yellow curve), the problem can be treated as if the sensor was point-like. In the other cases, the problem can still be treated as if the sensor was point-like provided that the force is located at a sufficient distance from the sensor center. For example, if T r 1 h  (magenta curve) in order to treat the sensor as point-like the force is to be located at a distance that is equal/larger than the sensor radius. Figure 4 is an important first achievement of this analysis, as it can be used as a practical tool for skin design. As a matter of fact, it allows one to optimize the system geometry (rT/h and sensor pitch) to achieve a desired skin resolution, which is related to the spatial concentration of the mechanical stress information around a single sensor.

Effect of Normal Force Distributed Over Circular Contact Area on Point-Like Sensor Response
Consider normal force acting on a circular contact area of radius a with Hertzian pressure distribution   qr [20], which realistically describes normal contact by a rigid spherical indenter ( Figure 5): The reader is warned that, unlike in Section 2.3.1, here r denotes a radial coordinate with origin at the center of the contact area (see sketch in Figure 5).  (6)) is plotted versus the distancer from sensor center (see Figure 2 for notations). Different curves are associated with different sensor sizes as measured by the dimensionless From an analysis of Figure 4, it follows that, for sufficiently small values of r T h (say 0.2, yellow curve), the problem can be treated as if the sensor was point-like. In the other cases, the problem can still be treated as if the sensor was point-like provided that the force is located at a sufficient distance from the sensor center. For example, if r T h = 1 (magenta curve) in order to treat the sensor as point-like the force is to be located at a distance that is equal/larger than the sensor radius. Figure 4 is an important first achievement of this analysis, as it can be used as a practical tool for skin design. As a matter of fact, it allows one to optimize the system geometry (r T /h and sensor pitch) to achieve a desired skin resolution, which is related to the spatial concentration of the mechanical stress information around a single sensor.

Effect of Normal Force Distributed Over Circular Contact Area on Point-Like Sensor Response
Consider normal force acting on a circular contact area of radius a with Hertzian pressure distribution q(r) [20], which realistically describes normal contact by a rigid spherical indenter ( Figure 5): The reader is warned that, unlike in Section 2.3.1, herer denotes a radial coordinate with origin at the center of the contact area (see sketch in Figure 5). = q(r)r dr dθ , Using Equation (7) into Equation (8) allows for retrieving 0 q as follows: The radius a of the imprint is related to the applied load F3 by the equation [20]: where R is the radius of the spherical indenter, E is the elastic modulus of the protective layer and ν its Poisson coefficient.
In the following, we describe the effect of this distributed force on a point-like sensor, which is aligned with the contact center. Considering other sensor positions is far from the scopes of this paper.
Having in mind that Equation (3) is the Green function of the problem, the normal stress T3 generated on a point-like sensor located at a depth h below the center of the contact area can be calculated by:  The overall contact force is the result of integrating q(r) over the contact area ( Figure 5b) Using Equation (7) into Equation (8) allows for retrieving q 0 as follows: The radius a of the imprint is related to the applied load F 3 by the equation [20]: where R is the radius of the spherical indenter, E is the elastic modulus of the protective layer and ν its Poisson coefficient. In the following, we describe the effect of this distributed force on a point-like sensor, which is aligned with the contact center. Considering other sensor positions is far from the scopes of this paper.
Having in mind that Equation (3) is the Green function of the problem, the normal stress T 3 generated on a point-like sensor located at a depth h below the center of the contact area can be calculated by: 2r h 3 dr By introducing new integration variable λ =ˆr h , defining a * = a h and using Equation (9) for q 0 we get: In Equation (12), the δ ( a h ) coefficient includes the double integral, which can be solved numerically as a function of a h , i.e., the contact radius (depending on F 3 through Equation (10)) normalized to the layer thickness. Note that this is the second important result of the analysis: indeed, δ is a measure of how the force F 3 (distributed with Hertzian pressure distribution q(r)) is transmitted to the point sensor through the elastic layer leading to the stress T 3 on the sensor. In Figure 6, δ is plotted as a function of a h .
In Equation (12), the () h  coefficient includes the double integral, which can be solved numerically as a function of a h , i.e., the contact radius (depending on F3 through Equation (10)) normalized to the layer thickness. Note that this is the second important result of the analysis: indeed, δ is a measure of how the force F3 (distributed with Hertzian pressure distribution   qr ) is transmitted to the point sensor through the elastic layer leading to the stress T3 on the sensor. In Figure 6,  is plotted as a function of a h . Figure 6. The proportionality coefficient δ between normal stress T3 on the sensor and overall contact force F3 (see Equation (12)) is plotted versus the imprint radius a (contact size) scaled by the elastomer thickness h (see Figure 5 for notations). Note that the point-like sensor is aligned with the contact center.

Combination of the Two Contributions: Effect of Distributed Normal Force on the Response of an Extended Sensor
Starting from the case study described in Section 2.3.1, if the force is distributed with Hertzian pressure distribution (Figure 7), Equation (6) can be written as: Note that the function , which weighs the contribution of each point force dF3 to the average normal stress T3 is a function of the radial distance of the point force dF3 from the projection of the sensor center on the outer surface. Recall that this distance coincides with the radial coordinate r employed to integrate over the contact area. Hence, the function γ must be included inside the integral, as shown by Equation (13).

By introducing new integration variable
r h  and defining a a* h  , we get: Figure 6. The proportionality coefficient δ between normal stress T 3 on the sensor and overall contact force F 3 (see Equation (12)) is plotted versus the imprint radius a (contact size) scaled by the elastomer thickness h (see Figure 5 for notations). Note that the point-like sensor is aligned with the contact center.

Combination of the Two Contributions: Effect of Distributed Normal Force on the Response of an Extended Sensor
Starting from the case study described in Section 2.3.1, if the force is distributed with Hertzian pressure distribution (Figure 7), Equation (6) can be written as: Note that the function γ (ˆr h , r T h ), which weighs the contribution of each point force dF 3 to the average normal stress T 3 is a function of the radial distance of the point force dF 3 from the projection of the sensor center on the outer surface. Recall that this distance coincides with the radial coordinater employed to integrate over the contact area. Hence, the function γ must be included inside the integral, as shown by Equation (13).
By introducing new integration variable λ =ˆr h and defining a * = a h , we get: T r T r a As above, for each chosen r T h , σ ( a h , r T h ) includes the double integral that must be solved numerically as a function of a h , i.e., the contact radius scaled by the layer thickness. This is the most important result of this analysis: indeed, σ is a measure of how the normal distributed force F 3 is transmitted to the sensor through the elastic layer leading to an average normal stressT 3 acting on the sensor.
In Figure 8, we plot σ ( a h , r T h ) as a function of a h . The reference black curve corresponds to the case of distributed force centered on point-like sensor and, as expected, it coincides with the function δ ( a h ). On the other hand, in the case of normal point force aligned with the sensor center (a = 0), we find on the y-axis are the same). From an analysis of Figure 8, it follows that, for sufficiently small values of r T h (say 0.2, yellow curve), the problem can be treated as if the sensor was point-like. In the other cases, the problem can still be treated as if the sensor was point-like provided the contact radius is sufficiently large. For example, if r T h = 1 (magenta curve) in order to treat the sensor as point-like the force is to be distributed over a contact area with radius that is at least twice the layer thickness.
Until this point, the specific type of sensor did not come into play. The previous analysis is completely independent of the specific sensor type and only describes the role of the elastomer layer in stress transmission. The specific sensor type (i.e., PVDF piezoelectric polymer) is now included into the picture. Using Equation (14), the charge response of a single extended sensor to distributed normal force centered on the sensor can be estimated as: For a given skin geometry Equation (15) allows one to calculate the effective piezoelectric coefficient d 33 of each PVDF sensor, once the charge Q 3 and the force F 3 have been measured. Comparison with the expected value of d 33 [21] allows one to validate sensor functioning and skin technology. Comparison with the expected value of d33 [21] allows one to validate sensor functioning and skin technology.  (14)) is plotted versus the imprint radius a (contact size) scaled by the elastomer thickness h (see Figure 7 for notations). Note that the applied force is centered on the sensor.  Figure 8. The proportionality coefficient σ between average normal stressT 3 on the sensor and overall (Hertzian) contact force F 3 (see Equation (14)) is plotted versus the imprint radius a (contact size) scaled by the elastomer thickness h (see Figure 7 for notations). Note that the applied force is centered on the sensor. Different curves are associated with different sensor sizes as measured by the dimensionless parameter r T h (yellow r T h = 0.2; green r T h = 0.5; black dotted r T h = 0.6-case study reported in the experimental session; blue r T h = 0.7; magenta r T h = 1; light blue r T h = 1.5; red r T h = 2).

FEM Simulations
We now remove the severe assumption whereby the elastomer may be modeled as an elastic half-medium. In other words, we consider an elastic incompressible medium consisting of a layer of finite thickness h, length l and width b (Figure 9). Length and width of the layer have been chosen arbitrarily, with the sole requirement of the elastomer sides being distant "enough" from the considered sensor such to justify the assumption that the lateral boundaries do not affect the stress field acting on the sensor significantly.
The free surface is assumed to be subject to an external Hertzian pressure distribution Equation (7), with q 0 small enough to lead to small amplitude deformations. The lower boundary is assumed to be rigid, while the side walls are free. The solution of the problem is sought numerically, with the help of the code COMSOL Multiphysics. Figure 9 illustrates an example of the results of FEM simulations of the elastomer subject to a normal force distributed over circular contact area (black line, radius a Equation (10)) with Hertzian pressure distribution. The employed value for the elastomer modulus E is the result of experimental characterization of the elastic layer of the reference e-skin as discussed in Section 3.2.1.
Simulations have been initially performed by assigning fixed values of ν and E (see Section 3.2.1). Moreover, for any given r T /h, the value of a/h has been changed by arbitrarily playing with the two parameters F 3 and R in Equation (10). The output of the simulation is reported in Figure 10, where it is compared with the simplified analytical solution based on Boussinesq's approach. Figure 10 shows that, under the conditions of the configuration investigated in numerical simulations, the numerical output follows a similar qualitative trend but it differs quantitatively from the analytical results discussed in the previous section. This is mainly due to the boundary condition imposed by the rigid substrate. However, note that if the substrate is moved at a depth higher than h, the two solutions do approach each other. This has been verified, performing further simulations where the substrate was set at depths of 10 mm and 30 mm. It turns out that-for r T /h = 0.2 and a/h = 0.289-the relative error of the analytical solution with respect to the numerical one decreases from about 40% to about 10% at 10 mm depth and 1% at 30 mm depth. This suggests that for sensors embedded into a thick protective layer the analytical solution would be quite appropriate. simulations where the substrate was set at depths of 10 mm and 30 mm. It turns out that-for rT/h = 0.2 and a/h = 0.289-the relative error of the analytical solution with respect to the numerical one decreases from about 40% to about 10% at 10 mm depth and 1% at 30 mm depth. This suggests that for sensors embedded into a thick protective layer the analytical solution would be quite appropriate.    It is also important to note that Figure 10 displays apparently irregular oscillations of the response for any given rT/h. This has prompted us to identify what causes these oscillations. Indeed, at a more careful examination, it turns out that, in the finite case, σ depends on an additional parameter, besides a/h and rT/h. This may be readily appreciated noting that, on physical ground, ignoring the effects of the sidewalls, the response of the system may be assumed to depend on the following dimensional quantities: h, rT, E, F and R. With the help of simple dimensional arguments, one may then envisage the following dimensionless relationship: It is also important to note that Figure 10 displays apparently irregular oscillations of the response for any given r T /h. This has prompted us to identify what causes these oscillations. Indeed, at a more careful examination, it turns out that, in the finite case, σ depends on an additional parameter, besides a/h and r T /h. This may be readily appreciated noting that, on physical ground, ignoring the effects of the sidewalls, the response of the system may be assumed to depend on the following dimensional quantities: h, r T , E, F and R. With the help of simple dimensional arguments, one may then envisage the following dimensionless relationship: Note that, in view of Equation (10), the parameter R/h is equivalent to the parameter a/h previously used. Below, we denote the additional parameter F/(R 2 × E) by L. The above argument suggests that the plot of Figure 10 must be modified such that each line corresponding to given r T /h be replaced by a strip of lines each associated with the same value of r T /h but a distinct value of L. The output of the simulations is therefore organized including the dependence of σ on the parameter L. Figure 11 illustrates the results for r T /h = 0.6, which corresponds to the geometry of the real e-skin prototype employed for experimental tests presented in the following section. Note that the figure confirms that distinct curves are associated with different values of L. Similar analysis can be extended to all values of r T /h and analogous strips of lines would be obtained. However, in Figure 11 we have restricted ourselves to a range of values of the parameter L of direct physical relevance for the tactile application.
L. Figure 11 illustrates the results for rT/h = 0.6, which corresponds to the geometry of the real e-skin prototype employed for experimental tests presented in the following section. Note that the figure confirms that distinct curves are associated with different values of L. Similar analysis can be extended to all values of rT/h and analogous strips of lines would be obtained. However, in Figure 11 we have restricted ourselves to a range of values of the parameter L of direct physical relevance for the tactile application. Figure 11. A comparison between Boussinesq's analytical model for the half-space (dotted line) and the numerical COMSOL simulations for the finite case (markers) is reported for rT/h = 0.6, corresponding to the geometry of the e-skin prototype employed in the experimental section. As in Figures 8 and 10, σ vs. a/h is plotted. The role of the parameter L is included and confirms that distinct curves are associated with different values of L.

Experimental Results and Discussion
A series of tests has been finally performed on an electronic skin based on arrays of PVDF transducers, applying a normal contact force on the e-skin surface by a rigid spherical indenter (Hertzian pressure distribution). For the sake of convenience, the normal contact force was aligned with the sensor center, hence a distinct run was needed for each sensor. Experiments allowed measurement of the system response function FRF. We recall that FRF corresponds to the ratio between the Fourier transform of the output charge and that of the input force.
Of course, FRF can also be predicted using the theoretical model. However, the latter prediction depends on the effective piezoelectric coefficient d33 of each PVDF sensor, a quantity that was a priori unknown in the present investigation. Hence, comparison between theoretical predictions and experimental observations has been employed to infer the value of the effective piezoelectric coefficient d33 that leads to best fit. In other words, at the present stage this model may be described as post-dictive. It will become pre-dictive once the effective piezoelectric coefficient d33 will be known a-priori.
A delicate issue arises because the PVDF piezoelectric polymer does not read static information, i.e., the model can only be used with dynamic contacts. A dynamic contact has been indeed employed, as explained below, relying on the expectation that such an approach is safe provided Figure 11. A comparison between Boussinesq's analytical model for the half-space (dotted line) and the numerical COMSOL simulations for the finite case (markers) is reported for r T /h = 0.6, corresponding to the geometry of the e-skin prototype employed in the experimental section. As in Figures 8 and 10, σ vs. a/h is plotted. The role of the parameter L is included and confirms that distinct curves are associated with different values of L.

Experimental Results and Discussion
A series of tests has been finally performed on an electronic skin based on arrays of PVDF transducers, applying a normal contact force on the e-skin surface by a rigid spherical indenter (Hertzian pressure distribution). For the sake of convenience, the normal contact force was aligned with the sensor center, hence a distinct run was needed for each sensor. Experiments allowed measurement of the system response function FRF. We recall that FRF corresponds to the ratio between the Fourier transform of the output charge and that of the input force.
Of course, FRF can also be predicted using the theoretical model. However, the latter prediction depends on the effective piezoelectric coefficient d 33 of each PVDF sensor, a quantity that was a priori unknown in the present investigation. Hence, comparison between theoretical predictions and experimental observations has been employed to infer the value of the effective piezoelectric coefficient d 33 that leads to best fit. In other words, at the present stage this model may be described as post-dictive. It will become pre-dictive once the effective piezoelectric coefficient d 33 will be known a-priori.
A delicate issue arises because the PVDF piezoelectric polymer does not read static information, i.e., the model can only be used with dynamic contacts. A dynamic contact has been indeed employed, as explained below, relying on the expectation that such an approach is safe provided that the forcing frequencies fall outside the range of any significant resonance. Resonances barely derive from the stimulated sensor, as they would result from longitudinal waves which have much higher characteristic frequencies (≈100 kHz) than those related to contact (<1 Hz-1 kHz). The order of magnitude of resonance frequencies f associated with longitudinal waves can be calculated from f = v λ , where v is the speed of sound in the elastomer layer (≈1000 m/s [22]) and λ 4 corresponds to the thickness of the elastomer layer (h = 2.5 mm). Resonances may also derive from the sensor itself, but we have previously proven [21] that its charge-force transfer function (d 33 (ω)) is approximately flat in the frequency range of interest for the tactile application (say, 1 Hz-1000 Hz). This notwithstanding, resonances may derive from a variety of additional causes (e.g., movable contacts, contact surface asperities, motor-induced vibrations), which cannot be a-priori controlled. The ultimate solution is then to identify a flat zone in the system response function and perform experiments only in that frequency range.

Experimental Setup
As mentioned in the introduction, the tactile sensing system we are considering basically consists of an array of stress sensors fixed on a rigid substrate and coated with an elastomer layer. An example of a real device of this kind, based on an array of piezoelectric PVDF sensors, is illustrated in [23]. This e-skin prototype integrates 64 ad-hoc screen-printed electrodes and tracks on both sides of a commercial piezoelectric polymer film (Measurement Specialties Inc., Hampton, VA, USA, http://www.meas-spec.com/default.aspx). Sensor radius equals 1.5 mm and pitch between adjacent sensors is 8 mm. As the flexible skin has been glued on a rigid substrate, the tactile stimulus applied on the outer surface of the elastomer (thickness h = 2.5 mm) is received by the sensors as a distribution of normal stresses (T 3 ), which is converted into charge.
It consists of a rigid frame with a lower fixed plate to which an electro-mechanical shaker (Brüel and Kjaer, Naerum, Denmark, Minishaker Type 4810 with Power Amplifier Type 2706) is assembled. A piezoelectric force transducer (Model 208C01, PCB Piezotronics, Depew, NY, USA) is fixed to the moving head of the shaker. The e-skin sample is faced down side and it is mechanically stimulated by the rigid spherical indenter coupled to the mini-shaker. The loading chain is formed, accordingly, by: shaker, force transducer, spherical indenter, and skin sample. All these elements have to be aligned before any test, for the indenter to be centered on a specific sensor. The mini-shaker is controlled through a graphical user interface (GUI) developed with NI LabVIEW on the host PC and NI DAQ data acquisition board is used. A swept sine signal is fed into the shaker (settable parameters: start and end frequencies, number of steps, amplitude). The output signals of sensor charge (response) and force transducer (stimulus) are continuously acquired by a 4-channel PCB Sensor Signal Conditioner 482C54 (ICP and charge mode) and processed in frequency (i.e., their Fourier transforms are calculated by LabVIEW software) to give the system response function FRF at each frequency step. Before running each test, a preload is applied to guarantee indenter-skin contact during the whole mechanical stimulation. It is this preload which is responsible for determining the contact radius a Equation (10), as for all tests the amplitude of the dynamic oscillation is maintained low enough not to affect significantly the contact area. This preload is controlled using a laser to measure the displacement of the flat plate coupled to the shaker and using shaker displacement-force calibration curves. The experimental setup ( Figure 12) is a slightly different version of that described in [21].

Results
In view of the analysis discussed in Section 2.4, it is appropriate to employ only the full numerical solution of the elastic problem for the actual e-skin prototype to fit experimental data.

Characterization of the Elastic Properties of the Elastomer
The compressive Young's modulus of the PDMS elastic layer was measured using an

Results
In view of the analysis discussed in Section 2.4, it is appropriate to employ only the full numerical solution of the elastic problem for the actual e-skin prototype to fit experimental data.

Characterization of the Elastic Properties of the Elastomer
The compressive Young's modulus of the PDMS elastic layer was measured using an electromechanical machine Zwick/Roell Z0.5 (maximum load 500 N, operating with TestXpert II software), on a cylindrical sample of 4 mm diameter and 2.5 mm length (Figure 13). The test was performed at a speed of 10 mm/min, in the displacement control mode.

Results
In view of the analysis discussed in Section 2.4, it is appropriate to employ only the full numerical solution of the elastic problem for the actual e-skin prototype to fit experimental data.

Characterization of the Elastic Properties of the Elastomer
The compressive Young's modulus of the PDMS elastic layer was measured using an electromechanical machine Zwick/Roell Z0.5 (maximum load 500 N, operating with TestXpert II software), on a cylindrical sample of 4 mm diameter and 2.5 mm length (Figure 13). The test was performed at a speed of 10 mm/min, in the displacement control mode. The elastic Young's modulus of the elastomer has been calculated as the slope of the first linear portion of the curve in Figure 13 and it corresponds to 16 [MPa]. This value has been employed for all numerical models in this paper. Note that non-linear elastomer behavior starts at stress values about 2 [MPa]. The elastic Young's modulus of the elastomer has been calculated as the slope of the first linear portion of the curve in Figure 13 and it corresponds to 16 [MPa]. This value has been employed for all numerical models in this paper. Note that non-linear elastomer behavior starts at stress values about 2 [MPa].

Frequency Selection
The first issue is identifying the appropriate range of frequencies of the mechanical stimulus to validate the electronic skin. As said in the introduction to Section 3, this dynamic approach proves safe provided that the range of forcing frequencies falls outside of any significant system resonance.
First tests have been performed to record the system response function (FRF) experimentally over the whole frequency range (1 Hz-1000 Hz). It turns out that: (i) resonances do exist and there characteristic frequencies depend on the preload and the size of the indenter; (ii) a fairly wide range of frequencies exists, where no significant resonance is detected; (iii) the imaginary part of the response function, which accounts for any viscoelastic component of the response, is roughly an order of magnitude smaller than the real (elastic) part. The latter statement is clarified in Figure 14, where the system response function is plotted in the non-resonant range (roughly 1-130 Hz).
Based on these results, hereafter the imaginary part of the response will be ignored and Re will be removed from the notation. In other words, the system is treated as purely elastic. Moreover, each run has been performed stimulating the skin over the non-resonant range and averaging the corresponding response.
where the system response function is plotted in the non-resonant range (roughly 1-130 Hz).
Based on these results, hereafter the imaginary part of the response will be ignored and Re will be removed from the notation. In other words, the system is treated as purely elastic. Moreover, each run has been performed stimulating the skin over the non-resonant range and averaging the corresponding response.

Indenter Selection
Two sets of tests have been performed over the same selected sensor, each set being characterized by a different indenter (radii equal to 4 mm and 10 mm, respectively). For each given indenter, the preload has been changed from a minimum (0.5 N) to a maximum (3 N) value.
As mentioned above, comparison between experimental and numerical values of FRF ( Figure  15), depends on the assumed value for d33. Best fit is obtained assuming that d33 =14 [pC/N].
Note that, with this choice, agreement is fairly satisfactory for the larger indenter, but differences are larger for the smaller indenter in the high preload range. This suggests the possibility that non-linearity of the response is responsible for such disagreement. Indeed, Figure 13 displays the presence of non-linear effects when the stress exceeds a threshold value of about 2 MPa, which is reached in the small indenter test for high preloads (recall that the effective contact radius is a

Indenter Selection
Two sets of tests have been performed over the same selected sensor, each set being characterized by a different indenter (radii equal to 4 mm and 10 mm, respectively). For each given indenter, the preload has been changed from a minimum (0.5 N) to a maximum (3 N) value.
As mentioned above, comparison between experimental and numerical values of FRF (Figure 15), depends on the assumed value for d 33. Best fit is obtained assuming that d 33 =14 [pC/N].
Note that, with this choice, agreement is fairly satisfactory for the larger indenter, but differences are larger for the smaller indenter in the high preload range. This suggests the possibility that non-linearity of the response is responsible for such disagreement. Indeed, Figure 13 displays the presence of non-linear effects when the stress exceeds a threshold value of about 2 MPa, which is reached in the small indenter test for high preloads (recall that the effective contact radius is a (Equation (10))). Furthermore, non-linearity of the sensor itself may contribute. Further work is needed to investigate this issue in depth. (Equation (10))). Furthermore, non-linearity of the sensor itself may contribute. Further work is needed to investigate this issue in depth. Anyhow, this problem has been further analyzed in the sequence of tests performed on the whole sensor array, using the small indenter (see the next section).

Response Function of the Sensor Array
Finally, the whole sensor array has been tested, by stimulating the e-skin surface with the smaller indenter (R = 4 mm) aligned with the center of each selected sensor. As said in Section 3.2.2, each run has been performed over the non-resonant range at small force amplitude (F_dyn = 0.09 N) and the corresponding FRF response has been averaged over that frequency range to get a single value of the response for each sensor. Anyhow, this problem has been further analyzed in the sequence of tests performed on the whole sensor array, using the small indenter (see the next section).

Response Function of the Sensor Array
Finally, the whole sensor array has been tested, by stimulating the e-skin surface with the smaller indenter (R = 4 mm) aligned with the center of each selected sensor. As said in Section 3.2.2, each run has been performed over the non-resonant range at small force amplitude (F_dyn = 0.09 N) and the corresponding FRF response has been averaged over that frequency range to get a single value of the response for each sensor.
Two sets of data have been obtained. The first extensive study corresponds to higher preload (=3 N) to study the non-linearity of the whole system over a large selection of sensors (41). The second study has been performed over a smaller number of sensors (12), to measure an average response function at smaller preloads (=1 N).
Note that the given preload affects the contact radius a (Equation (10)) while the amplitude of the dynamic swept sine force determines the PVDF charge. On the contrary, the dynamic component does not affect the computation of the contact radius, as the dynamic signal amplitude is negligible with respect to the preload. The effective average normal stress component T 3 affecting the response of the PVDF sensor area is the difference between the total stress (produced by the dynamic force amplitude added to the preload) and the stress due to the preload only. Note that-while the response of the PVDF sensor is not affected by static contacts-its piezoelectric coefficient d 33 (Equation (15)) is likely to be dependent on the preload and some non-linearities may be expected for high preload values.
The following parameters have been used to estimate the system transfer function from the numerical model: r T = 1.5 mm, h = 2.5 mm, F_dyn = 0.09 N. The radius of the Hertzian imprint has been calculated from Equation (10), using: ν = 0.5, E = 16 MPa (Section 3.2.1), R = 4 mm, preload = 1 N or 3 N.
However, predictions for the system transfer function depend on the effective piezoelectric coefficient d 33 of each PVDF sensor, a quantity that was a priori unknown in the present investigation as the PDMS elastic layer was already integrated on top. Some initial characterization of the commercial PVDF film was performed around 2010 and results are reported in [21]. However, a systematic measurement of the d 33 coefficient over the whole sensor array is not available. In any case, although such a characterization is done prior to sensor integration into the complex multilayer skin, long-term aging and fatigue may affect sensor behavior and degradation of sensor properties is expected over time. Being able to estimate the piezoelectric d 33 coefficient from the overall system response function is thus a useful tool to measure the reliability of the e-skin device over time, whenever embedded sensors are not accessible anymore for a direct characterization.
Theoretical predictions will therefore be employed below to infer the value of the effective piezoelectric coefficient d 33 that best fits the experimental observations.
The numerical model has been used to extract the dependence of the d 33 coefficient, which measures the piezoelectricity of the PVDF sensor, on the preload. Data associated with the higher preload (=3 N) are well fitted with a d 33 value equal to 22 [pC/N], while data corresponding to the lower preload (=1 N) yield a d 33 value of 14 [pC/N]. These lower d 33 values, associated with the latest runs, can be considered reasonable as the e-skin was subjected to a huge number of stress cycles in the past 4 years, and this has likely led to film degradation and consequent PVDF aging and fatigue. It is important to appreciate that, as already noted above, the proposed model becomes fully predictive with a preliminary measurement of the PVDF d 33 coefficient of each sensor, prior to integration of the elastic layer on top of the sensor array.
Some dispersion of sensor response is observed, which is a measure of the accuracy of indenter positioning and of the reliability of the whole e-skin fabrication technology. Dispersion in sensor behavior may be due to various factors, which can be considered intrinsic to the manufacturing process. These factors may include different point-to-point values for the sensor radius and/or for the local layer thickness, inhomogeneity in PVDF film polarization and quality of the skin integration technology, which may affect sensor working in pure compression mode (sensor bending may occur).
More sophisticated methods to precisely positioning the indenter are likely to be used in the future to further decrease the dispersion of measured sensor outputs. Results are reported in Figure 16.
integration of the elastic layer on top of the sensor array.
Some dispersion of sensor response is observed, which is a measure of the accuracy of indenter positioning and of the reliability of the whole e-skin fabrication technology. Dispersion in sensor behavior may be due to various factors, which can be considered intrinsic to the manufacturing process. These factors may include different point-to-point values for the sensor radius and/or for the local layer thickness, inhomogeneity in PVDF film polarization and quality of the skin integration technology, which may affect sensor working in pure compression mode (sensor bending may occur). More sophisticated methods to precisely positioning the indenter are likely to be used in the future to further decrease the dispersion of measured sensor outputs.

Conclusions
For the considered electronic skin, contact information is transmitted through an elastomer layer to PVDF sensors that only measure pressure.
In a preliminary stage, Boussinesq's approach has been used to solve the direct problem of how a distributed external normal force is transmitted to an extended PVDF sensor. Although this approach relies on the severe assumption to treat the elastomer as an elastic medium filling a half-space, it has allowed us to perform a comprehensive (albeit qualitative) analysis of the response of an extended sensor to a variety of forcing actions. A normal point force applied anywhere on the outer skin surface has been first investigated. The result of this study (γ curve) provides in itself a very useful practical tool for optimizing skin design to achieve a certain skin resolution. A specific validation of this analysis was not included among the scopes of the present work and will be the

Conclusions
For the considered electronic skin, contact information is transmitted through an elastomer layer to PVDF sensors that only measure pressure.
In a preliminary stage, Boussinesq's approach has been used to solve the direct problem of how a distributed external normal force is transmitted to an extended PVDF sensor. Although this approach relies on the severe assumption to treat the elastomer as an elastic medium filling a half-space, it has allowed us to perform a comprehensive (albeit qualitative) analysis of the response of an extended sensor to a variety of forcing actions. A normal point force applied anywhere on the outer skin surface has been first investigated. The result of this study (γ curve) provides in itself a very useful practical tool for optimizing skin design to achieve a certain skin resolution. A specific validation of this analysis was not included among the scopes of the present work and will be the topic of a forthcoming publication. Next, the sensor response to a Hertzian distributed force acting on the skin surface and centered on the considered sensor has been estimated. The sensor charge response is finally written in terms of the σ parameter, which depends on the sensor and contact radii (r T and a, respectively), both scaled by h.
The half-space assumption has then been relaxed. Numerical FEM simulations have therefore been used to validate the model on a layer of finite thickness, length and width. It turned out that the numerical output follows a qualitative trend similar to that of the previous simplified analysis, but it differs quantitatively mainly due to the boundary condition imposed by the rigid substrate. However, if the substrate is moved at increasing depths higher than the depth h at which the sensor is embedded, then the two solutions do approach each other. This suggests that for sensors embedded into a thick protective layer Boussinesq's analytical solution would be quite appropriate. Moreover, in the finite case, σ was found to depend on an additional parameter, L, which is a function of R, F 3 and E. This parameter did not appear in Boussinesq's analytical solution.
Finally, a real skin prototype based on a PVDF sensor array has been tested, using dynamic contacts over a non-resonant frequency range. First tests with different indenter sizes suggested the presence of some non-linear effects, attributed to both the elastomer layer and the PVDF sensor. Performing a systematic analysis over the whole sensor array, the numerical model has allowed estimation of the d 33 coefficient, which turns out to exhibit some dependence on the preload. Fairly low d 33 values are also found and appear to be due to PVDF aging and fatigue.
Although treated herein as a post-dictive model, the framework laid down in the present work contains all the fundamental ingredients of a fully pre-dictive model, suggesting a number of future developments potentially useful for skin design and validation of the fabrication technology. In its present form, the model is valid for incompressible protective layers (ν ≈ 0.5), but note that extension to compressible materials is feasible as a complete, although admittedly much more complex, solution valid for all ν is available [24]. On the other hand, the method is not constrained to the specific type of sensors, given that these sensors give a generic stress measurement n·T at the bottom of the elastomer layer.
Further future developments include extensions to distributed normal forces not necessarily aligned with the sensor center and to the case of multiple normal distributed loads. Modeling the effect of shear loads as well as of viscoelasticity of both the protective layer and the sensors will also deserve some attention.