This section discusses the optimisation-based procedures employed for the identification of the constitutive parameters of the finite strain viscoelastic material model coupled with gradient-enhanced continuum damage as introduced in
Section 3. Since the general framework can be applied to all kinds of energy functions, at first, considering the type of material and its behaviour, the chosen hyperelastic material model is briefly presented in
Section 4.1. Next, the applied damage function
is specified in
Section 4.2. The main idea is to first carry out a parameter identification based on homogeneous states of deformation as described in
Section 4.3. The obtained set of material parameters is then used as the initial value for the parameter identification based on inhomogeneous states of deformation, as elaborated in
Section 4.4. The gradient terms of the gradient-enhanced damage framework are activated only in the case of inhomogeneous deformations.
4.1. Hyperelastic Material Model
Considering the nearly incompressible, non-linear elastic material behaviour (see for example the stress-strain-response in
Figure 2b), the Yeoh hyperelastic material model was chosen [
25,
26]. In [
27], Yeoh and Flemming argued that Rivlin’s theory [
28], where the strain energy function depends on the first two invariants of the right Cauchy–Green tensor
, to be specific
runs into difficulties with respect to the identification of the material parameters since the contributions cannot be perfectly determined separately in the experiments. Here,
are the corresponding material parameters, and
and
are the first and second principal invariants related to the isochoric right Cauchy–Green tensor, respectively. Yeoh neglected the contribution of the second invariant assuming that no serious error would arise since the contribution would be sufficiently small with respect to the impact of the first invariant. In this work, however, the Yeoh energy potential
serves as the isochoric contribution to overall energy potential, where
denote the underlying material parameters. This phenomenological material model was developed for rubber elasticity and is applied in this work as a first choice with respect to the mentioned material behaviour. The volumetric energy contribution
is chosen in a standard manner as
where
K denotes the bulk modulus of the material. The related energy contributions
and
can, however, be straightforwardly replaced within the general constitutive framework elaborated in this work in order to account for different particular material characteristics.
4.2. Damage Function
The damage function mentioned in
Section 3 was chosen to be of an exponential type, following [
6,
29]. Thus, a damage initiation threshold, as well as a damage saturation rate can be used in the function
fulfilling the requirements for the damage function to be restricted to
. Here,
d denotes the classic damage variable with
indicating zero damage and
for
damage. The variable
represents a related internal damage variable, while
is the damage threshold parameter, and
introduces the damage saturation parameter. In addition,
represents the Macaulay brackets. Consequently, damage is obtained in the material if the local damage variable exceeds the damage threshold. The damage function for the isochoric energy contribution matches this chosen function except for the exponent
. Thus, a relation between both functions is still given and only adjusted via this material parameter. The evolution of the internal damage variable
is based on the associated form
where
denotes a proper Lagrange multiplier,
represents the damage condition, and
is the energy release rate. Hence, apart from the balance of linear momentum, the Euler–Lagrange equation governing the non-local damage variable needs to be solved simultaneously. For further information regarding this regularised damage model and its implementation, see Ostwald et al. [
6]. The implementation was performed in a user material subroutine (UMAT) in Abaqus following [
6].
4.3. Parameter Identification Based on Homogeneous States of Deformation
In the context of parameter identification procedures, it is advantageous to compute the Piola stress tensor as a quantity, facilitating the comparison of experimental and simulated material response; cf.
Section 2. The Piola stress tensor
is related to the Piola–Kirchhoff stress tensor
and the Kirchhoff stress tensor
via
where the sensitivity of the Piola stress with respect to the deformation gradient is denoted as
, in particular,
or, expressed in index notation,
Moreover,
relates to the spatial elasticity tensor
via
It is noted that the tangent operator, as typically required within implicit finite element formulations, respectively so-called constitutive drivers, is based on .
The Yeoh-type hyperelastic energy function was combined with the constitutive viscoelasticity model, presented in
Section 3, including two Maxwell elements. Apart from the Yeoh material parameters
,
and
, the viscoelastic parameters,
and
, denoting the relation of the stiffness of each Maxwell element with respect to the pure elastic Young’s modulus, and
,
, representing the ratio of the viscosity to the stiffness in each Maxwell element, need to be identified. Regarding the parameters, it has to be mentioned that the material was assumed to be nearly incompressible, and thus, Poisson’s ratio was fixed to
. Thus, seven material parameters needed to be identified.
Since these seven material parameters could be identified via experiments of homogeneous deformation states, strain- and stress-driven constitutive drivers were implemented in MATLAB instead of full finite element (FE) simulations in order to reduce the computational cost within each iteration of the parameter identification. The pseudo-codes of both types of constitutive drivers are depicted in Algorithms 1 and 2. The fminsearch-algorithm in MATLAB was used for the parameter identification, minimising the goal function with respect to the difference in the simulated and experimental reaction force over all load steps, i.e.,
where
,
, and
denote the weighting factors for the tensile, creep, and relaxation test, respectively. The additional weighting factors
,
, and
were introduced to emphasise specific time steps
t of each experiment.
In the following, the viscoelastic material parameters are identified via the three tests based on homogeneous deformation states, presented in
Section 2.1.
If the parameters of the Yeoh-model were solely fitted with respect to the tensile test, the simulated material response of the identified parameters matched the experimental data perfectly. Next, the relaxation parameters were obtained with respect to the relaxation and creep test, though the material behaviour highly depended on the previously identified Yeoh parameters. Thus, the experimental relaxation and creep response could not be sufficiently matched via optimisation of the relaxation parameters only. Consequently, all of the material parameters were identified simultaneously for all three different experiments in order to obtain the best parameter set reflecting the complete material behaviour. The comparison of the simulated and experimental material behaviour is shown in
Figure 6 for the optimised parameter set. The corresponding parameter set is presented in
Table 1.
In the graphs shown in
Figure 6, the experimental results of the relaxation and creep tests are nearly perfectly matched, and only slight differences in the diagrams can be seen at certain points. The simulated material response of the tensile test overestimated the experimental curve, but captured the general trend. The tensile test could be captured better by increasing the associated weight within the objective function (
33) at the cost of the accuracy with which creep and relaxation tests were captured. The focus here was set on the rate-dependent constitutive characteristics that were captured with high accuracy. Though beyond the scope of this work, a precise capturing of tensile, creep, and relaxation test at the same time could be achieved by introducing higher order energy functions and a larger number of independent Maxwell branches within the constitutive framework provided. The accompanying increase in parameter identification complexity due to a significantly increased number of constitutive parameters could then be dealt with using parameter correlation matrices.
Algorithm 1 Strain-driven constitutive driver (relaxation)—uni-axial stress state; see also Appendix A. |
- 1:
initialize internal variables . - 2:
get material parameters . - 3:
initialize . - 4:
for every time step do - 5:
given: the stretch . - 6:
while do - 7:
compute , , , , . - 8:
compute initial elastic stress Kirchhoff tensor . - 9:
compute algorithmic internal variables , . - 10:
compute according to ( A6). - 11:
compute the spatial tangent moduli . - 12:
compute the Piola–Kirchhoff stress tensor . - 13:
compute the tangent operator . - 14:
partition the stress tensor and the deformation gradient: - 15:
- 16:
. - 17:
partition the tangent operator: - 18:
. - 19:
update the transverse deformation gradient: - 20:
. - 21:
end while - 22:
assemble the deformation gradient and the stress tensor - 23:
update internal variables . - 24:
end for
|
Considering the curve of the tensile test, especially the values and the signs of the Yeoh parameters are important, since the first parameter weights the linear behaviour of the stress strain-response, whereas scales the quadratic term of the isochoric energy contribution resulting in the decreasing slope. The third parameter corresponds to the third-order polynomial contribution of the energy function and thus yields the increasing slope at the end of the loading path. Thus, different values for the Yeoh parameters would probably fit the tensile response way better; the dependence on the relaxation parameters, however, would result in a worse fit of the relaxation and creep behaviour. Higher order contributions to the Helmholtz free energy, in particular the isochoric contribution, would therefore provide a potential in order to capture the tensile test response better.
Algorithm 2 Stress-driven constitutive driver (creep)—uni-axial stress state; see also Appendix A. |
- 1:
initialize internal variables . - 2:
get material parameters . - 3:
initialize . - 4:
for every time step do - 5:
given: Piola–Kirchhoff stress tensor . - 6:
while do - 7:
compute , , , , . - 8:
compute initial elastic Kirchhoff stress tensor . - 9:
compute algorithmic internal variables . - 10:
compute according to ( A6). - 11:
compute the spatial tangent moduli . - 12:
compute Piola–Kirchhoff stress tensor . - 13:
compute tangent operator . - 14:
compute . - 15:
update the deformation gradient - 16:
- 17:
end while - 18:
update internal variables - 19:
end for
|
The identified viscoelastic material parameters were used as fixed values for the subsequent optimisation of the damage-related parameters.
4.4. Parameter Identification Based on Inhomogeneous States of Deformation
In the following, the previously identified Yeoh and relaxation parameters are used to identify the damage-related material parameters
,
,
, and
via the tensile test with inhomogeneous deformation states presented in
Section 2.2. In contrast to the identification of the viscoelastic parameters, the constitutive driver is not sufficient enough to capture the material behaviour of the gradient-enhanced damage model. Thus, a finite element (FE) formulation is required. Therefore, the material model was implemented into a UMAT in Abaqus, as already mentioned in
Section 4.2.
In addition, a parameter identification tool was implemented in Python, using the Nelder–Mead-Simplex algorithm of the scipy-package. In contrast to the calibration of the Yeoh and relaxation parameters, it did not suffice to solely consider the difference in the Piola stress
or the (clamping) displacement
, depending on the experiment, for the goal function
of the optimisation. Since the gradient terms—capturing the mesh-objectivity of damage—of the regularised damage framework are activated by inhomogeneous deformation states, the difference in the displacement field between the simulated and experimental material response needs to be added to the goal function apart from the difference in the reaction force
where
denotes the number of node points considered,
the weighting factor for the displacement contribution, and
the weighting factor for the impact of the reaction force. The number of node points was the element nodes on the surface of the specimen in the FE simulations. The experimental data was interpolated onto those node points via a 2D-interpolation scheme prior to the parameter identification following Kleuter [
13].
For the purpose of reducing the computational cost of the FE simulation within each iteration of the parameter identification, the symmetry properties of the sample were used; cf.
Figure 7a. Furthermore, considering the material properties of the soft polymer, only the segment shown in the figure was used to exclude boundary effects of the clamping jaws. In order to still use the experimental boundary conditions in the simulations, the experimentally measured displacements, taken via the CCD camera system, were applied to the right boundary of the specified segment. As can be seen from the photos of the specimen during the experiment (cf.
Figure 5) and from the sketch of the segment, the displacements were not uniform over the boundary of the chosen segment of the sample.
The initial guess for the damage-related model parameters was
,
,
, and
and resulted in no damage initiation in the material; cf.
Figure 8a. In
Figure 7b, the comparison of the load-displacement curves is presented. Since the stress-strain path of the homogeneous tensile test already overestimated the experimental curve, the simulated response for the tensile test with inhomogeneous deformation states lied above the experimental result as well. Furthermore, considering the large total stretch values present in combination with the identified Yeoh material parameters, the increasing response for the reaction force was triggered by the parameters
and
weighting the hyperelastic energy contributions of second and third order.
Since parameter
strongly influences the simulated deformation behaviour of the specimen in each iteration,
was fixed to
in a first step, thereby neglecting a different damage contribution of the volumetric and isochoric part. The corresponding optimised parameter set, denoted as Final #1,
,
, and
, was rounded to five decimal places. Considering the load-displacement curves, the response of the final set lied closer to the experimental behaviour than the initial guess; see
Figure 7b. The stress distribution of the initial guess and the final set was comparable; the magnitude of the final set, however, was less than half of the stress of the initial set; cf.
Figure 9. The von Mises stress distribution within the specimen—with the maximum stress value obtained at the upper surface, corresponding to the region dominated by the initially circular notch—was in line with observations made by, e.g., Kleuter [
13].
The contour plot of the damage value
d in
Figure 8 compares the damage value obtained for the initial parameter set with the damage value obtained for the Final #1 set. Initially, no damage was initiated, and with the Final #1 set, a damage value of nearly
occurred in the symmetry plane of the loading direction.
In
Figure 7b and
Figure 10, the influence of the different damage functions for the volumetric and isochoric contributions is visible. The set Final #2 included
, and Final #3 considered
. A higher value for
flattens the load-displacement curve, while a lower value provides an increased path. Another important impact of the parameter is indicated in
Figure 10. Apart from the difference in the stress magnitude, the deformation of the sample was already different for slight changes in the parameter with regard to the necking of the sample at the symmetry plane of the loading direction.
The advantage of the regularised damage framework is the mesh objectivity of the simulated material response. To demonstrate this feature, the boundary value problem was calculated with three different meshes for the optimised parameter set. In addition to the mesh containing 1548 elements, used for the results in
Figure 7,
Figure 8,
Figure 9 and
Figure 10, a coarser mesh with 1092 and a finer mesh with 2613 elements were used in order to analyse the mesh sensitivity of the results. As can be seen in
Figure 11, only marginal differences in the contour plots of the damage function values and the von Mises stress, as well as the load-displacement curves are visible. Thus, the gradient-enhanced damage model is working properly. In contrast, it is noted that the local damage model diverged at different load steps, depending on the mesh discretisation. In the case of the fine mesh, the solution diverged at a displacement of
mm, for the basis mesh at
mm and for the coarse mesh at
mm. The displacement was taken from the bottom node of the right boundary in the symmetry plane.