Materials that are self-sensing via the piezoresistive effect have been widely explored for embedded sensing in aerospace structural composites, pressure sensing in touch pads, diagnostics in biomedical implant technology, and many other applications. In this approach, changes in electrical transport of the material are used as an indicator of stress/strain, pressure, or damage. However, engineers and other users of self-sensing materials are typically not directly interested in the electrical state of the material. Rather, they want to know the underlying mechanical state of the material that gives rise to an observed electrical change. Recovering material condition from electrical observations is referred to as the self-sensing inverse problem (SSIP). Prior work by the author has shown that the SSIP can be solved using an electrical impedance tomography (EIT)-generated conductivity map as an input, but this is undesirable because it requires solving a second inverse problem (i.e., the SSIP) on top of the EIT inverse problem. To that end, a direct formulation for the SSIP is herein presented. In this approach, voltage-current data is directly inverted to find the displacement field without using EIT as an intermediate step. Additionally, the direct SSIP is solved within the primal-dual interior point (PDIPM) framework such that an ℓ1-norm can be used on the regularization term, which promotes sparsity in the solution space. This approach is applied to a representative piezoresistive nanocomposite using the Laplace matrix as regularization to promote a spatially smooth solution. From the displacement field, strains and stresses are calculated and compared to a commercial finite element solution with good accuracy.
]]>Materials with self-sensing properties have been explored in a wide range of applications including pressure or tactile sensing [1][2][3], robotics [4][5][6], strain and damage sensing in polymeric and cementitious structural materials [7][8][9][10][11], sensing skins applied to other non-self-sensing materials [12][13][14], energetic materials [15][16][17], and even biomedical materials [18][19][20][21][22]. The basic idea of self-sensing is that external effects-deformation, pressure, temperature, chemicals, etc.-change the electrical conductivity of the material. Monitoring electrical transport through the material therefore gives insight into its condition. Because of this effect, self-sensing materials have often been paired with conductivity imaging modalities such as electrical impedance tomography (EIT) [23][24][25][26][27].
But self-sensing materials have a problem. That is, users of self-sensing materials are typically not interested in the electrical state of the material. Instead, they want to know the underlying condition of the material that gives rise to some observed change in conductivity. The process of estimating material condition from electrical data is referred to as the self-sensing inverse problem (SSIP). The SSIP was first postulated, solved, and experimentally demonstrated by Tallman and colleagues. Initial work was computational and assumed a known conductivity distribution [28]. Later work experimentally validated the process in a carbon nanofiber (CNF)-modified polyurethane (PU) nanocom- posite for distributed pressure characterization [29] and for failure prediction near a stress concentration in CNF-modified epoxy [30,31]. Deep learning was also used to solve the failure prediction SSIP [32].
Prior work on the SSIP required knowing the deformationinduced conductivity change distribution. This could be obtained (or at least approximated) experimentally via EIT. However, needing to know the conductivity change distribution is undesirable because EIT is itself an ill-posed inverse problem. This means that the SSIP, by using the EIT solution as an input, is necessarily subject to the same assumptions and limitations of EIT. It is instead highly desirable to deduce the mechanics of the material directly from raw voltage-current data without using EIT as an intermediate step. As a nomenclature clarification, the phrase indirect SSIP will herein be used to describe the process of solving the SSIP with EIT as an intermediate step, whereas the phrase direct SSIP will be used to describe the process of solving the SSIP directly from the raw voltage-current data (i.e., without EIT). Thus, the goal of this work is to derive a direct SSIP formulation. Further, this direct formulation will be experimentally validated on CNF-PU pressure sensing data and regularized by a smoothness-promoting prior, the discrete Laplace operator. An ℓ 1 -norm will be used on the regularization term in the direct SSIP in order to promote a sparse solution. In the experimental validation, the direct SSIP is solved to find the full-field displacements of the material. From the displacements, strains are calculated, and elastic properties are used to calculate stresses (i.e., contact pressures on the sensor). These predictions are compared to a commercial finite element code.
The remainder of this manuscript is organized as follows. First, the experimental details are summarized. Second, forward models used in the SSIP are described. Third, the direct SSIP is derived, discretized for application to a finite element mesh, and solved via the primal-dual interior point method (PDIPM). Fourth, results are presented, compared to a commercial finite element code, and discussed. And lastly, a summary and conclu-sions are presented.
For the sake of brevity, only brief summaries of experimental details are provided in this section. Interested readers are directed to the original publications of this data [29,33] for additional information.
CNF-modified PU was used as a representative piezoresistive pressure sensor. To calibrate the piezoresistivity forward model used in this work (see section CNF/PU Piezoresistivity Model), several CNF/PU volume fraction variations had to be produced. This was done at 7.5, 10.5, 12.5, and 15% XT-HHT CNFs (Applied Sciences) by volume. CNFs were initially hand-mixed into ReoFlex 20 PU (SmoothOn), which resulted in a dough-like mixture. The dough-like mixture was next shear mixed in a two-roll mill and spread into molds measuring 25.4 mm × 76.2 mm × 2.5 mm. The mixture cured at room temperature for 16 hours. After the room-temperature cure, it was post-cured at 65 • C for an additional six hours.
Prismatic strips measuring 25.4 mm × 6.35 mm × 2.5 mm were cut from each volume fraction specimen for electrical testing. Electrodes were attached using silver paste, and four-point methods were used to measure the resistance of each strip. From the dimensions of the strips and the measured resistances, the conductivity of each volume fraction was calculated. This data is shown to the left of Figure 1. The red curve in this plot is a curve fit of the piezoresistivity model shown in equation (1) to the data. The fitting process is how the model is calibrated (see section CNF/PU Piezoresistivity Model). Because stresses are calculated, the elastic properties of the CNF/PU are needed. The 10.5% CNF/PU is used in the direct SSIP, so only this volume fraction was tested mechanically. Normal stresses, normal strain, and transverse strain were measured tensile testing such that Young's modulus and Poisson's ratio could be determined. These results are shown in the middle and right of Figure 1. Curve fitting to this data (the red curves in Figure 1) was used to find that the Young's modulus is E = 1.099 MPa and the Poisson's ratio is ν = 0.4875.
Three contact points were applied to the top of a 25.4 mm × 25.4 mm section of 10.5% CNF-modified PU, as shown in Figure 2. Pressure was applied by resting a 1.25 kg mass atop the marbles resulting in a total force of 11.76 N acting over all three marbles. The mass was allowed to rest for several hours before collecting electrical measurements so that any viscoelastic effects had time to dissipate. The CNF/PU square was adhered to an acrylic base during testing to prevent it from moving (i.e., a fixed bottom surface).
Electrical measurements were taken following an adjacent EIT injection-measurement scheme. That is, 2.5 mA DC was injected between adjacent or neighboring electrode pairs, and the resulting voltage was measured on the remaining electrodes with respect to ground. Adjacent inter-electrode voltage differences, which are used in the direct SSIP, were calculated offline. Current was injected using a Keithley 6221 current supply, and voltages were measured using a National Instruments 6368-PXIe data acquisition system.
Forward models are used to predict the electro-mechanical response of the system computationally. That is, the direct SSIP ultimately seeks a displacement field that, when supplied to the forward models, gives rise to the same voltage-current response as observed experimentally. The SSIP requires two forward models: a model to predict the conductivity change of the CNF/PU for a given state of deformation and a model to predict the electrode voltages for a given conductivity distribution and current injection. To better focus on the novel aspects of this work, the derivation of a direct SSIP formulation, only brief overviews of these forward models are given. However, references are provided so that interested readers can acquire more details.
This work makes use of a nanocomposite piezoresistivity modeled developed previously by Tallman and Wang [34] and that is based on a nanocomposite conductivity model developed by Takeda et al. [35]. These references should be consulted for more information these models. The basic idea of this piezoresistivity model is to identify the strain-dependent parameters in the conductivity model developed by Takeda et al. [35], update those parameters for a given strain state, and then recalculate the conductivity for the given strain state. The electrical conductivity of nanocomposite, σ c , is given below by equation (1).
In equation ( 1), σ m is the PU conductivity, σ f is the CNF conductivity (10 5 S/m), P is the percolation probability, v is the filler volume fraction, l f is the CNF length (11.5 µm), d f is the CNF diameter (0.1 µm), λ is the CNF waviness ratio, A is the projected area of overlap between neighboring CNFs (assumed to be equal to the CNF cross-sectional area), t is the average inter-CNF distance between neighboring CNFs, h is Planck's constant, e is the charge of an electron, and ψ is the barrier height felt by tunneling electrons. Properties for the CNFs are acquired from reference [33]. Baseline (i.e., strain-free) values of percolation probability and average inter-filler distance can be found by fitting equation (1) to experimental conductivity versus volume fraction data. This curve fit is shown by the red curve in the left-most plot of Figure 1. Deformation-dependent conductivity is incorporated into equation ( 1) by expressing the filler volume fraction, percolation probability, and average inter-filler distance as functions of strain. Further, these dependencies, originally described in reference [34], are modified for finite strains in terms of the deformation gradient. The explicit forms of these finite strain dependencies are not included in this manuscript, but interested readers are directed to reference [29].
The voltage-current relationship at the electrodes is found by solving the steady-state diffusion equation subject to the complete electrode model (CEM) boundary conditions. The diffusion equation, assuming no internal current sources, is shown below in equation (2).
In this equation, σ c is the conductivity of the nanocomposite and φ is the domain potential. The CEM boundary conditions enforce a voltage drop between the domain and the electrodes due to contact impedance, that charge is conserved through the electrodes, and that current only enters or exits the domain through the electrodes. These conditions written in equations ( 3), (4), and (5), respectively.
In these equations, z l is the contact impedance under the lth of L electrodes, n is an outward-pointing normal vector, V l is the voltage on the lth electrode, E l is the size of the lth electrode, and ∂ Ω is the boundary of the domain. The finite element method is used to solve these equations. All results in this manuscript were generated using code written by the author in Matlab and linear hexahedral elements. The discretized forms of these equations are not shown here for the sake of brevity, but they are thoroughly treated elsewhere [23,36].
The time-difference direct SSIP attempts to find the difference in strain between two times. Such difference methods are often necessary in experimental inverse problems in order to remove systemic errors or noise common to both times. Define the voltage difference vectors y and W as shown below.
Above, V m (•) is a vector of the experimentally measured electrode voltages at time t 1 or t 2 . F (•) is the electrode voltages predicted by the forward model for some conductivity distribution. Further, it is assumed that the conductivity depends on the deformation gradient, F i j , and the two times differ by a change in the deformation gradient, δ F i j . The goal is to find δ F i j to minimize the difference between y and W . To achieve this goal, linearize F (σ c (F i j + δ F i j )) using a Taylor series expansion and the chain rule of differentiation.
Next, plug equation (8) into equation (7). Upon simplification, the following minimization can be stated. In this minimization δ F * i j is a deformation gradient change that satisfies the minimization, R (•) is a regularization term, α is a scalar hyperparameter controlling the amount of regularization. Note that the ℓ 1 -norm is used on the regularization term. This norm is used because it promotes sparsity in the solution. Recalling that only three pressure points are sought, it is physically intuitive that the direct SSIPpredicted displacement field should be sparse (i.e., zero almost everywhere, except where pressure is being applied).
People familiar with EIT will recognize J σ c as the EIT sensitivity matrix. Herein, this derivative is referred to as the conductivity Jacobian, and it is defined in references [23,36]. The other derivative arising from the chain rule of differentiation, the deformation Jacobian, J F i j , is calculated numerically by the twopoint secant method.
To solve the direct SSIP, equation ( 9) needs to be discretized. This allows the direct SSIP to be recast in terms of nodal displacements. To do this discretization, recall the deformation gradient as shown in equation (10). Recall also the finite elementdiscretized form of the displacement field, u ei , as shown in equation (11). Here, u ei corresponds to the displacement of the eth element in the i-direction. δ i j is the Kronecker-delta, N A is the Ath interpolation function (recall linear hexahedral elements used in this work, so A runs up to eight), and d A ei is displacement of the Ath node of the eth element in the i-direction.
Next, the volumetric average of ∂ u ei /∂ x j is calculated for each element as shown in equation (12). Bars are used to denote volumetric averages. The volumetric average is used to eliminate intra-element variations when solving the direct SSIP. The rightmost equality of equation ( 12) leverages the fact that the integral of the sum is equal to the sum of the integrals and the fact that nodal displacements are constant.
To proceed, consider J F i j δ F i j . δ F i j can be replaced with equation ( 12) by noting that the Kronecker-delta terms will cancel when taking the difference between two deformation gradients and substituting d A ei for δ d A ei (i.e., the change between nodal displacements), as shown for the eth element in equation (13).
Equation ( 13) can be simplified more by substituting the derivative form J e F i j into equation ( 13) and moving it inside the summation for the eth element. These changes are shown in equation (14). Examining the right-most part of equation ( 14) (and remembering to sum over repeated indices), it can be seen that this is just a scalar value for each element. The summation over A and the repeated i index can therefore be replaced with linear algebra as the product of a row vector and column vector, such as shown in the lower portion of equation (14).
Elemental contributions,
where the j index is summed over, can assembled into a global system matrix. Call this global assembly of g A ei contributions as G. Further, define the global nodal displacement change vector as δ d. With these changes in mind, the direct SSIP can be rewritten as shown in equation (15). Degrees of freedom corresponding to known Dirichlet boundary conditions are removed before solving equation (15).
Solving the Direct SSIP via the PDIPM The direct SSIP with ℓ 1 regularization requires the PDIPM to solve because of the non-differentiability of the regularization term. The PDIPM introduces a second or dual optimization problem to be solved concurrently with the primal problem (i.e., the direct SSIP) [37][38][39]. The PDIPM is also an iterative process. Therefore, define x = δ d to avoid using δ (δ d) terms. Finding x is the primal problem, whereas finding v is the dual problem. These are shown in equations ( 16) and ( 17), respectively. R (x) is now written as Rx in anticipation of discretizing the regularization term. The discrete form of the Laplace operator, R (δ d) = ∇ 2 (δ d), is used in this work in order to promote a spatially smooth solution.
x and v are updated at each iteration as
. Updates to x and v for the n + 1 iteration are found by solving the system of equations shown in equation (18). In equations (19) and (20), R i is the ith row of the regularization matrix, v i is the ith entry of the dual vector, and β is a small scalar term.
The direct SSIP was solved on a hexahedral mesh with 578 tri-linear elements and 972 nodes. The L-curve method was used to select the regularization hyperparameter, and a convergence criterion of the relative residual ≤ 10 -3 was used. Figure 3 shows the direct SSIP-predicted displacement field and third principal stress distribution along with the solution from a commercial finite element solver, ANSYS Workbench 14.5.
To clarify the intent of comparing the direct SSIP predictions to a commercial finite element package, the previously described forward models were used by the SSIP to reconstruct displacement field from experimental voltage-current data (i.e., the SSIP does not simulate the mechanics of the material; it only seeks to find a state of deformation that gives rise to an electrical response matching the experimental measurements). These models were all written in-house by the author in Matlab. ANSYS was used to simulate the loading conditions and obtain displacement and stress values against which SSIP predictions can be compared. In ANSYS, the elastic data presented in Figure 1 was used to define a neo-Hookean material. The marble indentors were simulated as rigid spheres subject to downward forces of 3.92 N each. The ANSYS mesh used 11,520 hexahedral elements and had five elements in the thickness direction. ANSYS was used to perform a frictional contact simulation between the marbles and the CNF/PU. In the absence of data for the coefficient of friction between CNF/PU and glass, values ranging from 0.25 to 1.25 were all considered, but it was found that solution is very insensitive to this parameter. A coefficient of friction of 1.0 was consequently used to in Figure 3.
In Figure 3, can be seen that the SSIP accurately predicts the peak displacement magnitude. Further, the direct SSIP predicts virtually zero displacement away from the indentors, which is expected both physically and mathematically due to the use of the ℓ 1 -norm on the regularization term.
From the displacement field, the left Cauchy-Green tensor was calculated as B i j = F ik F jk after projecting the gradient of the displacement field to the nodes using standard ℓ 2 -projection methods. From the strains and using the elastic properties of the CNF/PU, the Cauchy stress for a quasi-incompressible neo-Hookean elastic solid were calculated using equation ( 21) [40]. In equation ( 21), µ = E/2 (1 + ν) is the shear modulus and K = E/3 (1 -2ν) is the bulk modulus. E and ν are the elastic modulus and Poisson's ratio of the CNF/PU. J is the determinant of the deformation gradient.
From the stress tensor, the third principal stress was calculated (i.e., the left-most principal stress on a number line) since this represents the stress of compression acting on the sensor (i.e., the pressure load). Comparing the direct SSIP-predicted third principal stress and the ANSYS-predicted third principal Bottom row: The direct SSIP-predicted displacement magnitude (left) and third principal stress (right). Insets show the domain slightly rotated so that the indentation curvature can be seen. Note that displacements are multiplied by a factor of five ease of visibility. stress in Figure 3, it can be seen that the direct SSIP again does a good job-the peak stress magnitudes are very close and virtually zero stress exists away from the indentors.
In summary, this manuscript derived and discretized the direct SSIP, used the PDIPM to solve the direct SSIP, and experimentally validated this approach for pressure sensing in CNF/PU. This work was motivated by the broad potential of selfsensing materials but the limitation that electrical effects are generally not of interest. That is, it is much more useful to know the underlying condition of the material that gives rise to an observed change in conductivity. Further motivation was found in eliminating the EIT step from the SSIP. The direct SSIP-predicted displacements and stresses were compared to a commercial fi-nite element solution with very good agreement.
Based on the results presented in this manuscript, a couple of conclusions can be drawn. First, it is indeed possible to eliminate the EIT step from the SSIP. This is beneficial because it frees the SSIP from the assumptions and limitations inherent to EIT. And second, regularization terms and norms are very important to the SSIP. Even though only one combination of regularization term and norm (i.e., the discrete Laplace operator with the ℓ 1 -norm) was considered in this manuscript, these were selected with the expectation that they would favor the expected solution-a sparse and smoothly varying displacement field.
Several recommendations for future work can also be made. First, extension of these techniques to more complicated shapes (especially shapes representative of structural components) is needed. Second, mechanics-based constraints and/or regulariz-ers should be integrated into the SSIP. Currently, the inversely predicted displacement field only needs to satisfy the electrical data. Physically, however, elasticity derived constraints must also be satisfied (e.g., stress equilibrium). And third, these techniques should also be extended to material systems of more structural relevance such as carbon fiber composites. Doing so will require suitable piezoresistivity relations, which currently do not exist for highly anisotropic and/or layered materials.
Copyright © 2025 by ASME