- Email: [email protected]
Methodologies for constitutive model parameter identification for strain locking materials
Mechanics of Materials 134 (2019) 30–37
Contents lists available at ScienceDirect
Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat
Research paper
Methodologies for constitutive model parameter identification for strain locking materials
T
Nicholas Payne, Kishore Pochiraju
⁎
Department of Mechanical Engineering, Stevens Institute of Technology, Hoboken, NJ 07030, United States
ARTICLE INFO
ABSTRACT
Keywords: Strain locking materials Limited strain Optimization Parameter identification Nonlinear material characterization Natural materials Engineered materials Carotid Thoracic Arterial tissue
Strain locking materials exhibit a pronounced stiffening above an apparent critical strain. In practice, both natural materials, such as biological tissues, and human-made composites with flexible or undulating microstructures, can exhibit this type of response. This paper describes several methodologies for identifying the material constants required for a thermodynamically consistent constitutive model, that appropriately simulates the strain-locking behavior. Attempts to fit the material constants by minimizing an RMS error function between model behavior and experimental stress-strain curves reveal a non-smooth surface of the stress-based error function, which lead to difficulties in convergence to a material constant set that has a minimum error. An analysis of the constitutive model based on the relationships between the parameters of the constitutive model and physical behavior they govern lead to the construction of a smoother error function and minimization of an objective function with reduced dimensionality. This technique yielded parameters values that accurately represented the experimental data for several material systems. Two optimization methods (both gradient-based and direct) were investigated and their effectiveness in converging to the global minimum solution of the error function was compared. Selected composites and biological materials with strain-locking behaviors were analyzed, and the material constants required for the constitutive model were successfully determined.
1. Introduction Some materials exhibit a dramatic stiffening behavior as they are stretched. The stiffening can be pronounced that the material has the appearance of locking at a particular critical strain limit. This stiffening effect is due to progressive reorientation of the microstructure as stretching occurs. Specifically, stiffer constituents within the microstructure that initially have a curved, or disordered configuration, gradually straighten and align with the direction of the applied load to cause a substantial increase in stiffness. Fig. 1 shows a schematic of a typical compliance curve of the material with the strain on the vertical axis and stress on the horizontal axis. Illustrations of the microstructure realignment are shown in the inset bubbles, where the structures responsible for stiffening are initially wavy at zero stress, but gradually straighten as the material is stretched. As the stress increases, the slope of the curve approaches horizontal as the strain reaches a critical locking value. Examples of materials that have a strain locking response include both engineered and natural materials. These include some biological tissues composed of collagen fibrils, which produce strain locking behavior as the fibrils straighten and align with the applied load (Billiar and Sacks, 2000; Freed and Rajagopal, 2016). ⁎
Corresponding author. E-mail address: [email protected] (K. Pochiraju).
https://doi.org/10.1016/j.mechmat.2019.04.004 Received 29 August 2018; Received in revised form 18 February 2019 Available online 06 April 2019 0167-6636/ © 2019 Published by Elsevier Ltd.
Additionally, composite materials can be fabricated with particular architectures which produce strain locking behavior. These are composites which utilize a very compliant matrix material with embedded wires or meshes which are much stiffer than the matrix material (Ma et al., 2017). The bulk response of the composite is dominated by the highly compliant matrix material at low strains, but exhibit strain locking as the embedded inclusions realign and progressively take on more of the applied load (Ma et al., 2016). Such composite materials have numerous potential applications, particularly in the design of flexible electronic substrates. Layers of strain locking material can be used to allow flexibility at low strains, yet stiffen drastically to protect embedded electronic circuitry from experiencing larger strains which may cause damage. Additionally, wearable and implantable devices could be enabled by such composite materials which are tailored to have a similar mechanical response to human skin or other internal tissues (Wang et al., 2015). The purpose of tailoring the response is to match the behavior of the human body and allow the implanted device to function without impeding or altering the natural mechanics of the body. Developing strain locking materials for biomedical or other uses first requires a knowledge of both the mechanical environment where the device will
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju
Fig. 2. Stress-strain behavior of the material defined by α = 0.1, β = 0.001 1 ,
Fig. 1. Strain vs. stress plot of a strain locking material.
γ = 0.1 1 , ι = 0.001 Pa
be employed and accurate descriptions of the constitutive behavior of both the strain limiting composite and biological tissues. In this communication, a previously derived constitutive model for a strain locking material is presented which is constructed in such a way as to include a limiting strain directly in the mathematics of a strain energy function for the model. The model is defined by four material constants which characterize the stress/strain response of the material. Attempts to fit these constants to experimental data, by way of minimizing an error function with optimization algorithms, fail to converge on global minimum solutions due to a highly irregular error surface. Analysis of the constitutive model is conducted in this study which establishes relationships between the material constants and physically meaningful quantities. The error surface is transformed and smoothed by incorporating these relationships in the constitutive model which enable the optimization methods to converge to global minimum solutions reliably.
Pa2
Pa
. Strain asymptotically locks at about 30%.
stored strain energy within the material. (2)
W = W( ) The model establishes the form of W to be:
W (I1, I2) =
I1
1
ln(1 + I1) +
1 + 2 I2
(3)
1 tr( 2) . 2
At zero stress, the with the stress invariants: I1 = tr( ), I2 = slope of the energy function with respect to stress is zero, but as stress increases, the slope increases and attains a constant value. This feature of the strain energy function produces the strain limit. The strain function which is derived by inserting Eq. (3) into Eq. (1), which yields:
=
1
1 I+ (1 + I1)
(4)
1 + 2 I2
A corresponding plot of strain vs. stress for the material constant values of α = 0.1, β = 0.001 1 , γ = 0.1 1 , and ι = 0.001 1 2 is shown in Pa Pa Pa Fig. 2. The set of material constants, [α, β, γ, ι] control locking strain and shape of the strain energy curve. By examining the strain equation, physical significance for the material constants can be assigned. Alpha is a unitless, non-dimensional scaling factor; beta, being attached to the first stress invariant, I1 = tr( ) , is a bulk compliance constant with the unit 1 (where stress signifies force2 in whichever set of consistent
2. Constitutive model for strain locking materials Different hyperelastic models describe the behavior of materials with limited extensibility. The worm-like chain model (WLC) is one such model which is based on the relative force required to extend long chains of molecules which exist in a series of misaligned segments (hence the term worm) that align upon the application of a force (Ogden et al., 2006). This force-extension relationship defines a strain energy function, W(σ), which for a hyperelastic material is related to the strain function by = W / . The Gent model is another such hyperelastic model with strain locking that is also based on limited polymer chain extensibility (Horgan, 2015). The strain energy function and strain function are described by a stress singularity when a limit to extension is reached in the material. Originally, much of the work on constitutive modeling of inextensible hyperelastic materials was done with an eye towards representing incompressibility behavior of rubberlike materials and biological materials comprised of networks of proteins (Horgan and Saccomandi, 2003; Puglisi and Saccomandi, 2016). The concept of constructing implicit elastic strain locking constitutive models whose mathematics contain a limiting strain but are also thermodynamically consistent (Bridges and Rajagopal, 2015) was introduced by Rajagopal (2011, 2003). This model is designed to describe material systems with a limiting strain. The particular model that is used in this study comes from the work of OrtizBernardin et al. (2014) and Montero et al. (2016). A comprehensive discussion of implicit constitutive laws and strain locking models are discussed by Bulíček et al. (2014). First, beginning with the premise that strain is described as a function of stress,
= g( ) = W/
1
stress
length
units); gamma is another compliance constant with the unit 1 which stress is coupled with the stress tensor, σ; finally, iota is a squared compliance constant with the unit 1 2 attached to the second stress invariant, stress
I2 = 2 tr( 2) . These four parameters must be characterized to describe the stress-strain response of a strain locking material based on experimental data. 1
3. Mechanical stability and parameter constraints Further analysis provides insight into the constraints on the relative values which [α, β, γ, ι] may take and still yield a physically admissible constitutive response. A stiffness tensor derived for a linearized version = 0 (see Ortizof the weak form of the governing equation Bernardin et al., 2014) can be produced from:
C( ) =
(5)
where the linearized weak form is: ij
ij kl (k 1)
(k ) kl d
=
t
ui t^i d
(k 1) d ij ij
(6)
With (k) correspond to the step number in an incremental numerical solution scheme. The corresponding compliance tensor can be found by
(1)
noting that Eq. (5) can be also be written as C( ) =
moreover, there exists a scalar energy function which represents the 31
=
( ) g( )
1
with
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju
= g ( ) , so that its inverse produces the compliance tensor:
C 1( ) =
(7)
After substituting Eq. (1) to Eq. (1), the compliance tensor is found to be: ij
=
(1 + I1)2
kl
ij kl
3
(1 + 2 I2) 2
ij kl
+
1
2(1 + 2 I2) 2
(
ik jl
+
il jk )
(8) After taking the limit as the stress approaches zero, the compliance for the uniaxial case at zero stress is:
limC 1 ( ) = 0
(
)
Fig. 3. Experimental stress–strain relationships obtained from the literature (Chirita and Ionescu, 2011) for porcine carotid and thoracic artery tissue.
(9)
For thermodynamic stability, it is assumed that W must be positive definite (Malvern, 1969), therefore the compliance tensor C 1 must be positive definite. The physical significance of this stability criteria is that for any mechanical work done to the system, the strain energy in the material must increase and thus it is required that a positive stress increment result in a positive strain increment and vice versa. This condition is met when the eigenvalues of C 1 are all positive values. Here, the eigenvalues for the stress-free case are found to be: 1
=
(
);
2
=
5 4
;
3
=0
useful to fit a reduced subset of parameters initially and then fix those values and fit the remaining parameters in subsequent optimization steps (Chen and Kam, 2007; Payne and Pochiraju, 2017). The uniaxial tensile stress-strain experimental datasets of porcine carotid and thoracic (longitudinal and transverse orientation) artery tissue from Chirita and Ionescu (2011) are used as example datasets for identification of parameters of the constitutive relationship. The character of each of these data sets, which are plotted in Fig. 3, display a prominent strain locking phenomenon. Here, a root-mean-squared error (RMSE) function is selected which has the form:
(10)
First, requiring that α, β, γ, ι ≥ 0 and recognizing that λ1 in Eq. (10) is the compliance at zero stress determined in Eq. (9), the following constraint is obtained:
>
error =
(11)
1 N
N
(yexp
ymodel ) 2
1
(12)
where yexp represents the strain values in the experimental data set and ymodel is the strain value given by the constitutive model and N is the total number of data points. Fitting the constitutive model to the test data is done by employing an optimization algorithm which seeks to find the values of the parameter set [α, β, γ, ι] which results in the smallest possible value of the error function. A diagram of this process is shown in Fig. 4, where the dashed arrow indicates a data input that is supplied to the error function and the solid arrows represent the optimization algorithm, which is employed to search the parameter space for a minimum of the error function.
4. Identification of the material constants In this study, the process to fit the strain locking material model by determining the values of [α, β, γ, ι] is carried out by using numerical optimization. The optimization process finds a set of values for the material constants that minimize the value of an error function which represents the degree to which the constitutive model represents the actual response of the material measured in the experiment. This type of approach has a history of use in characterizing hyperelastic material models (Ogden et al., 2004; Twizell and Ogden, 1983) where an error function is selected which is constructed of the squared 2-norm of the difference between the stress defined by the material model and the measured stress response from the experiment. There are also other applications of optimization methods to characterize material constants based on experimental data that use measures other than stress response. For example, in the characterization of composite laminate materials, vibrational frequencies can be measured in the experiment and an error function constructed from the measured eigenfrequencies and the eigenfrequencies produced by a numerical model defined by the trial set of material parameters (Araujo et al., 1996). Additionally, strain measured during an experiment can be used as reported by Kam et al. (2009) which loaded composite beam specimens in 3-point bending configuration and constructed an error function to minimize the difference between the measured strains and the strain produced by a closed form solution for a composite beam loaded in an identical configuration. In other instances, error functions which compare experimental data to the results of finite element analyses are used within optimization algorithms. Sasso et al. (2008) performed the dynamic compression of ductile metal cylinders obtained by the split Hopkinson Bar pressure test and compared the compressive stress results of a dynamic finite element analysis, where the error function that was minimized was based on the measured stress during experiment and resulting stress of the finite element simulation. Multiple optimization steps can also be used to fit material parameters, where it is more
4.1. Initial four parameter optimization The optimization algorithm first selected to perform the search on the parameter set [α, β, γ, ι] is sequential quadratic programming (SQP) (for an introduction and overview of the algorithm, see Nocedal and Wright, 2006). Since parameter constraint given by Eq. (11) has been determined which bound the sets of parameters which are physically admissible, it is advantageous to incorporate these into the optimization to reduce the search area of the SQP algorithm and ensure physically valid results. Additionally, we require that the bounds of the search space be constrained to values that are between 0 and positive infinity in consideration of the physical meaning of the material
Fig. 4. Optimization scheme to fit the parameters of a constitutive model to experimental data. 32
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju
Fig. 5. 4-Parameter optimization results for carotid artery data. The start point used is [0.01 0.1 0.01 0.001]; (a-d) Final converged error function value vs. final parameter value. (e) Final model fit corresponding to the lowest error value. (f) The plot of the residuals for the fit. (g) Strain vs. stress plots of consecutive attempts to fit the data using the same start point [0.01 0.1 0.01 0.001].
constants which correspond to the compliance response of the material as previously discussed. A multi-start point approach is used where an array of 500 randomized points is generated each within ± 10% of a user-supplied start point. Each random point is executed separately in the SQP algorithm, and the resulting parameter values and the error value at convergence for each run are collected and compared. The results for this multi-start point optimization study for the carotid artery material are shown in Fig. 5. The non-smoothness of the error surface is evident in the widely scattered distribution of converged values for each parameter shown in Fig. 5 plots (a–d). In each plot, the parameter set corresponding to the lowest error value reported at convergence is plotted with a solid black square marker. It is not graphically demonstrated that this solution represents a global minimum
value and furthermore when observing the fit graphed along with the experimental data in the plot (e), the constitutive law charted by the blue line poorly represents the data both in the initial region and as the strain progresses to the locking value. Additionally, the locking value produced by the fit is about 0.5 which is significantly higher than what is evident in the experimental data, which is about 0.45. The plot of residuals in Fig. 5 plot (f) depicts the poor representation of the test data. Furthermore, the inaccuracy and unrepeatability of this fitting method are seen in a series of subsequent fittings of the same data. In Fig. 5 plot (g), the strain vs. stress results are plotted for four consecutive fitting attempts that are overlaid onto the experimental data for comparison. As can be seen in the wide variation between the response curves corresponding to each attempt, the optimization 33
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju
Fig. 6. 3-Parameter optimization results for carotid artery data. The start point used is [0.01, 0.1, 0.05]; lock = 0.45 . (a–d) Final converged error function value vs. the final parameter value. (e) Final model fit corresponding to the lowest error value. (f) Plot of the residuals for the fit. (g) Strain vs. stress plots of consecutive attempts to fit the data using the same start point [0.01 0.1 0.05].
algorithm using the array of randomized start points converges to a different solution each time. This behavior indicates the existence of a large number of local minima which populate the error surface and present difficulty for the optimization algorithm to converge on common solutions among even minutely differing start points. To demonstrate this behavior, consider the two start points defined by the two sets of values of [0.01, 0.1, 0.01, 0.0001] and [0.01, 0.1, 0.01, 0.000101]. In the second set, α, β, and γ are identical to the first. However, ι differs by just 1%. After running the SQP algorithm solely for each start points, two very different solutions are produced for each. The first fitted solution set is [1.593E−02, 0, 3.147E−03, 2.636E−08] while the second is [1.668E−02, 1.162E−04, 4.392E−04, 1.190E−10] respectively. These are two very different result sets with
parameter values that are orders of magnitude different in the case of ι. There is such a large degree of sensitivity to the starting point for the 4parameter optimization approach that it becomes computationally intractable to adequately sample the parameter space with initial start points to ensure that a global solution is found. 4.2. Modification of the optimization approach using ɛlock parameter relationship Due to the difficulties encountered in performing an optimization search of all four parameters of this constitutive model, a transformation of the error surface is sought to reduce the multitude of local minima. A modification to the approach is made to increase the 34
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju
Fig. 7. Pattern search optimization results for carotid artery data. The start point used is [0.01, 0.1, 0.05]. (a) The plot of the model fit. (b) The plot of the residuals for the fit.
smoothness of the error surface and accomplish a more stable model fitting process. The modification entails defining an initial value of a measurable, physically meaningful quantity and using this value definition to remove a parameter from the optimization search space and instead determine its value in a manner which is fully consistent with the physics of the constitutive model. In the following implementation, the parameter ι is removed from the search space by utilizing a parameter relationship to the locking strain, ɛlock. Consider the following, that an expression for the locking strain can be determined by the limit as the stress goes to infinity: lock
= lim ( ) =
1
Fig. 8. (a) Repeated triangular unit within the triangular “horseshoe” lattice with initially wavy sides (Jang et al., 2015). θ measures the initial angle of the arc. (b) The measured strain vs. stress response of a triangular "horseshoe" (120° arc length) lattice composite architecture and the fitted strain locking model using the 4-parameter approach. (c) The measured strain vs. stress response of a triangular "horseshoe" lattice composite architecture and the fitted strain locking model using the 3-parameter approach.
(13)
where, an additional requirement for physical admissibility is revealed from the consideration that ɛlock > 0. The following physical constraint is obtained from Eq. (13) by requiring it to be positive: 2
by noting that the residual plot (f) in the 3-parameter solution is much more focused about y = 0 with a resulting R2 measure of 0.99 as opposed to 0.95 in the 4-parameter result obtained in Fig. 5. Additionally, the final error value for the 3-parameter solution is 0.0085 as opposed to 0.028 obtained with the 4-parameter approach which is roughly a 70% reduction in the error measure. Finally, in Fig. 5 plot (g), a series of four repeated attempts at fitting using the multi-startpoint approach each with 500 randomized initial startpoints is again conducted. Now, the spread of fitting results observed in the consecutive attempts for the 4-parameter approach show in Fig. 5 plot (g) has collapsed to a single fitting solution which is reliably found in each consecutive fitting attempt. The transformation made to the error surface has increased the effectiveness of the SQP algorithm employed in this case with the array of random startpoints which were generated in a region whose bounds were manually set prior to searching. However, the error surface transformation presented in this study need not be limited to being used with this specific startpoint approach and could be used in conjunction with more sophisticated schemes for sampling the search space including adaptive Monte Carlo methodology (Kan and Timmer, 1984; Patel et al., 1989) and the OQNLP global search method (Ugray et al., 2007).
(14)
>
Additionally, Eq. (13) is rearranged to produce the following parameter relationship:
=
2 2 /(
+
lock )
2
(15)
where ɛlock is the locking strain and can is predefined from experimental data a priori of the fitting process. The derived Eq. (15) is the constraint equation which allows ι to be omitted from the search space and have a value fixed consistently with the constitutive model and observed locking limit in the experimental data. After examining the carotid artery data in Fig. 3, the locking strain is selected to be 0.45 for that particular data set. Now the multi-start point approach is repeated using a reduced parameter search space of [α, β, γ]. The results of this 3parameter version of the approach for the carotid artery data are shown in Fig. 6. The error value verses the parameter values take on a smoother, parabolic distribution (with the exception of β which converges to 0 for all executed runs) with the minimum graphically demonstrated at the point indicated with the solid black marker. Visually comparing the plots (e) of the fit in Figs. 5 and 6, the 3-parameter solution matches the test data much more closely which is also confirmed 35
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju
characterization of the material constants is conducted on the strain locking of the data corresponding to the response between 0.5 and 1.06 strain. Again, it is observed that the original 4-parameter search approach fails to produce a global minimum. However the transformed 3parameter SQP and pattern search methods both produce matching solutions ([α, β, γ, ι] = [0.34, 1.4E−4, 1.1E−4, 1.7E−9]).
Table 1 Summary of material constants obtained for various materials. Method
Carotid
SQP: 4-param SQP: 3-param Pattern Search: 3param
α
β
γ
ι
RMS error
0.00371 0.00276 0.00276
0.00231 0 1.60E−12
0.02774 0.00342 0.00342
4.95E−03 4.34E−10 4.34E−10
0.014 8.51E−03 8.51E−03
4.4. Characterization of strain locking composite materials The strain locking material model characterized in this study can also describe the response of a composite material whose architecture has been designed to produce a locking behavior as the material is progressively stretched. The data from the experiment by Jang et al. (2015), where a composite specimen composed of a compliant silicon substrate with an embedded initially wavy copper wire mesh was tested in uniaxial tension is taken and fitted with the 3parameter fitting method. The geometry of the unit cell of the wire mesh is shown in the diagram Fig. 8 plot (a). It is composed of segments of circular arcs with the shape defined by the initial angle measure θ and is sometimes termed a “horseshoe” design. As the composite is stretched, the wires that comprise each unit cell elongate and progressively become straight which produce a drastic stiffening effect. The tensile test data for this composite is plotted, and the strain locking model characterized using the 4-parameter method and shown in Fig. 8 plot (b). Then, the characterization is repeated using the proposed 3parameter method, and the results are shown plotted in Fig. 8 plot (c). The solution determined in by the former method produced an iota value of 0 which results in an approximately linear solution with a final residual error value of 0.0382. In comparison, the latter method produced a fit with a much lower residual error value of 0.0018 and follows the nonlinear character of the measured strain vs. stress response (Table 1).
Thoracic (transverse)
SQP: 4-param SQP: 3-param Pattern Search: 3param
α
β
γ
ι
RMS error
0.0063 0.0455 0.0454
1.26E−03 3.14E−05 3.13E−05
1.49E−03 1.16E−04 1.16E−04
4.03E−10 4.25E−10 4.25E−10
0.011 3.42E−03 3.42E−03
Thoracic (longitudinal)
SQP: 4-param SQP: 3-param Pattern Search: 3param
α
β
γ
ι
RMS Error
1.73E−05 0.348 0.347
5.77E−01 1.44E−04 1.44E−04
5.77E−01 1.10E−04 1.11E−04
4.83E−10 1.75E−09 1.75E−09
0.218 0.012 0.012
Strain Locking Composite (“horseshoe” lattice)
SQP: 4param SQP: 3param Pattern Search: 3-param
α
β
γ
ι
RMS error
4.548E−08
0.215
1.387
0
0.038
0.017
0
7.875E−06
2.897E−13
4.50E−03
2.344E−04
5.92E−04
5.92E−04
3.79E−13
5.15E−03
5. Concluding remarks A strain-limiting material model with four parameters was characterized using experimental data of biological tissue and a flexible composite material which exhibit strain locking response. A typical approach to fitting the material model using an optimization approach which employed sequential quadratic programming (SQP) initially failed to produce model fits that were adequately representative of the experimental data nor demonstrated to be global minimum solutions. The results from these initial attempts suggested a highly non-smooth error surface which contained a large number of local minima that made converge of the optimization algorithm difficult. To reduce this complexity in the error surface, we performed an analysis of the constitutive model that produced a relationship between the locking strain and the material constants which allowed the number of parameters included in the search space was reduced from 4 to 3. When incorporated into the optimization algorithm, the parameter relationship produced a transformation of the error surface which enabled global minimum solutions to be attained for both a gradient-based and direct optimization method.
4.3. Pattern search There are additional optimization methods which can be used to corroborate the solution obtained by the gradient-based SQP method in the previous section. Among these are optimization methods that do not rely on computing a numerical gradient to facilitate the search procedure (Hooke and Jeeves, 1961; Torczon, 1997). A pattern search method is one such algorithm that employs a grid of sample points within the parameter space on which values of the error function are computed. Based on the relative values of the error at each point, the mesh grid is translated, dilated or contracted as the algorithm searches for a minimum solution. These gradient-free methods are frequently called direct search algorithms. They often are advantageous in problems where the error surface is non-smooth, and additionally, they tend to be more effective at converging on a global minimum, bypassing local minimum solution points. However, the advantages often come at the expense of increased computational cost. In this study, optimization using the pattern search method is supplied with the same bounds and parameter constraints as the 3-parameter SQP approach outlined previously. The result plots for fitting the carotid artery data are shown in Fig. 7 and demonstrate that the same solution resulted. In the case of characterizing the constants of the longitudinal thoracic tissue data shown in Fig. 3, the experimental data show that the response of the material has an initial portion between 0.1 and 0.5 strain that has a relatively consistent low modulus. After 0.5 strain, the response begins to show the characteristic strain locking phenomenon as the other two data sets. This material would be more aptly described with a piecewise constitutive response: between 0 and 0.5 strain, the behavior can be approximated as linear, and above 0.5 strain the strain locking model is applicable. Therefore for this present study, the
Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.mechmat.2019.04.004. References Araujo, A., Soares, C.M., De Freitas, M.M., 1996. Characterization of material parameters of composite plate specimens using optimization and experimental vibrational data. Compos. Part B 27 (2), 185–191. Billiar, K.L., Sacks, M.S., 2000. Biaxial mechanical properties of the natural and glutaraldehyde treated aortic valve cusp - Part I: experimental results. J. Biomech. Eng. Trans. ASME 122 (1), 23–30.
36
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju Bridges, C., Rajagopal, K., 2015. Implicit constitutive models with a thermodynamic basis: a study of stress concentration. Zeitschrift für angewandte Mathematik und Physik 66 (1), 191–208. Bulíček, M., et al., 2014. On elastic solids with limiting small strain: modelling and analysis. EMS Surv. Math. Sci. 1 (2), 283–332. Chen, C., Kam, T., 2007. Elastic constants identification of symmetric angle-ply laminates via a two-level optimization approach. Compos. Sci. Technol. 67 (3–4), 698–706. Chirita, M., Ionescu, C., 2011. Models of Biomimetic Tissues for Vascular grafts, in on Biomimetics. InTech. Freed, A.D., Rajagopal, K.R., 2016. A promising approach for modeling biological fibers. Acta Mech. 227 (6), 1609–1619. Hooke, R., Jeeves, T.A., 1961. `` Direct Search'' solution of numerical and statistical problems. J. ACM 8 (2), 212–229. Horgan, C., Saccomandi, G., 2003. A description of arterial wall mechanics using limiting chain extensibility constitutive models. Biomech. Model. Mechanobiol. 1 (4), 251–266. Horgan, C.O., 2015. The remarkable Gent constitutive model for hyperelastic materials. Int. J. Non Linear Mech. 68, 9–16. Jang, K.-I., et al., 2015. Soft network composite materials with deterministic and bioinspired designs. Nat. Commun. 6, 6566. Kam, T., Chen, C., Yang, S., 2009. Material characterization of laminated composite materials using a three-point-bending technique. Compos. Struct. 88 (4), 624–628. Kan, A.R., Timmer, G.T., 1984. Stochastic methods for global optimization. Am. J. Math. Manag. Sci. 4 (1–2), 7–40. Ma, Y., et al., 2017. Design and application of ‘J-shaped'stress–strain behavior in stretchable electronics: a review. Lab Chip 17 (10), 1689–1704. Ma, Y.J., et al., 2016. Design of strain-limiting substrate materials for stretchable and flexible electronics. Adv. Funct. Mater. 26 (29), 5345–5351. Malvern, L.E., 1969. Introduction to the Mechanics of a Continuous Medium. Montero, S., Bustamante, R., Ortiz-Bernardin, A., 2016. A finite element analysis of some boundary value problems for a new type of constitutive relation for elastic bodies. Acta Mech. 227 (2), 601–615. Sequential Quadratic Programming, in Numerical Optimization. In: Nocedal, J., Wright,
S.J. (Eds.), Sequential Quadratic Programming, in Numerical Optimization. Springer New York, New York, NY, pp. 529–562. Ogden, R.W., Saccomandi, G., Sgura, I., 2004. Fitting hyperelastic models to experimental data. Comput. Mech. 34 (6), 484–502. Ogden, R.W., Saccomandi, G., Sgura, I., 2006. On worm-like chain models within the three-dimensional continuum mechanics framework. Proc. Math. Phys. Eng. Sci. 462 (2067), 749–768. https://royalsocietypublishing.org/doi/10.1098/rspa.2005.1592. Ortiz-Bernardin, A., Bustamante, R., Rajagopal, K.R., 2014. A numerical study of elastic bodies that are described by constitutive equations that exhibit limited strains. Int. J. Solids Struct. 51 (3–4), 875–885. Patel, N.R., Smith, R.L., Zabinsky, Z.B., 1989. Pure adaptive search in Monte Carlo optimization. Math. Program. 43 (1–3), 317–328. Payne, N., Pochiraju, K., 2017. A methodology for characterization of material constants for strain-locking materials. In: Proceedings of the American Society for Composites—Thirty-Second Technical Conference. Puglisi, G., Saccomandi, G., 2016. Multi-scale modelling of rubber-like materials and soft tissues: an appraisal. Proc. R. Soc. Math. Phys. Eng. Sci. 472 (2187), 20160060. Rajagopal, K.R., 2003. On implicit constitutive theories. Appl. Math. 48 (4), 279–319. Rajagopal, K.R., 2011. Conspectus of concepts of elasticity. Math. Mech. Solids 16 (5), 536–562. Sasso, M., Newaz, G., Amodio, D., 2008. Material characterization at high strain rate by Hopkinson bar tests and finite element optimization. Mater. Sci. Eng. 487 (1–2), 289–300. Torczon, V., 1997. On the convergence of pattern search algorithms. SIAM J. Optim. 7 (1), 1–25. Twizell, E., Ogden, R., 1983. Non-linear optimization of the material constants in Ogden's stress-deformation function for incompressinle isotropic elastic materials. J. Aust. Math. Soc. Ser. B Appl. Math. 24 (04), 424–434. Ugray, Z., et al., 2007. Scatter search and local NLP solvers: a multistart framework for global optimization. INFORMS J. Comput. 19 (3), 328–340. Wang, S., Huang, Y., Rogers, J.A., 2015. Mechanical designs for inorganic stretchable circuits in soft electronics. IEEE Trans. Compon. Packag. Manuf. Technol. 5 (9), 1201–1218.
37
Contents lists available at ScienceDirect
Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat
Research paper
Methodologies for constitutive model parameter identification for strain locking materials
T
Nicholas Payne, Kishore Pochiraju
⁎
Department of Mechanical Engineering, Stevens Institute of Technology, Hoboken, NJ 07030, United States
ARTICLE INFO
ABSTRACT
Keywords: Strain locking materials Limited strain Optimization Parameter identification Nonlinear material characterization Natural materials Engineered materials Carotid Thoracic Arterial tissue
Strain locking materials exhibit a pronounced stiffening above an apparent critical strain. In practice, both natural materials, such as biological tissues, and human-made composites with flexible or undulating microstructures, can exhibit this type of response. This paper describes several methodologies for identifying the material constants required for a thermodynamically consistent constitutive model, that appropriately simulates the strain-locking behavior. Attempts to fit the material constants by minimizing an RMS error function between model behavior and experimental stress-strain curves reveal a non-smooth surface of the stress-based error function, which lead to difficulties in convergence to a material constant set that has a minimum error. An analysis of the constitutive model based on the relationships between the parameters of the constitutive model and physical behavior they govern lead to the construction of a smoother error function and minimization of an objective function with reduced dimensionality. This technique yielded parameters values that accurately represented the experimental data for several material systems. Two optimization methods (both gradient-based and direct) were investigated and their effectiveness in converging to the global minimum solution of the error function was compared. Selected composites and biological materials with strain-locking behaviors were analyzed, and the material constants required for the constitutive model were successfully determined.
1. Introduction Some materials exhibit a dramatic stiffening behavior as they are stretched. The stiffening can be pronounced that the material has the appearance of locking at a particular critical strain limit. This stiffening effect is due to progressive reorientation of the microstructure as stretching occurs. Specifically, stiffer constituents within the microstructure that initially have a curved, or disordered configuration, gradually straighten and align with the direction of the applied load to cause a substantial increase in stiffness. Fig. 1 shows a schematic of a typical compliance curve of the material with the strain on the vertical axis and stress on the horizontal axis. Illustrations of the microstructure realignment are shown in the inset bubbles, where the structures responsible for stiffening are initially wavy at zero stress, but gradually straighten as the material is stretched. As the stress increases, the slope of the curve approaches horizontal as the strain reaches a critical locking value. Examples of materials that have a strain locking response include both engineered and natural materials. These include some biological tissues composed of collagen fibrils, which produce strain locking behavior as the fibrils straighten and align with the applied load (Billiar and Sacks, 2000; Freed and Rajagopal, 2016). ⁎
Corresponding author. E-mail address: [email protected] (K. Pochiraju).
https://doi.org/10.1016/j.mechmat.2019.04.004 Received 29 August 2018; Received in revised form 18 February 2019 Available online 06 April 2019 0167-6636/ © 2019 Published by Elsevier Ltd.
Additionally, composite materials can be fabricated with particular architectures which produce strain locking behavior. These are composites which utilize a very compliant matrix material with embedded wires or meshes which are much stiffer than the matrix material (Ma et al., 2017). The bulk response of the composite is dominated by the highly compliant matrix material at low strains, but exhibit strain locking as the embedded inclusions realign and progressively take on more of the applied load (Ma et al., 2016). Such composite materials have numerous potential applications, particularly in the design of flexible electronic substrates. Layers of strain locking material can be used to allow flexibility at low strains, yet stiffen drastically to protect embedded electronic circuitry from experiencing larger strains which may cause damage. Additionally, wearable and implantable devices could be enabled by such composite materials which are tailored to have a similar mechanical response to human skin or other internal tissues (Wang et al., 2015). The purpose of tailoring the response is to match the behavior of the human body and allow the implanted device to function without impeding or altering the natural mechanics of the body. Developing strain locking materials for biomedical or other uses first requires a knowledge of both the mechanical environment where the device will
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju
Fig. 2. Stress-strain behavior of the material defined by α = 0.1, β = 0.001 1 ,
Fig. 1. Strain vs. stress plot of a strain locking material.
γ = 0.1 1 , ι = 0.001 Pa
be employed and accurate descriptions of the constitutive behavior of both the strain limiting composite and biological tissues. In this communication, a previously derived constitutive model for a strain locking material is presented which is constructed in such a way as to include a limiting strain directly in the mathematics of a strain energy function for the model. The model is defined by four material constants which characterize the stress/strain response of the material. Attempts to fit these constants to experimental data, by way of minimizing an error function with optimization algorithms, fail to converge on global minimum solutions due to a highly irregular error surface. Analysis of the constitutive model is conducted in this study which establishes relationships between the material constants and physically meaningful quantities. The error surface is transformed and smoothed by incorporating these relationships in the constitutive model which enable the optimization methods to converge to global minimum solutions reliably.
Pa2
Pa
. Strain asymptotically locks at about 30%.
stored strain energy within the material. (2)
W = W( ) The model establishes the form of W to be:
W (I1, I2) =
I1
1
ln(1 + I1) +
1 + 2 I2
(3)
1 tr( 2) . 2
At zero stress, the with the stress invariants: I1 = tr( ), I2 = slope of the energy function with respect to stress is zero, but as stress increases, the slope increases and attains a constant value. This feature of the strain energy function produces the strain limit. The strain function which is derived by inserting Eq. (3) into Eq. (1), which yields:
=
1
1 I+ (1 + I1)
(4)
1 + 2 I2
A corresponding plot of strain vs. stress for the material constant values of α = 0.1, β = 0.001 1 , γ = 0.1 1 , and ι = 0.001 1 2 is shown in Pa Pa Pa Fig. 2. The set of material constants, [α, β, γ, ι] control locking strain and shape of the strain energy curve. By examining the strain equation, physical significance for the material constants can be assigned. Alpha is a unitless, non-dimensional scaling factor; beta, being attached to the first stress invariant, I1 = tr( ) , is a bulk compliance constant with the unit 1 (where stress signifies force2 in whichever set of consistent
2. Constitutive model for strain locking materials Different hyperelastic models describe the behavior of materials with limited extensibility. The worm-like chain model (WLC) is one such model which is based on the relative force required to extend long chains of molecules which exist in a series of misaligned segments (hence the term worm) that align upon the application of a force (Ogden et al., 2006). This force-extension relationship defines a strain energy function, W(σ), which for a hyperelastic material is related to the strain function by = W / . The Gent model is another such hyperelastic model with strain locking that is also based on limited polymer chain extensibility (Horgan, 2015). The strain energy function and strain function are described by a stress singularity when a limit to extension is reached in the material. Originally, much of the work on constitutive modeling of inextensible hyperelastic materials was done with an eye towards representing incompressibility behavior of rubberlike materials and biological materials comprised of networks of proteins (Horgan and Saccomandi, 2003; Puglisi and Saccomandi, 2016). The concept of constructing implicit elastic strain locking constitutive models whose mathematics contain a limiting strain but are also thermodynamically consistent (Bridges and Rajagopal, 2015) was introduced by Rajagopal (2011, 2003). This model is designed to describe material systems with a limiting strain. The particular model that is used in this study comes from the work of OrtizBernardin et al. (2014) and Montero et al. (2016). A comprehensive discussion of implicit constitutive laws and strain locking models are discussed by Bulíček et al. (2014). First, beginning with the premise that strain is described as a function of stress,
= g( ) = W/
1
stress
length
units); gamma is another compliance constant with the unit 1 which stress is coupled with the stress tensor, σ; finally, iota is a squared compliance constant with the unit 1 2 attached to the second stress invariant, stress
I2 = 2 tr( 2) . These four parameters must be characterized to describe the stress-strain response of a strain locking material based on experimental data. 1
3. Mechanical stability and parameter constraints Further analysis provides insight into the constraints on the relative values which [α, β, γ, ι] may take and still yield a physically admissible constitutive response. A stiffness tensor derived for a linearized version = 0 (see Ortizof the weak form of the governing equation Bernardin et al., 2014) can be produced from:
C( ) =
(5)
where the linearized weak form is: ij
ij kl (k 1)
(k ) kl d
=
t
ui t^i d
(k 1) d ij ij
(6)
With (k) correspond to the step number in an incremental numerical solution scheme. The corresponding compliance tensor can be found by
(1)
noting that Eq. (5) can be also be written as C( ) =
moreover, there exists a scalar energy function which represents the 31
=
( ) g( )
1
with
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju
= g ( ) , so that its inverse produces the compliance tensor:
C 1( ) =
(7)
After substituting Eq. (1) to Eq. (1), the compliance tensor is found to be: ij
=
(1 + I1)2
kl
ij kl
3
(1 + 2 I2) 2
ij kl
+
1
2(1 + 2 I2) 2
(
ik jl
+
il jk )
(8) After taking the limit as the stress approaches zero, the compliance for the uniaxial case at zero stress is:
limC 1 ( ) = 0
(
)
Fig. 3. Experimental stress–strain relationships obtained from the literature (Chirita and Ionescu, 2011) for porcine carotid and thoracic artery tissue.
(9)
For thermodynamic stability, it is assumed that W must be positive definite (Malvern, 1969), therefore the compliance tensor C 1 must be positive definite. The physical significance of this stability criteria is that for any mechanical work done to the system, the strain energy in the material must increase and thus it is required that a positive stress increment result in a positive strain increment and vice versa. This condition is met when the eigenvalues of C 1 are all positive values. Here, the eigenvalues for the stress-free case are found to be: 1
=
(
);
2
=
5 4
;
3
=0
useful to fit a reduced subset of parameters initially and then fix those values and fit the remaining parameters in subsequent optimization steps (Chen and Kam, 2007; Payne and Pochiraju, 2017). The uniaxial tensile stress-strain experimental datasets of porcine carotid and thoracic (longitudinal and transverse orientation) artery tissue from Chirita and Ionescu (2011) are used as example datasets for identification of parameters of the constitutive relationship. The character of each of these data sets, which are plotted in Fig. 3, display a prominent strain locking phenomenon. Here, a root-mean-squared error (RMSE) function is selected which has the form:
(10)
First, requiring that α, β, γ, ι ≥ 0 and recognizing that λ1 in Eq. (10) is the compliance at zero stress determined in Eq. (9), the following constraint is obtained:
>
error =
(11)
1 N
N
(yexp
ymodel ) 2
1
(12)
where yexp represents the strain values in the experimental data set and ymodel is the strain value given by the constitutive model and N is the total number of data points. Fitting the constitutive model to the test data is done by employing an optimization algorithm which seeks to find the values of the parameter set [α, β, γ, ι] which results in the smallest possible value of the error function. A diagram of this process is shown in Fig. 4, where the dashed arrow indicates a data input that is supplied to the error function and the solid arrows represent the optimization algorithm, which is employed to search the parameter space for a minimum of the error function.
4. Identification of the material constants In this study, the process to fit the strain locking material model by determining the values of [α, β, γ, ι] is carried out by using numerical optimization. The optimization process finds a set of values for the material constants that minimize the value of an error function which represents the degree to which the constitutive model represents the actual response of the material measured in the experiment. This type of approach has a history of use in characterizing hyperelastic material models (Ogden et al., 2004; Twizell and Ogden, 1983) where an error function is selected which is constructed of the squared 2-norm of the difference between the stress defined by the material model and the measured stress response from the experiment. There are also other applications of optimization methods to characterize material constants based on experimental data that use measures other than stress response. For example, in the characterization of composite laminate materials, vibrational frequencies can be measured in the experiment and an error function constructed from the measured eigenfrequencies and the eigenfrequencies produced by a numerical model defined by the trial set of material parameters (Araujo et al., 1996). Additionally, strain measured during an experiment can be used as reported by Kam et al. (2009) which loaded composite beam specimens in 3-point bending configuration and constructed an error function to minimize the difference between the measured strains and the strain produced by a closed form solution for a composite beam loaded in an identical configuration. In other instances, error functions which compare experimental data to the results of finite element analyses are used within optimization algorithms. Sasso et al. (2008) performed the dynamic compression of ductile metal cylinders obtained by the split Hopkinson Bar pressure test and compared the compressive stress results of a dynamic finite element analysis, where the error function that was minimized was based on the measured stress during experiment and resulting stress of the finite element simulation. Multiple optimization steps can also be used to fit material parameters, where it is more
4.1. Initial four parameter optimization The optimization algorithm first selected to perform the search on the parameter set [α, β, γ, ι] is sequential quadratic programming (SQP) (for an introduction and overview of the algorithm, see Nocedal and Wright, 2006). Since parameter constraint given by Eq. (11) has been determined which bound the sets of parameters which are physically admissible, it is advantageous to incorporate these into the optimization to reduce the search area of the SQP algorithm and ensure physically valid results. Additionally, we require that the bounds of the search space be constrained to values that are between 0 and positive infinity in consideration of the physical meaning of the material
Fig. 4. Optimization scheme to fit the parameters of a constitutive model to experimental data. 32
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju
Fig. 5. 4-Parameter optimization results for carotid artery data. The start point used is [0.01 0.1 0.01 0.001]; (a-d) Final converged error function value vs. final parameter value. (e) Final model fit corresponding to the lowest error value. (f) The plot of the residuals for the fit. (g) Strain vs. stress plots of consecutive attempts to fit the data using the same start point [0.01 0.1 0.01 0.001].
constants which correspond to the compliance response of the material as previously discussed. A multi-start point approach is used where an array of 500 randomized points is generated each within ± 10% of a user-supplied start point. Each random point is executed separately in the SQP algorithm, and the resulting parameter values and the error value at convergence for each run are collected and compared. The results for this multi-start point optimization study for the carotid artery material are shown in Fig. 5. The non-smoothness of the error surface is evident in the widely scattered distribution of converged values for each parameter shown in Fig. 5 plots (a–d). In each plot, the parameter set corresponding to the lowest error value reported at convergence is plotted with a solid black square marker. It is not graphically demonstrated that this solution represents a global minimum
value and furthermore when observing the fit graphed along with the experimental data in the plot (e), the constitutive law charted by the blue line poorly represents the data both in the initial region and as the strain progresses to the locking value. Additionally, the locking value produced by the fit is about 0.5 which is significantly higher than what is evident in the experimental data, which is about 0.45. The plot of residuals in Fig. 5 plot (f) depicts the poor representation of the test data. Furthermore, the inaccuracy and unrepeatability of this fitting method are seen in a series of subsequent fittings of the same data. In Fig. 5 plot (g), the strain vs. stress results are plotted for four consecutive fitting attempts that are overlaid onto the experimental data for comparison. As can be seen in the wide variation between the response curves corresponding to each attempt, the optimization 33
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju
Fig. 6. 3-Parameter optimization results for carotid artery data. The start point used is [0.01, 0.1, 0.05]; lock = 0.45 . (a–d) Final converged error function value vs. the final parameter value. (e) Final model fit corresponding to the lowest error value. (f) Plot of the residuals for the fit. (g) Strain vs. stress plots of consecutive attempts to fit the data using the same start point [0.01 0.1 0.05].
algorithm using the array of randomized start points converges to a different solution each time. This behavior indicates the existence of a large number of local minima which populate the error surface and present difficulty for the optimization algorithm to converge on common solutions among even minutely differing start points. To demonstrate this behavior, consider the two start points defined by the two sets of values of [0.01, 0.1, 0.01, 0.0001] and [0.01, 0.1, 0.01, 0.000101]. In the second set, α, β, and γ are identical to the first. However, ι differs by just 1%. After running the SQP algorithm solely for each start points, two very different solutions are produced for each. The first fitted solution set is [1.593E−02, 0, 3.147E−03, 2.636E−08] while the second is [1.668E−02, 1.162E−04, 4.392E−04, 1.190E−10] respectively. These are two very different result sets with
parameter values that are orders of magnitude different in the case of ι. There is such a large degree of sensitivity to the starting point for the 4parameter optimization approach that it becomes computationally intractable to adequately sample the parameter space with initial start points to ensure that a global solution is found. 4.2. Modification of the optimization approach using ɛlock parameter relationship Due to the difficulties encountered in performing an optimization search of all four parameters of this constitutive model, a transformation of the error surface is sought to reduce the multitude of local minima. A modification to the approach is made to increase the 34
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju
Fig. 7. Pattern search optimization results for carotid artery data. The start point used is [0.01, 0.1, 0.05]. (a) The plot of the model fit. (b) The plot of the residuals for the fit.
smoothness of the error surface and accomplish a more stable model fitting process. The modification entails defining an initial value of a measurable, physically meaningful quantity and using this value definition to remove a parameter from the optimization search space and instead determine its value in a manner which is fully consistent with the physics of the constitutive model. In the following implementation, the parameter ι is removed from the search space by utilizing a parameter relationship to the locking strain, ɛlock. Consider the following, that an expression for the locking strain can be determined by the limit as the stress goes to infinity: lock
= lim ( ) =
1
Fig. 8. (a) Repeated triangular unit within the triangular “horseshoe” lattice with initially wavy sides (Jang et al., 2015). θ measures the initial angle of the arc. (b) The measured strain vs. stress response of a triangular "horseshoe" (120° arc length) lattice composite architecture and the fitted strain locking model using the 4-parameter approach. (c) The measured strain vs. stress response of a triangular "horseshoe" lattice composite architecture and the fitted strain locking model using the 3-parameter approach.
(13)
where, an additional requirement for physical admissibility is revealed from the consideration that ɛlock > 0. The following physical constraint is obtained from Eq. (13) by requiring it to be positive: 2
by noting that the residual plot (f) in the 3-parameter solution is much more focused about y = 0 with a resulting R2 measure of 0.99 as opposed to 0.95 in the 4-parameter result obtained in Fig. 5. Additionally, the final error value for the 3-parameter solution is 0.0085 as opposed to 0.028 obtained with the 4-parameter approach which is roughly a 70% reduction in the error measure. Finally, in Fig. 5 plot (g), a series of four repeated attempts at fitting using the multi-startpoint approach each with 500 randomized initial startpoints is again conducted. Now, the spread of fitting results observed in the consecutive attempts for the 4-parameter approach show in Fig. 5 plot (g) has collapsed to a single fitting solution which is reliably found in each consecutive fitting attempt. The transformation made to the error surface has increased the effectiveness of the SQP algorithm employed in this case with the array of random startpoints which were generated in a region whose bounds were manually set prior to searching. However, the error surface transformation presented in this study need not be limited to being used with this specific startpoint approach and could be used in conjunction with more sophisticated schemes for sampling the search space including adaptive Monte Carlo methodology (Kan and Timmer, 1984; Patel et al., 1989) and the OQNLP global search method (Ugray et al., 2007).
(14)
>
Additionally, Eq. (13) is rearranged to produce the following parameter relationship:
=
2 2 /(
+
lock )
2
(15)
where ɛlock is the locking strain and can is predefined from experimental data a priori of the fitting process. The derived Eq. (15) is the constraint equation which allows ι to be omitted from the search space and have a value fixed consistently with the constitutive model and observed locking limit in the experimental data. After examining the carotid artery data in Fig. 3, the locking strain is selected to be 0.45 for that particular data set. Now the multi-start point approach is repeated using a reduced parameter search space of [α, β, γ]. The results of this 3parameter version of the approach for the carotid artery data are shown in Fig. 6. The error value verses the parameter values take on a smoother, parabolic distribution (with the exception of β which converges to 0 for all executed runs) with the minimum graphically demonstrated at the point indicated with the solid black marker. Visually comparing the plots (e) of the fit in Figs. 5 and 6, the 3-parameter solution matches the test data much more closely which is also confirmed 35
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju
characterization of the material constants is conducted on the strain locking of the data corresponding to the response between 0.5 and 1.06 strain. Again, it is observed that the original 4-parameter search approach fails to produce a global minimum. However the transformed 3parameter SQP and pattern search methods both produce matching solutions ([α, β, γ, ι] = [0.34, 1.4E−4, 1.1E−4, 1.7E−9]).
Table 1 Summary of material constants obtained for various materials. Method
Carotid
SQP: 4-param SQP: 3-param Pattern Search: 3param
α
β
γ
ι
RMS error
0.00371 0.00276 0.00276
0.00231 0 1.60E−12
0.02774 0.00342 0.00342
4.95E−03 4.34E−10 4.34E−10
0.014 8.51E−03 8.51E−03
4.4. Characterization of strain locking composite materials The strain locking material model characterized in this study can also describe the response of a composite material whose architecture has been designed to produce a locking behavior as the material is progressively stretched. The data from the experiment by Jang et al. (2015), where a composite specimen composed of a compliant silicon substrate with an embedded initially wavy copper wire mesh was tested in uniaxial tension is taken and fitted with the 3parameter fitting method. The geometry of the unit cell of the wire mesh is shown in the diagram Fig. 8 plot (a). It is composed of segments of circular arcs with the shape defined by the initial angle measure θ and is sometimes termed a “horseshoe” design. As the composite is stretched, the wires that comprise each unit cell elongate and progressively become straight which produce a drastic stiffening effect. The tensile test data for this composite is plotted, and the strain locking model characterized using the 4-parameter method and shown in Fig. 8 plot (b). Then, the characterization is repeated using the proposed 3parameter method, and the results are shown plotted in Fig. 8 plot (c). The solution determined in by the former method produced an iota value of 0 which results in an approximately linear solution with a final residual error value of 0.0382. In comparison, the latter method produced a fit with a much lower residual error value of 0.0018 and follows the nonlinear character of the measured strain vs. stress response (Table 1).
Thoracic (transverse)
SQP: 4-param SQP: 3-param Pattern Search: 3param
α
β
γ
ι
RMS error
0.0063 0.0455 0.0454
1.26E−03 3.14E−05 3.13E−05
1.49E−03 1.16E−04 1.16E−04
4.03E−10 4.25E−10 4.25E−10
0.011 3.42E−03 3.42E−03
Thoracic (longitudinal)
SQP: 4-param SQP: 3-param Pattern Search: 3param
α
β
γ
ι
RMS Error
1.73E−05 0.348 0.347
5.77E−01 1.44E−04 1.44E−04
5.77E−01 1.10E−04 1.11E−04
4.83E−10 1.75E−09 1.75E−09
0.218 0.012 0.012
Strain Locking Composite (“horseshoe” lattice)
SQP: 4param SQP: 3param Pattern Search: 3-param
α
β
γ
ι
RMS error
4.548E−08
0.215
1.387
0
0.038
0.017
0
7.875E−06
2.897E−13
4.50E−03
2.344E−04
5.92E−04
5.92E−04
3.79E−13
5.15E−03
5. Concluding remarks A strain-limiting material model with four parameters was characterized using experimental data of biological tissue and a flexible composite material which exhibit strain locking response. A typical approach to fitting the material model using an optimization approach which employed sequential quadratic programming (SQP) initially failed to produce model fits that were adequately representative of the experimental data nor demonstrated to be global minimum solutions. The results from these initial attempts suggested a highly non-smooth error surface which contained a large number of local minima that made converge of the optimization algorithm difficult. To reduce this complexity in the error surface, we performed an analysis of the constitutive model that produced a relationship between the locking strain and the material constants which allowed the number of parameters included in the search space was reduced from 4 to 3. When incorporated into the optimization algorithm, the parameter relationship produced a transformation of the error surface which enabled global minimum solutions to be attained for both a gradient-based and direct optimization method.
4.3. Pattern search There are additional optimization methods which can be used to corroborate the solution obtained by the gradient-based SQP method in the previous section. Among these are optimization methods that do not rely on computing a numerical gradient to facilitate the search procedure (Hooke and Jeeves, 1961; Torczon, 1997). A pattern search method is one such algorithm that employs a grid of sample points within the parameter space on which values of the error function are computed. Based on the relative values of the error at each point, the mesh grid is translated, dilated or contracted as the algorithm searches for a minimum solution. These gradient-free methods are frequently called direct search algorithms. They often are advantageous in problems where the error surface is non-smooth, and additionally, they tend to be more effective at converging on a global minimum, bypassing local minimum solution points. However, the advantages often come at the expense of increased computational cost. In this study, optimization using the pattern search method is supplied with the same bounds and parameter constraints as the 3-parameter SQP approach outlined previously. The result plots for fitting the carotid artery data are shown in Fig. 7 and demonstrate that the same solution resulted. In the case of characterizing the constants of the longitudinal thoracic tissue data shown in Fig. 3, the experimental data show that the response of the material has an initial portion between 0.1 and 0.5 strain that has a relatively consistent low modulus. After 0.5 strain, the response begins to show the characteristic strain locking phenomenon as the other two data sets. This material would be more aptly described with a piecewise constitutive response: between 0 and 0.5 strain, the behavior can be approximated as linear, and above 0.5 strain the strain locking model is applicable. Therefore for this present study, the
Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.mechmat.2019.04.004. References Araujo, A., Soares, C.M., De Freitas, M.M., 1996. Characterization of material parameters of composite plate specimens using optimization and experimental vibrational data. Compos. Part B 27 (2), 185–191. Billiar, K.L., Sacks, M.S., 2000. Biaxial mechanical properties of the natural and glutaraldehyde treated aortic valve cusp - Part I: experimental results. J. Biomech. Eng. Trans. ASME 122 (1), 23–30.
36
Mechanics of Materials 134 (2019) 30–37
N. Payne and K. Pochiraju Bridges, C., Rajagopal, K., 2015. Implicit constitutive models with a thermodynamic basis: a study of stress concentration. Zeitschrift für angewandte Mathematik und Physik 66 (1), 191–208. Bulíček, M., et al., 2014. On elastic solids with limiting small strain: modelling and analysis. EMS Surv. Math. Sci. 1 (2), 283–332. Chen, C., Kam, T., 2007. Elastic constants identification of symmetric angle-ply laminates via a two-level optimization approach. Compos. Sci. Technol. 67 (3–4), 698–706. Chirita, M., Ionescu, C., 2011. Models of Biomimetic Tissues for Vascular grafts, in on Biomimetics. InTech. Freed, A.D., Rajagopal, K.R., 2016. A promising approach for modeling biological fibers. Acta Mech. 227 (6), 1609–1619. Hooke, R., Jeeves, T.A., 1961. `` Direct Search'' solution of numerical and statistical problems. J. ACM 8 (2), 212–229. Horgan, C., Saccomandi, G., 2003. A description of arterial wall mechanics using limiting chain extensibility constitutive models. Biomech. Model. Mechanobiol. 1 (4), 251–266. Horgan, C.O., 2015. The remarkable Gent constitutive model for hyperelastic materials. Int. J. Non Linear Mech. 68, 9–16. Jang, K.-I., et al., 2015. Soft network composite materials with deterministic and bioinspired designs. Nat. Commun. 6, 6566. Kam, T., Chen, C., Yang, S., 2009. Material characterization of laminated composite materials using a three-point-bending technique. Compos. Struct. 88 (4), 624–628. Kan, A.R., Timmer, G.T., 1984. Stochastic methods for global optimization. Am. J. Math. Manag. Sci. 4 (1–2), 7–40. Ma, Y., et al., 2017. Design and application of ‘J-shaped'stress–strain behavior in stretchable electronics: a review. Lab Chip 17 (10), 1689–1704. Ma, Y.J., et al., 2016. Design of strain-limiting substrate materials for stretchable and flexible electronics. Adv. Funct. Mater. 26 (29), 5345–5351. Malvern, L.E., 1969. Introduction to the Mechanics of a Continuous Medium. Montero, S., Bustamante, R., Ortiz-Bernardin, A., 2016. A finite element analysis of some boundary value problems for a new type of constitutive relation for elastic bodies. Acta Mech. 227 (2), 601–615. Sequential Quadratic Programming, in Numerical Optimization. In: Nocedal, J., Wright,
S.J. (Eds.), Sequential Quadratic Programming, in Numerical Optimization. Springer New York, New York, NY, pp. 529–562. Ogden, R.W., Saccomandi, G., Sgura, I., 2004. Fitting hyperelastic models to experimental data. Comput. Mech. 34 (6), 484–502. Ogden, R.W., Saccomandi, G., Sgura, I., 2006. On worm-like chain models within the three-dimensional continuum mechanics framework. Proc. Math. Phys. Eng. Sci. 462 (2067), 749–768. https://royalsocietypublishing.org/doi/10.1098/rspa.2005.1592. Ortiz-Bernardin, A., Bustamante, R., Rajagopal, K.R., 2014. A numerical study of elastic bodies that are described by constitutive equations that exhibit limited strains. Int. J. Solids Struct. 51 (3–4), 875–885. Patel, N.R., Smith, R.L., Zabinsky, Z.B., 1989. Pure adaptive search in Monte Carlo optimization. Math. Program. 43 (1–3), 317–328. Payne, N., Pochiraju, K., 2017. A methodology for characterization of material constants for strain-locking materials. In: Proceedings of the American Society for Composites—Thirty-Second Technical Conference. Puglisi, G., Saccomandi, G., 2016. Multi-scale modelling of rubber-like materials and soft tissues: an appraisal. Proc. R. Soc. Math. Phys. Eng. Sci. 472 (2187), 20160060. Rajagopal, K.R., 2003. On implicit constitutive theories. Appl. Math. 48 (4), 279–319. Rajagopal, K.R., 2011. Conspectus of concepts of elasticity. Math. Mech. Solids 16 (5), 536–562. Sasso, M., Newaz, G., Amodio, D., 2008. Material characterization at high strain rate by Hopkinson bar tests and finite element optimization. Mater. Sci. Eng. 487 (1–2), 289–300. Torczon, V., 1997. On the convergence of pattern search algorithms. SIAM J. Optim. 7 (1), 1–25. Twizell, E., Ogden, R., 1983. Non-linear optimization of the material constants in Ogden's stress-deformation function for incompressinle isotropic elastic materials. J. Aust. Math. Soc. Ser. B Appl. Math. 24 (04), 424–434. Ugray, Z., et al., 2007. Scatter search and local NLP solvers: a multistart framework for global optimization. INFORMS J. Comput. 19 (3), 328–340. Wang, S., Huang, Y., Rogers, J.A., 2015. Mechanical designs for inorganic stretchable circuits in soft electronics. IEEE Trans. Compon. Packag. Manuf. Technol. 5 (9), 1201–1218.
37