Topology optimization of self-sensing nanocomposite structures with designed boundary conditions

Controlling volume fractions of nanoparticles in a matrix can have a substantial influence on composite performance. This paper presents a topology optimization algorithm that designs nanocomposite structures for objectives pertaining to stiffness and strain sensing. Local effective properties are obtained by controlling local volume fractions of carbon nanotubes (CNTs) in an epoxy matrix, which are assumed to be well dispersed and randomly oriented. The method is applied to the optimization of a plate with a hole structure. Several different allowable CNT volume fraction constraints are examined, and the results show a tradeoff in preferred CNT distributions for the two objectives. It is hypothesized that the electrode location plays an important role in the strain sensing performance, and a surrogate model is developed to incorporate the electrode boundary as a set of additional design variables. It is shown that optimizing the topology and boundary electrode location together leads to further improvements in resistance change.


Introduction
Recent developments in advanced materials have led to the emergence of multifunctional structures that 'combine the functional capabilities of one or more subsystems with that of the load bearing structure' [1]. One of these capabilities is self-sensing, in which a structure is able to directly collect information about its operating environment and relay that information to pilots, testing engineers, and maintenance engineers [2,3].
The inability to embed a traditional sensor (such as a strain gauge) in the structure is a siginificant limitation for composites, in which cross-sectional or interlaminate failures may not be observable at the surface [4]. This motivates the investigation of multifunctional structures in which the sensing material is dispersed throughout the structure. Of course, this sensing structure must still perform its role as a load carrying member. Of the candidate materials for use in creating self-sensing structures, carbon nanotubes (CNTs) are the subject of much attention [5]. CNT based composite strain sensors have been shown to have higher sensitivities than classic strain gauges at the macroscale [6] and exhibit strain sensing through several mechanisms [7].
Paralleling the rise in multifunctional materials is the blooming field of additive manufacturing. Specifically, composite additive manufacturing continues to be a hot topic in labs, where 'Additive manufacturing holds strong potential for the formation of a new class of multifunctional nanocomposites through embedding of nanomaterials [8].' It is even possible to additively manufacture CNT/polymer composites with finely tailored microstructures using liquid deposition [9].
Recently, much has been done to apply topology optimization to the design of multifunctional materials. Maute et al [10] used level set topology optimization to design a set of printable shape memory polymer (SMP)-elastic matrix composites to match a specified deformed shape once actuated. This two-material system was able to closely match a deformed shape once actuated, and showed that there is indeed benefit to combining advanced manufacturing, multifunctional materials, and topology optimization. Pertaining more specifically to sensing structures in topology optimization, Rubio [11] investigated topology optimization of a piezoresistive patch in a compliant mechanism in which orientation of a monolithic Wheatstone bridge was optimized in addition to the topology of the compliant structure. Gusti, Mello, and Silva [12,13] optimized the topology of a piezoresistive membrane that was stretched over a structure to maximize the sensing capability and the stiffness. They were able to show over 150% increases in measured potential difference due to the piezoresistive effect.
However, a limitation of most published work is the focus single or two-material systems. Given the advances in additive manufacturing, it is possible to envision a finely controlled structural system with a locally varying, or graded microstructure. This paper presents an algorithm capable of designing such a structure. This is done via an application of topology optimization, in which the design variables are set to be the local CNT volume fractions of a CNT/Epoxy composite structure and the objectives are measures of structural stiffness and electrical resistance change due to strain (strain-sensing).
A plate with a hole, loaded in tension in the vertical direction, is selected as the test structure. Constraints are imposed on both the local and global CNT volume fractions as representations of manufacturing and cost constraints, respectively. Micromechanics models are used to obtain element effective properties and are functions of the local CNT volume fraction, and finite element analysis uses these local effective properties to solve for the global objectives. Sensitivities of the objectives are analytically derived and used to fuel a gradient-based, Sequential Quadratic Programming (SQP) optimization scheme.
It has also been seen that the simultaneous optimization of sensor and structure can highly depend on the selection of the electrode location [14]. While a structural loading environment is often not at the discretion of the engineer to prescribe, it is much more likely that one may choose where to locate the electrodes on a structure, and in doing so may improve sensing or decouple the trade-off between stiffness and sensing.
The test structure is first optimized with a set of fixed electrodes. Then a surrogate model is developed to incorporate the discrete electrode variables (start and end nodes within a FEA mesh) within the continuous optimization scheme.
Section 2 details the problem statement, design space, and solution algorithm. Section 3 follows with the relevant micromechanics, finite element equations, and sensitivities of the objectives. Section 4 presents results for the optimization of the fixed-electrode structure, commenting on differences in performance for various volume fraction constraints as well as interpreting what makes designs optimal in one or multiple objectives. Section 5 introduces the boundary condition surrogate model and relevant equations and sensitivities. Section 6 then solves the combined topology and boundary condition optimization problem and compares the results to the fixedelectrode case. Conclusions are presented in section 7.

Problem statement
Topology optimization seeks to design a structure by first discretizing the design space and then driving the local material volume fractions in each element of that space to their optimal values. The general problem is formulated as follows: v Here the set of design variables are designated as the vector v, and may correspond to a 'relative density' of material or a phase volume fraction. The relative density in each element is denoted by v e . The objective function F may be multiobjective, and be formed from single objective functions f 1 and f 2 . For this , the resistance change due to strain, and v f U 2 = ( ), the strain energy. In pseudo-density methods such as SIMP or RAMP [15] the design variables are used within one or more material interpolation schemes, which govern the effective material properties of the corresponding element. It is the aim of the optimization to find the set of v that minimize the objectives. Classically these methods include some penalty term that drives v e to either 0 or 1, representing either material or void. For more detailed reviews of topology optimization the reader is encouraged to examine [15] and [16].
Rather than considering a single material property and driving this property to 1 (on) or 0 (off) via a SIMP-like method, one may instead consider a micromechanics model, such as a rule of mixtures [17,18], inverse rule of mixtures [17,19], or a method that makes use of the Eshelby solution [20], such as the Mori-Tanaka method [21]. These models relate effective material properties to the volume fraction of an inhomogeneity in a matrix in a continuous manner. In the case of this paper, that inhomogeneity is the volume fraction of randomly oriented and well-dispersed CNT in each element. As no artificial penalization is added, what results is no longer an 'on' or 'off' design, but rather a distributed system of CNT-epoxy nanocomposite in which each element may have a different material composition, and different effective properties.
The design space for the plate with a hole is introduced in figure 1, making note of symmetry boundary conditions that are used to reduce the design space to a single quadrant. The left and bottom edge contain this symmetry condition, the top edge is loaded with a uniform tensile load, and the blue and red lines mark the locations of the electrodes. Constraints are placed on both the amount of CNT available to a single element and on the total amount of CNT in the cross-section; v p and V p , respectively. It is desirable that the structure have some measure of stiffness so that it can perform its structural application. Strain energy is chosen as an objective to capture the stiffness. It is also necessary that a measure of the sensing signal be maximized. This is acheived via maximization of Resistance change in the presence of strain [22,23]. These objectives will be shown to be competing for a limited amount of CNT, with the stiffness optimization wanting to place material in locations that may be disadvantageous for sensing, and vice versa. The problem is solved using an epsilon-constraint optimization, in which the strain energy objective is rewritten as a constraint [24][25][26], leading to the optimization problem where v e is the CNT volume fraction of the eth element.
v v is the resistance change between the strained and unstrained cross-section ΔR, normalized by the unstrained resistance, R 0 . U is the strain energy, and U * is a prescribed strain energy constraint. By changing U * one can trade relative importance of stiffness versus sensing in the design. However, care must be taken in the selection of U * to ensure the constraint is feasible. This problem is solved using SQP optimization within Matlab's fmincon optimization suite [27]. Default constraint, objective function decrease, and optimality tolerances are used unless otherwise specified.

Micromechanics
Micromechanics laws relate the design variables to local Young's Modulus, resistivity, and piezoresistive constant in a given element. It is assumed that within an element the CNT are well dispersed and randomly oriented, giving effectively isotropic properties. . This behavior is seen to be nonlinear even at low volume fractions, requiring use of an inverse rule of mixtures model [19]. Effective resistivity is isotropic and given as Piezoresistivity is a property that dictates how changes in strain influence resistivity. A piezoresistive constant, sometimes called a normalized gage factor, can be used to measure this property. The piezoresistive constant is denoted as the variable g, and the local effective piezoresistive constant of the eth element is g e . Depending on the percolation threshold of a given CNT-Epoxy composite, the piezoresistive behavior can exhibit an almost discrete on/off behavior [29]. Below the percolation threshold the piezoresistivity is small, and at the percolation threshold the piezoresistivity is maximized. Continuing to add CNT beyond percolation CNT will reduce the piezoresistive constant [23,30]. The percolation threshold of CNT-Epoxy composites can be as low as 0.0025% CNT volume fraction [31] but it is most common that this threshold is between 1.5% and 4.5% [29,32], depending on type and processing method of the CNT. 2% volume fraction was chosen for the percolation threshold. The effective piezoresistive constant is small before 2%, peaks at percolation, and decreases for larger volume fractions. A curve fit model is used to approximate this behavior, and is modeled after figure 8 in [33].
The sensitivity of the piezoresistive constant to the CNT volume fraction is dg dv In equations (7) and (8), the constants A 1 -A 3 , B 1 -B 2 , and C 1 -C 3 are selected to ensure that the curve is continuous and has a continuous first derivative. These parameters may be altered to tune the piezoresistive model to fit a specific manufacturing process and/or available experimental data. Table 1 shows the values identified for these constants in this paper.

Material properties.
Material properties for CNT and Epon 862 are presented in table 2 [32,34]. It should be noted that the Poisson's ratio of the nanocomposite was assumed to be a constant ν=0.3. Effective Poisson's ratios of CNT-Epon composites were modeled using a Mori-Tanaka method in [35,36], where it was found that for aligned CNT the composite effective properties were ν 12 =0.377, ν 23 = ν 13 =0.263. For randomly oriented nanotubes it can be assumed that these values may be averaged, resulting in an effective Poisson's ratio of 0.3. The micromechanics equations are plotted against CNT volume fraction in figure 2.

Strain energy
The structure's strain energy under prescribed loading is chosen as a measure of stiffness. A finite element model using 4 node bi-linear quadrilateral elements is used to solve both the mechanical and the electrostatic problems. For the mechanics, from which strain energy is calculated, each node has two degrees of freedom, u x and u y . The element equilibrium equation for the 2D mechanics is e e e e = ( ) ( ) where K e is the element stiffness matrix, v e the element CNT volume fraction, u e the element displacement vector, and f e the element load vector. The CNT volume fractions affect the stiffness matrix by modifying the constitutive matrix, D e .
where constants E mod ) . The element stiffness matrix can be written as where B k is the mechanical strain-displacement matrix and J e | | is the determinant of the element Jacobian matrix. K e is computed for each element and assembled into a global stiffness matrix, K. The symmetry boundary conditions specify which u x and u y are set to 0, and the problem is solved for the remaining displacements. The strain energy is then computed as The sensitivity of stiffness-based topology optimization problems is well studied, and the sensitivity of the strain energy is shown via a self-adjoint solution [15, ?] where the sensitivity of the stiffness matrix is dependent only on the sensitivity of the constitutive matrix.
and as the element Young's modulus has been factored out of D e the sensitivity may be easily computed.

Resistance change due to strain
Gage factor is defined as the change in resistance between the strain and unstrained cross section divided by the unstrained resistance and the strain, and is a standard measure of sensing.
Maximizing the resistance change between the strained and unstrained structure leads to an increase in signal-to-noise ratio in strain sensing. The resistance change maximization problem is formulated based on figure 1. A set of electrodes, denoted by the red and blue bars, are located on the boundary of the structure and are used to prescribe a voltage difference. Two solutions of an electrostatics finite element problem are required to obtain the resistance change, one for the base unstrained structure, and one once strain has been applied. The finite element solution is used to obtain electrical currents through the boundary electrode in both cases, which may be related to the relevant resistances through Ohm's law.
3.3.1. Electric current. The electrostatics continuity equation states that the divergence of the current density (Y) is 0.
Current density is related to electric conductivity (s) and the electric field (E) via Ohm's law as The electric field is the negative of the gradient of the potential. Substituting this into equation (16) gives In the 2D case the electric potential varies in the z and the y directions, z y , f f = ( ). Conductivity may also change in both directions, z y , s s = ( ). Rewriting the equation gives z z y z y z y z y z y y , , , 19 The governing equations are discretized via the finite element method, resulting in the algebraic equations or, for a given element where C e is the element electrostatic 'stiffness' matrix, e f is the element electric potential vector, and f e is the element current vector. The electrostatic version of the stiffness matrix depends on the conductivity matrix s.
e e T e e ò ò s where B is the gradient matrix, J e | | is the determinant of the element Jacobian, and e s is the element conductivity, which will vary between the strained and unstrained problems.
Equation (20) is divided into submatrices based on which degrees of freedom are constrained. The subscript u denotes degrees of freedom which are unspecified, but on the boundary. The subscript s indicates these degrees of freedom are part of the boundary condition, and have their electric potential specified. This represents specifying the placement of electrodes on the structure. Finally, the subscript i indicated degrees of freedom on the interior of the structure. Here the symbol II is used to represent the identity matrix.
Total current, I bc , is measured as the summation of the nodal currents across a boundary electrode. The vector q is created to aid in the summation. q has a value of 1 for degrees of freedom on the boundary electrode to be summed over, and is 0 for the degrees of freedom on the other boundary electrode.
For the unstrained calculation an uncoupled adjoint method is used to obtain the sensitivity of the current. For this is is convenient to rearrange equation (24).
where the sensitivities of the resistivity are provided by the micromechanics equations. Adding the piezoresistive term alters the conductivity matrix, as a piezoresistive term is added to the resistivity as g ii ¢ and 0 r ¢ indicate derivatives of the micromechanics equations with respect to the local volume fraction. In literature [7,35] the shear terms of the CNT-polymer composite piezoresistivity were seen to be small, and thus for the representative model used herein it is assumed that g 11 =g 12 =g and g 66 =0, where g is provided via equation (7).
General forms of coupled adjoint sensitivities are adapted for the sensitivity of the coupled piezoresistive problem [37]. The equilibrium equations are labeled The objective functions are F 1 =U for strain energy and F 2 =I bc for the strained current out of the boundary electrode. The strained current is considered here as it is the term in the resistance change that includes the coupling. In matrix form The total derivative is given as where the subscript, u and, f indicate derivatives of the residuals and objectives with respect to that state variable. If only the strained current sensitivity is required (as the strain energy sensitivity has been solved via self-adjoint), this reduces to Similarly, partitioning equation (44) to only consider the strained current objective results in an updated version of equation (31).
The adjoint variables and the sensitivities of the stiffness and conductivity matrices to volume fraction changes must now be obtained. The sensitivities of the residuals with respect to the states are The sensitivities of the objectives with respect to the states are and the adjoint equation is reposed as As p q C T = , and q is just a selection vector of 1ʼs and 0ʼs, the only remaining term to solve is R u 2, in equation (50).
This may be computed on an element-wise basis and then assembled into a global sensitivity matrix, similar to the finite element assembly of the stiffness and conductivity matrices. r is the local unstrained resistivity, g ē is the local piezoresistive matrix, B ke is the element strain-displacement matrix, B e is the element electrostatic-gradient matrix (the electrostatic analog to the strain-displacement matrix), e ŝ is the element conductivity matrix, and J e | | is the determinant of the element Jacobian.

Optimal topologies with a fixed electrode
The coupled optimization is performed using the epsilonconstraint method introduced in equation (2). Single objective optimization with a uniform CNT distribution is used to obtain utopia points that inform the bounds on the strain energy epsilon-constraint. Global volume fractions of 2% and 5% are considered, as well as a case without a global volume fraction constraint. Pareto Fronts for the plate with a hole are plotted in figure 3. The top left subplot compares Pareto Fronts across all three volume fraction constraint cases, and the remaining plots isolate a single constraint case.
The vertical axis in figure 3 marks the stiffness performance and the horizontal axis marks the sensing performance, with higher values of both being preferred. Figure 4 plots the optimal stiffness and sensing values against the volume fraction constraint. Note that because the side constraints on each element restrict the local volume fraction to be less than or equal to 10%, the unconstrained case can have a max global volume fraction of 10%.
The optimal stiffness is dominated by the volume fraction constraint, with the unconstrained case being 40% more stiff than the 5% constrained case, which was itself twice as stiff as the 2% constrained case. The topologies that perform the best in stiffness are plotted for each constraint level in figure 5 alongside their volumetric strain fields.
As adding more CNT will always increase the local Young's modulus, the unconstrained optimum for stiffness  maximization is a topology with maximum CNT in each element. Once the volume fraction constraint is activated and begins to restrict CNT usage, the stiffness is maximized by placing higher volume fractions of CNT near the right edge of the hole, minimizing the stress concentration. Another common feature in both the 2% and 5% topologies is a stiffening arc leading up and around the net-negative strain region. Away from the hole the strains are relatively uniform, and the optimizer has little preference as to volume fraction or specific topology in this region.
When comparing sensing performance the 5% constrained and unconstrained design both dominate the 2% constrained designs, but there is little to no increase in the sensing objective when relaxing the constraint beyond 5%. The optimal sensing topologies for each constraint level are plotted in figure 6 alongside their local resistivity changes due to strain.
Sensing is optimized by placing highly piezoresistivenear 2% CNT volume fraction-material near the highly strained right edge of the hole. A conductive path connecting this area to the electrode on the right vertical edge is also common across all of the designs. A region of low CNT volume fraction material in the center of the design, common across all three constraint levels, seems to be a preferred feature of the topology. The optimizer is manipulating the strain field to concentrate the load in the highly piezoresistive region, resulting in a higher sensing performance. Of course, carving out a large piece of the design and dumping more load into a relatively compliant section of the structure may not make for a good structure, but it makes for the best sensor. Figure 7 characterizes the transition from optimal structure to optimal sensor. This figure shows four of the Pareto Optimal topologies from the 5% constraint case along with their respective stiffness and sensing performance. As the stiffness epsilon-constraint is relaxed, the region of low CNT in the center of the design grows, trading stiffness for improved sensing.

A surrogate model for the designed electrode
In order to obtain measurable resistance change due to strain there must be a path of conductive CNT linking the electrodes. The location of the electrodes dictates where this path must form, and so it is of interest to allow the optimizer to tailor the electrode in addition to the CNT distribution. It has been seen that this can serve to both increase sensing performance and partially decouple stiffness and sensing objectives [14]. This presents unique challenges in that the prescribed electrodes exist as discrete degrees of freedom in a finite element mesh. Including the electrode placement within the optimization necessitates either the use of a mixed-variable optimizer or a way of converting the discrete electrode variables into pseudo-continuous variables.
Surrogate models, covered in depth in [38], use curve fitting and statistics to interpolate the behavior of a function between discrete evaluation points. Unless remeshing is performed at every optimization iteration, the boundary nodes on a mesh are fixed points and the optimal location of an electrode may well lie between nodes on the mesh. A surrogate model allows for interpolation of the resistance change performance for electrodes that end between nodes. A quadratic response surface model is developed for this purpose. = [ ] where b 1 and b 2 mark the starting and ending location of the variable electrode. The electrode along the circular edge was chosen as the variable electrode, as this area has been seen to be important in the results presented in the previous section. All nodes on the specified edge that fall between b 1 and b 2 are considered part of the electrode. The performance and sensitivity of the surrogate model is formulated as follows.
First, the continuous electrode variables are rounded to the nearest discrete value, b b round = ( ). The remainder, r b b = -¯, will also be used. As b is discrete, it is used within the established analysis and sensitivity to compute the resistance change at the rounded electrode location.
The resistance change objective at the rounded electrode values is

Results for simultaneous optimization of topology and boundary electrode
The Pareto Fronts for designed electrodes and topology are plotted alongside the optima with only the designed topology in figure 8. The 'x' indicates an optima obtained with the designed electrode, the 'o' is reprinted from the fixed electrode results in section 4. As the stiffness objective is independent of electrode placement, the stiffness-optimal designs show little to no improvement between fixed and designed electrode cases. As the strain energy epsilon-constraint is relaxed the design that includes the variable electrode dominates its fixed electrode counterpart. Table 3 compares sensing performance for select fixed and variable electrode designs. Each comparative case shows the sensing performance for a matching stiffness performance. Optimizing the electrode offers a significant sensing across all constraint levels, at least a 1.46 times increase in sensing for the values shown. This increase is more significant, at least a 1.81 times increase for 5% and a 1.88 times increase for the unconstrained case, as the volume fraction constraint is relaxed. Figure 9 plots two topologies with the same stiffness for a constraint volume fraction of 5%. The first is optimized with    a fixed electrode, and the second is optimized with a variable inner electrode. The purple bar indicates the electrode location along the inner edge for both designs. The topology with the designed electrode places stiff material (red elements) around the edge of the hole, next to the piezoresistive material. This is a feature not shared by the topology with the fixed electrode, and allows for a larger cutout in the center-right of the design while still satisfying the strain energy epsilon-constraint. Red elements are both stiff and conductive, and as such red elements attached to the inner electrode form a path of least resistance that would bypass the sensing elements in the bottom right of the hole section. With a full length electrode, there is no way to bring this stiff material down to the circular edge without sacrificing sensing performance. Another advantage of the designed electrode is that is is concentrated only at the region of highest sensing. This electrode design is consistent across the best sensing topologies for all volume fraction constraints, which are shown in figure 10.  . Select Pareto optimal topologies for the 5% constrained plate w/hole structure and the designed electrode. Top left to bottom right transitions from optimal structure to optimal sensor. In all cases, once the stiffness requirement is sufficiently reduced the designed electrode optima dominate the fixed electrode optima. This stiffness threshold corresponds to when the optimizer no longer needs stiff, poor sensing material to occupy the location of the stress concentration around the hole. For designs with intermediate stiffness requirements shifting the electrode allows for unique topologies that satisfy stiffness without subverting the sensing-optimal conductive path. Even when the stiffness constraint is completely relaxed the optimized electrode works with the topology to force the current to flow through the best sensing regions of the design.
The development of the topology with a designed electrode as the stiffness epsilon-constraint is relaxed is shown in figure 11. The development of the optimal sensing features mentioned in the preceding paragraphs are seen to develop as the stiffness requirement decreases.

Conclusions
This paper presented a method for optimal distribution of a limited amount of CNTs within an epoxy matrix to provide Pareto-Optimal designs for stiffness and strain sensing objectives. Analytic analysis and sensitivity equations based on micromechanics and FEA were developed and used within an SQP optimization routine. Designs were first optimized with a fixed electrode, and then a surrogate model was used to incorporate the electrode into the design. Main results of this work are presented as: 1. Stiffness is maximized by placing high volume % CNT elements around the stress concentrations. 2. Sensing is maximized by placing highly piezoresistive elements around the stress concentrations, forming conductive paths between electrodes, and manipulating the load path to concentrate loads in the best sensing region. 3. Stiffness monotonically increases with available CNT, sensing is less depending on material available after a threshold, around 5% for the cases presented here. 4. Adding the electrode to the optimization allows for tailored conductive paths, increasing sensing for all volume fraction cases once the stiffness requirement is low enough to allow location of piezoresistive material around the designed electrode.