Thermodynamic compatibility conditions of a new class of hysteretic materials

The thermodynamic compatibility defined by the Drucker postulate applied to a phenomenological hysteretic material, belonging to a recently formulated class, is hereby investigated. Such a constitutive model is defined by means of a set of algebraic functions so that it does not require any iterative procedure to compute the response and its tangent operator. In this sense, the model is particularly feasible for dynamic analysis of structures. Moreover, its peculiar formulation permits the computation of thermodynamic compatibility conditions in closed form. It will be shown that, in general, the fulfillment of the Drucker postulate for arbitrary displacement ranges requires strong limitations of the constitutive parameters. Nevertheless, it is possible to determine a displacement compatibility range for arbitrary sets of parameters so that the Drucker postulate is fulfilled as long as the displacement amplitude does not exceed the computed threshold. Numerical applications are provided to test the computed compatibility conditions.


Introduction
Uniaxial constitutive models describing hysteretic load-response behaviors play a significant role in several engineering applications. In general, such kinds of response depend on the past history of the applied loads, while the first derivative of the action can influence (for rate-dependent) or not (for rate-independent) the current response value [6]. Within this framework, hysteretic models are involved in several applications of structural engineering such as ductile [15] and fragile materials [2] as well as seismic devices [5,12,13,22].
An important aspect for many applications consists in the fact that it is desirable to adopt models defined by means of smooth relationships so that the response and its first derivative are continuous, an essential aspect in a large variety of applications [1,4,9,16].
To this end, several constitutive models have been developed by adopting analytical functions such as the famous material proposed by Bouc [3] and later extended by Wen in [24] and [23] whose popularity is due to its capability of reproducing complex and continuous hysteretic loop shapes by means of a single function and a limited number of parameters.
The Bouc-Wen material belongs to a wider class of differential constitutive models whose formulations are based on analytical relationships involving both the response and its derivatives. Thus, their main drawback consists in the necessity of computing the actual response by solving a nonlinear differential equation. This Communicated by Andreas Öchsner.
Research funded by University of Naples Federico II and Compagnia di San Paolo -Progect MuRA, CUP E69C21000250005.

S. Sessa (B)
Università degli Studi di Napoli Federico II, Via Claudio 21, 80125 Napoli, Italy E-mail: salvatore.sessa2@unina.it is usually done by means of computationally demanding iterative numerical approaches, being closed-form solutions available for a few structural problems.
In order to overcome such an issue and within the framework of rate-independent, nondegrading symmetric hysteresis, a new class of uniaxial constitutive models has been recently proposed by Vaiana et al. [20]. Such materials are capable of modeling a broad variety of hysteretic behaviors by means of smooth curves whose expression is analytically determined [19,21].
From a qualitative point of view, such a class provides hysteresis loops quite similar to the Bouc-Wen model, at least in terms of hysteresis shape and of hardening-softening phenomena, although it is essential to remark that their analytic background is completely different. Implementation of models belonging to such a class is very convenient since they do not require to solve any differential equations; hence, they are, in general, faster and less computationally demanding than Bouc-Wen-based models. Moreover, the possibility of expressing the response by means of explicit nonlinear equations makes such a class very convenient to address applications relying on closed-form solutions, such as dynamic viscoelastic problems [8], and especially to determine the response sensitivity in closed form [18].
A further benefit of the models belonging to the class presented in Vaiana et al. [20], and particularly of most versatile one, the algebraic material developed by Vaiana et al. [21], consists in the fact that their constitutive parameters can be directly related to quantities deduced from experimental data so that their identification is straightforward, as discussed in Sessa et al. [18].
In order to ensure physical significance of the computed response, it is important to investigate the thermodynamic compatibility of the model and its stability.
Thermodynamic compatibility consists in the fulfillment of the thermodynamic principles by means of the postulate formulated by Drucker [7]. Specifically, such a condition prescribes that the first derivative of the energy dissipated by the hysteretic model is always non-negative. Moreover, the stability postulate introduced by Hill [10] is a stronger condition for which it is necessary that the incremental internal energy of the system must be nondecreasing. For generic hysteresis loops, both the conditions usually depend on the adopted constitutive parameters as well on the strain/displacement range.
Materials that fulfill both these criteria are generally well-suited for numerical analysis, while models that fail to satisfy such criteria are likely to present drawbacks, such as nonuniqueness or singularity, during the solution process.
It is worth being emphasized that Hill's stability criterion is never fulfilled in presence of softening, although modern integration algorithms are capable of determining stable responses even for materials not fulfilling it. Obviously, if softening is introduced, the designer must be aware of the instability of the model and should account for it in the response computation. On the contrary, thermodynamic compatibility represents a fundamental condition in order to avoid physically impossible energy balances of the analyzed material.
The research hereby presented introduces simple analytical conditions that must be fulfilled in order to ensure the thermodynamic compatibility of rate-independent hysteretic models and investigates their specialization for the case of the algebraic material presented in Vaiana et al. [21] and belonging to the class formulated in Vaiana et al. [20].
In particular, a general purpose condition for the fulfillment of the Drucker's postulate is discussed. This can be adopted for any of the models belonging to the class Vaiana et al. [20] and even for different materials, such as the celebrated Bouc-Wen. Such a condition is based on the computation of the incremental energy of a generic hysteresis loop and on the evaluation of its reversible part. This permits to develop a conservative analytical condition based on the computation of the tangent and secant stiffness of the analyzed model. Such a general condition is therefore specialized for the case of the algebraic material proposed in [21] and the determination of analogous compatibility conditions defined by means of the constitutive parameters.
The computed conditions are finally tested by numerical applications discussed before a few conclusions.

Thermodynamic compatibility of a generic hysteretic material
Let us consider a generic state of a uniaxial hysteretic material defined by a current strain value ε, supposed positive for simplicity, and the associated stress value σ (ε). We emphasize that strain and stress can be interpreted as generalized quantities in order to refer to displacement and forces.
The first derivative dσ (ε)/dε of the stress σ with respect to the strain ε is the tangent stiffness represented by the slope of the hysteretic curve at the considered point, as shown in Fig. 1.  Moreover, it is possible to define the unloading curve associated with the state [ε, σ (ε)], represented in red in Fig. 1, and the residual strain ε 0 (ε) corresponding to the intercept of the unloading curve with the axis σ = 0 and associated with the tangent stiffness dσ (ε 0 )/dε.
The condition defined by the Drucker's postulate [7] states that the increment of inelastic work must be non-negative. In particular, making reference to the case of plasticity as example, such a relationship is defined by means of the plastic strain: where the · operator represents the scalar product, being stress and strain tensorial quantities. Moreover, dε pl is the increment of plastic strain ε pl , i.e., the difference between the total strain ε and the residual strain ε 0 (ε). From a physical point of view, Equation (1) states that the second degree plastic work is always positive. Such a condition is less limiting than the Hill's postulate [10] but is far more significant so that it is often addressed as thermodynamic compatibility. In fact, plastic work corresponds to the energy which is dissipated by the constitutive model. In order to fulfill the thermodynamic principles, it is necessary that such a quantity is never negative. Whenever a stress-strain state does not fulfill Equation (1), the system presents the absurd (and physically impossible) condition for which its free energy turns out to be greater than the energy introduced by the external action. Therefore, while the Hill's stability postulate represents a useful condition to obtain stable materials from a computational point of view, the thermodynamic compatibility condition (or Drucker's postulate) of Equation (1) represents a fundamental constraint to be fulfilled in order to reproduce physically meaningful constitutive behaviors.
To investigate the fulfillment of the thermodynamic compatibility condition it is convenient to define such a relationship by means of inelastic work. In fact, computational models belonging to the class proposed in Vaiana et al. [20] have phenomenological nature; thus, inelastic phenomena are not necessarily related to plasticity.
Such an inelastic work can be computed as the difference between the total work W tot associated with the stress-strain state and the work associated with the unloading curve starting at the current state [ε, σ (ε)], which is represented by the yellow area in Fig. 2.
This latter quantity is the energy that the material releases if the external action is suppressed. For simplicity, such a contribution will be addressed as elastic energy below.
Computation of the inelastic work increment at a stress-strain state [ε, σ (ε)] depends on the algebraic sign of the strain increment dε. In particular, the following subsections investigate the cases of negative and positive ratios separately.  It is worth being emphasized that, because of the phenomenological nature of the considered models, the uniaxial relationship may refer to stress and strain as well as to forces and displacements. For this reason, from now on, strain-like quantities will be addressed as generalized displacements.

Inelastic work increment for dε < 0
In case of negative displacement increments, the inelastic work increment is trivially zero. In fact, the total work increment associated with dε, represented in yellow in Fig. 3, turns out to be negative and it is equal to the released elastic energy since the material follows the unloading curve. Hence, for dε < 0 the thermodynamic compatibility is always fulfilled.

Inelastic work increment for dε > 0
In case of positive displacement increments dε > 0, the response follows the blue line. The total work increment, computed by means of its second-order Taylor series expansion, is: where dσ + (ε)/dε denotes the tangent operator computed for positive displacement increments on the loading curve, represented by the blue line in Fig. 4 where the quantity dW tot is represented by the blue area. Analogously, the increment of elastic energy, represented by the green area in Fig. 5 and reported in detail in Fig. 6, can be computed by means of a Taylor series expansion. Specifically, the load increment corresponds to the height of the green area and turns out to be: while the quantity δh is evaluated by means of the tangent stiffness computed along the unloading curve: where dσ − (ε)/dε represents the tangent operator computed on the unloading curve (red line of Fig. 5). Hence, the elastic energy increment turns out to be: and, finally, the inelastic work increment is computed as the difference of the quantities evaluated by Equations (2) and (5) as: We recall that the displacement increment is a positive infinitesimal quantity dε → 0 + ; thus, the following infinitesimal quantity: turns out to be of higher order with respect to the remaining terms within the curly brackets in Equation (6). Thus, the first derivative of the inelastic work turns out to be: To determine the algebraic sign of dW in /dε, it is necessary to determine (ε 0 (ε) − ε) which yields: wherek is the average stiffness computed along the unloading curve: The latter quantity can be interpreted as the slope of the regression line matching the unloading curve; i.e., a segment which subtends equal areas with respect to the unloading curve, as shown in Fig. 7.
Recalling Equation (9), the thermodynamic compatibility condition of Equation (8) becomes: It is worth being emphasized that such a condition must be fulfilled for all possible load-displacement states which are, eventually, infinite. With this respect, the quantity dσ (ε)/dε represents the tangent stiffness at an arbitrary point of the hysteretic curve, whilek can span within all possible values that the secant stiffness k s of the unloading curve can attain. Hence, the thermodynamic compatibility condition can be generalized as: The latter relationship aims to avoid thermodynamically not-compatible behaviors such the one represented in Fig. 8 where, provided an equilibrium path denoted by the blue line, we consider the work increment between the states ε and ε + ε. In this particular case, the increment of the reversible work W el , represented by the area in yellow underneath the unloading curve in red, turns out to be greater than the total work increment W tot , represented by the cyan area.
Along with the previous one, a further compatibility condition should be checked. In fact, in presence of softening material, the hysteresis loop may intercept the horizontal axis determining states for which even the total work W tot turns out to be negative, as for the example reported in Fig. 9 where the negative work increment is represented by the cyan area.
Condition of Eq. (12) is a general purpose relationship since it does not refer to peculiar formulations of the constitutive model. Hence, it can be used to investigate any hysteretic material, as long as the computation of the tangent and secant stiffness is possible, and can be applied to all models belonging to the class formulated in Vaiana et al. [20].

Thermodynamic compatibility conditions of the algebraic model
The algebraic material Vaiana et al. [21] is a constitutive model capable of reproducing hysteresis loops with complex shapes. It is defined by means of five constitutive parameters. In particular, k a and k b are stiffness coefficients relevant to the elastic and inelastic regions, respectively, while parameter α determines a smooth transition between elastic and inelastic states. Two further parameters, namely β 1 and β 2 , introduce hardening/softening behaviors, pinching and local minima/maxima of the response.
Within the framework of a dynamic or pseudo-dynamic structural analysis, the values of the generalized displacement, velocity and response at the last computed time (or pseudo-time) step are denoted by u c , s c and f c , respectively. In order to compute the response f t associated with a trial generalized displacement u t , it is necessary to introduce the following auxiliary parameters: and the history variable: with: and: Hence, the trial response is: with: and: while the tangent operator, i.e., the derivative of the response with respect to the generalized displacement, is: It is important to recall the physical interpretation of the constitutive parameters since it represents a pivotal issue of the compatibility conditions presented in Subsect. 3.1. Specifically, if β 1 = β 2 = 0, parameter k a represents the tangent stiffness at the beginning of each loading curve, while k b is the asymptotic tangent, as shown in Fig. 10a. Moreover, coefficient α rules the transition between the loading branch and the inelastic region. Such a parameter is not directly represented by the hysteresis loop although it determines the auxiliary parametersf and u 0 of Equations (13) and (14). The quantityf represents the half-distance between the two asymptotic branches, while u 0 is a half of the displacement for which the nonlinear response reaches the asymptotic branch, as shown in Fig. 10a.
Parameters β 1 and β 2 determine the nonlinear hardening behavior, as shown in Fig. 10b where a set of cycles with different values of such parameters have been represented.

Thermodynamic compatibility conditions
In order to investigate the thermodynamic compatibility of such a constitutive model by the condition of Equation (12), it is necessary to compute the maximum of the tangent stiffness of Equation (22) and the minimum of the secant stiffness.
Referring to the original paper [21] for the details, we recall that the maximum of the tangent stiffness is always greater than parameter k a . In fact, such parameter represents a local maximum reached nearby the regions of the hysteresis loops where the velocity changes its algebraic sign. Further local maxima of the tangent stiffness can be achieved within the inelastic region of the loop for which the last addend of Equation (22) tends to zero. In such a case, local maxima can be computed as the stationary points of the function: where ε denotes the generalized displacement. The secant stiffness is easily computed by dividing the response of Equation (19) by the trial displacement. Assuming that velocity is positive, it yields: The latter equation can be simplified recalling that the minimum of the secant stiffness is attained for high values of the displacements. In such a case, the last term of the right-hand side of Eq. (24) is negligible so that The stationary points of the tangent and secant stiffness depend upon the values of the constitutive parameters. Recalling that k a > 0 and k b < k a , it is possible to determine typical cases depending on β 1 and β 2 .

Condition for the case β
If the two modulating parameters β 1 and β 2 are both zero, the secant and tangent stiffness present the trend reported in Fig. 11 so that their minimum and maximum values turn out to be: hence, the compatibility condition depends on k a and k b only and is not dependent by the displacement range: Whenever the latter condition is not fulfilled, it is still possible to determine a range of the generalized displacement for which the thermodynamic compatibility is ensured. Specifically, being k a the maximum of the tangent stiffness and on account of Equations (12) and (25), the displacement range ε is: 3.1.2 Condition for the case β 1 = 0 and β 2 < 0 If the modulating parameter β 2 is negative, the tangent stiffness decreases monotonically and its maximum is k a , as reported in Figs. 12a and b. Nevertheless, the minimum of the secant stiffness is −∞; thus, it is not possible to determine conditions about the values of the parameters for which the thermodynamic condition is always fulfilled. Yet, it is possible to determine a compatibility range ε of the displacement by accounting for Equations (12) and (25) which yields: In such a case, the displacement compatibility range ε must be determined numerically because the latter equation has no closed-form solutions. An approximate root of Equation (30) can be determined by observing that, if the inelastic region of the hysteresis has been reached, thenf /ε ≤ k a . In such a case, recalling that such ratio is always positive, a stronger condition is: so that the displacement range ε is determined by the smallest real root of the following: It is worth being emphasized that the latter equation always provides at least one real root. In fact, since β 2 is negative and |β 1 | < β 1 2 − 4β 2 (k b + k a /2) , the argument: is always positive.

Condition for the case β 1 < 0 and β 2 > 0
The third analyzed case assumes β 1 < 0 and β 2 > 0. Both the secant and tangent stiffness tend to diverge with positive values, as shown in Fig. 13a. The minimum of the secant stiffness cannot be computed in closed form although it is greater than the minimum of the tangent stiffness fulfilling the condition: and attained for ε = √ −3β 1 /10β 2 .
Thus, the minimum secant stiffness can be approximated as: Besides the global maximum of +∞, the tangent stiffness presents the local maximum: k a attained at ε = 0. Hence, a first compatibility check consists in the following: . 13 Trends of the tangent and secant stiffness for β 2 > 0 associated with the local minimum of k(ε). Moreover, it is possible to determine the displacement range by fulfilling the condition: which yields: with: We emphasize that, for peculiar values of the parameters, Equation (38) provides complex roots. In such a case, the displacement compatibility range ε must be determined by numerically.  Fig. 12b. Nevertheless, the minimum of the secant stiffness cannot be determined in closed form since a condition analogous to Equation (34) would have a complex root.

Condition for
A conservative approximation of the secant stiffness consists in assuming: Being k a a local maximum of the tangent stiffness at ε = 0, a first condition analogous to Eq. (28) yields: If such a condition is fulfilled, accounting for Eq. (23) a further compatibility condition for ε > 0 yields: and the displacement range turns out to be: where, being β 2 > 0 and k b > 0 because of Eq. (41), the condition provides at least one real root.
S. Sessa A further case assumes β 1 < 0 and β 2 = 0. As shown in Fig. 14a, both the tangent and secant stiffness are monotonically decreasing. Hence, their maximum and minimum turn out to be: Again, a displacement range can be computed by recalling Equations (44) and (25), and the condition of Eq. (12) becomes:f which has three solutions in the complex space. The sole real solution determines the displacement range which turns out to be: with: 3.1.6 Condition for the case β 1 > 0 and β 2 = 0 The last analyzed case is relevant to β 1 > 0 and β 2 = 0, and the trends of the tangent and secant stiffness are reported in Fig. 14b. Hence, the minimum of the secant stiffness turns out to be: and is attained for ε min = 3 f /2β 1 .
The tangent stiffness has a local maximum of k a and the global maximum at ∞. Thus, as long as the following condition yields: the displacement range, computed by recalling Equations (23) and (12), turns out to be:

Softening thermodynamic condition
Thermodynamically not-consistent conditions may arise whenever the softening branch of the hysteresis crosses the horizontal axis so that increments of the total work are negative, as in the example of Fig. 9.
Since all parameters but k a and α can assume negative values, the algebraic material may present loops with noncompatible softening branches. Hence, it is necessary to perform a further compatibility check.
The algebraic constitutive model is analytically formulated so that the response has two asymptotic curves c u and c l representing the upper and the lower bound, respectively, of all the possible equilibrium states. A scheme of such boundaries is reported in Fig. 15.
Asymptotic curves are symmetrical and defined by fifth-degree polynomials, specifically: The softening compatibility conditions are defined as: hence, in general, their solution must be determined numerically. As a matter of fact, such a condition is far easier to be checked with respect to the thermodynamic compatibility introduced in Sect. 2 and investigated in Sect. 3.1. In fact, besides checking the shape of the loop, a simple control on the algebraic sign of generalized force and displacement is sufficient to avoid non admissible responses.

Numerical applications
To investigate the thermodynamic compatibility of the algebraic material, seven set of parameters have been used to determine the structural response and the energetic quantities of interest.
Specifically, five sets of parameters, denoted by AM1-AM5 and not associated with real materials, have been chosen in order to be significant of all the conditions exploited in Sect. 3.1 with exception of those with β 2 > 0 which are investigated by two further models based on experimental loops.
The first one has been calibrated on the experimental response of a fiber-reinforced rubber bearer (FREB) tested by Kelly and Takhirov [11], while the second one (test PW16010L) corresponds to the longitudinal experimental response of a wire-rope isolator type PWHS 16010 tested in Vaiana et al. [22]. The numerical values of the parameters relevant to such experimental materials have been calibrated by the identification procedure proposed in Sessa et al. [18]. Analyses have been performed by using the finite element framework OpenSees [14] in which the algebraic material has been implemented. A release of the framework can be downloaded at Sessa et al. [17] where example input files of some analyzed materials are also available.
All the sets of parameters considered in the numerical applications are reported in Table 1 whose last column reports the maximum displacement range ε for which the thermodynamic compatibility is fulfilled.
For each set of parameters, the load-displacement response has been computed by performing static analyses consisting in four cycles with increasing amplitude. The maximum displacement associated with  Figs. 16a, 17a, 18a, 19a,20a, 21a, 22a.
The work of the static response associated with each step of the analyses is the total work which has been introduced in the mechanical system. Moreover, the energy dissipated by the model has been computed as the difference between the total energy and the work associated with a monotonic unloading curve. These are reported in Figs. 16b, 17b, 18b, 19b,20b, 21b, 22b as functions of the displacement. Moreover, in these figures, the displacement boundary ε has been represented by dotted lines.
As a matter of fact, such inelastic work turns out to be nondecreasing for all the considered models and, despite the fact that the displacement overcomes its boundary for models AM1 FREB and PW16010L, the load-displacement response does not present any instability or violation of the thermodynamic compatibility since the inelastic work remains nondecreasing.
Such an aspect confirms the conservative nature of the bounds computed in Sect. 3 consisting in the fact that conditions reported in Equations (29), (32) and (38) are sufficient but not necessary to ensure thermodynamic compatibility.   Displacement paths adopted in the applications of the previous section were regular and are relevant to quasi-static tests. In order to exploit the compatibility conditions in a common practice context, time history analyses have been performed and are hereby presented.

Dynamic analysis
In particular, two single degree-of-freedom oscillators have been considered. These adopt the FREB and PW16010L constitutive relationships with mass of 3000 Kg and 1000 Kg, respectively. No additional damping has been introduced. Such models have been selected for the time history analysis since they are relevant to real seismic devices. The remaining models reported in Table 1, while useful to test the compatibility conditions, may be not effective in dynamic analyses adopting real accelerograms.
The considered base excitation is the longitudinal main shock of the Kobe 1995 earthquake, reported in Fig. 23a lasting about 40 seconds. Analyses have been performed in OpenSees by a central difference explicit algorithm.
The displacement of both the oscillator remains smaller than its compatibility range ε reported in Table  1 during all the time history analysis as shown in Fig. 23b reporting the displacement as function of time.
Moreover, force-displacement responses, reported in Figs. 24a and 24b, show that the oscillators are subject to load paths resulting far more irregular than the quasi-static ones of the previous examples.
As a matter of fact, the outcome of the time history analyses confirms the effectiveness of the compatibility conditions provided in Sect. 3.1. Indeed, while the total work, reported as function of time in Fig. 25a, presents decreasing regions, the inelastic work, reported in Fig. 25b, is a monotonically increasing curve.

Conclusions
The thermodynamic compatibility of the algebraic material presented by Vaiana et al. [21], intended as the fulfillment of such postulate formulated by Drucker [7], has been investigated. In particular, it is necessary that the energy dissipated by the model, denoted for convenience as inelastic work, must be nondecreasing.
In general, such a condition is not fulfilled by any possible set of constitutive parameters since it is proved that it depends on both the maximum tangent stiffness and the minimum secant stiffness that the model can attain within an interval.
If the values of the constitutive parameters do not ensure the fulfillment of the Drucker postulate for any possible value of the displacement, a compatibility range can be defined as function of the constitutive parameters. This ensures that the Drucker postulate is not violated as long as the displacement does not exceed such a boundary.
Numerical applications have been presented in order to test the compatibility conditions defined in closed form. In particular, seven constitutive models have been analyzed by means of quasi-static numerical analyses in which the total and inelastic work have been determined.
Results show that all the models present nondecreasing inelastic work regardless of the attainment of the displacement boundary. Thus, the compatibility displacement threshold represents a strongly conservative condition so that the behavior of the constitutive model does not necessarily become unstable if the boundary is exceeded.
Two further time history analyses, adopting the Kobe 1995 earthquake acceleration as base excitation, have been performed in order to test the compatibility conditions in case of irregular and strongly oscillating loads consistent with real environments. Again, the inelastic work turns out to be monotonically increasing, thus confirming the outcomes of the quasi-static tests.
The computed conditions represent a useful tool to assess the compatibility of the algebraic material [21] when it is adopted to perform dynamic analyses in which the displacement path is not defined a priori.