Dissipation energy based parameter identification of anisotropic linear viscoelastic composites

The current work presents a relaxation analysis based procedure to identify effective material parameters of the multiaxial generalized Maxwell model (GMM) by a numerical homogenization of the microscopic dissipation energy density for anisotropic linear viscoelastic composites. The employed GMM enables the derivation of a thermodynamically consistent constitutive law and a function of the dissipation energy density for direction‐dependent viscoelastic materials. The identity of this energy function to the microstructure's homogenized dissipation energy density is utilized for the identification of effective relaxation times. Furthermore, the identified relaxation times enable a simple determination of the remaining stiffness parameters. Finally, the presented procedure is demonstrated and evaluated for a randomly endless fibre‐reinforced plastic with a polymer matrix exhibiting a significant viscoelastic behaviour.

The current work presents a relaxation analysis based procedure to identify effective material parameters of the multiaxial generalized Maxwell model (GMM) by a numerical homogenization of the microscopic dissipation energy density for anisotropic linear viscoelastic composites. The employed GMM enables the derivation of a thermodynamically consistent constitutive law and a function of the dissipation energy density for direction-dependent viscoelastic materials.
The identity of this energy function to the microstructure's homogenized dissipation energy density is utilized for the identification of effective relaxation times. Furthermore, the identified relaxation times enable a simple determination of the remaining stiffness parameters. Finally, the presented procedure is demonstrated and evaluated for a randomly endless fibre-reinforced plastic with a polymer matrix exhibiting a significant viscoelastic behaviour.

K E Y W O R D S
anisotropic viscoelasticity, dissipation energy function, homogenization, parameter identification, relaxation analysis

INTRODUCTION
For the description of the homogeneous anisotropic linear viscoelastic material behaviour of a heterogeneous microstructure, a tensorial stress-strain relation by means of a convolution integral including a multiaxial relaxation function has been proved viable [1,2]. Herein the relaxation modulus consists of anisotropic stiffness tensors and scalar relaxation times. For the identification of these material parameters, a multiaxial relaxation function is required, which can be determined by the evaluation of homogenized time-dependent stress and strain components of a considered microstructure. Moreover, multiple optimizations including the time-dependent coefficients of the relaxation function have to be performed in order to identify the stiffness tensors and relaxation times. The number of optimizations depends on the class of material symmetry, whereas the scalar relaxation times appear redundantly and a separate identical identification of them is not a requirement [3]. Another possibility is their prior estimation, what results in a good adjustment of the stiffness coefficients for special cases (see [4,5]). However, for fibre-reinforced plastics (FRP) exhibiting a high degree of anisotropy and considerable viscoelastic behaviour on short time scales, own studies revealed, that existing parameter identification methods are inapplicable. Within an optimization algorithm including the scalar relaxation times according to [3], they did not match for different directions and their estimation or adoption from the matrix material according to [4] or [5] led to inaccurate adjustments of the stiffness coefficients. Since the reinforcement of a viscoelastic matrix with fibres changes the relaxation times [6], it is obvious that the resulting viscoelastic behaviour depends on the fibres' direction. In order to find relaxation times, which are suitable for each direction, their identification within the componentwise consideration of the relaxation function or even estimation is greatly insufficient in general. For that reason a method for the accurate determination of the scalar relaxation times by the adjustment of a dissipation energy density function for anisotropic viscoelastic materials to a corresponding scalar homogenized response, covering the direction-dependencies of an arbitrary composite under a complex deformation, is presented. Herein the equivalence between the utilized energy expression, derived from a multiaxial generalized Maxwell model (GMM) and the volume-averaged dissipation energy density according to [7] is fundamental. With the determined relaxation times, the optimization of the remaining stiffness coefficients with respect to homogenized stress and strain curves leads to a precise fitting of the entire relaxation function.
This work is structured as follows. First, Section 2 gives a theoretical background on linear multiaxial viscoelastic behaviour, mechanical mean field homogenization and parameter identification of anisotropic viscoelastic material constants. Thereby, a discussion on existing parameter identification methods for homogenized viscoelastic behaviour is presented as well. Next, the proposed algorithm for parameter identification is introduced in Section 3.1, whereas it is applicable for any composite exhibiting a homogeneous iso-or anisotropic linear viscoelastic behaviour and has been implemented by the use of the finite element software ABAQUS. An example including an evaluation of the procedure is given in Section 3.2. The current work is based on infinitesimal strains, linear viscoelastic material behaviour, constant temperatures and microstructures with ideal adherence of the phases. Finally, the reader is referenced to appendix "Nomenclature/Notation" for a better understanding of the applied notation.

Linear anisotropic and isotropic viscoelastic relations
The relation between the Cauchy stress tensor ( ) and the time derivative of the strain tensoṙ( ) for a linear viscoelastic material can be written by means of the convolution integral [8,9] Within a relaxation analysis (constant strain applied at = 0) Equation (1) becomes with the relaxation function ( ) as a positive-definite fourth-order tensor with time-dependent coefficients. The total number of independent coefficients is given by the class of material symmetry and restricted to 21, due to equivalent symmetry properties as an elastic stiffness tensor. For the case of general anisotropy, [1] proves the following choice for ( ) as viable what represents Maxwell-type relaxation behaviour as a special case of the general relaxation function ( ) given in Equation (1). Within Equation ( with ∞ and ∞ being the bulk and shear modulus, respectively, for the long-term response. Moreover, and are volumetric and deviatoric relaxation strengths, respectively, conjugated with the corresponding relaxation times and . A more detailed discussion of Equation (4) and its application to polymers can be found in [10]. Moreover, [7] reveals that the GMM enables the thermodynamic consistent formulation of the constitutive law according to Equation (1) including the relaxation function given by Equation (3). With this phenomenological model, [7] derives the corresponding dissipation power density as follows For a constant strain applied at = 0 (relaxation analysis) it follows The time integration of Equation (6) then yields the dissipation energy density for Since the fourth-order tensor is positive-definite (see [11]), the scalar-valued mapping = ⋅⋅ ⋅⋅ is greater than or equal zero for any ≠ 0. Thus, ( ) ≥ 0 holds for any ≠ 0.

Mechanical mean-field homogenization
In the following, the basic equations for the mechanical mean-field homogenization according to [12], [13] and [14] are given. For that reason, a representative volume element (RVE) of a microscopically heterogeneous structure, described as a region  with the volume and the boundary , is considered.  consists of several phases, each being homogeneous on the microscale and exhibiting a volume fraction ( ) = ( ) , where the superscript ( ) denotes the th phase. Within a mean-field homogenization of a viscoelastic composite, the following initial boundary value problem for  needs to be solved with given constitutive equations on microscale for each phase Here, Equation (8a) is the balance of linear momentum (no body forces and zero acceleration). Furthermore, Equation (8b) represents periodic displacement boundary conditions (PDBcs), with + and − being microscopic local vectors on opposite surface points. Moreover, ( , + ) and ( , − ) are corresponding microscopic displacement vectors and ( ) expresses the homogeneous strain tensor. The solution of the initial boundary value problem yields the time-dependent microscopic stress and strain tensor fields, which may be homogenized to effective stress and strain tensors by the volume averaging These homogenized results are assigned to a homogeneous substitute material, what allows the indirect identification of corresponding material parameters. To this end, the Hill-Mandel lemma integrates a physical plausibility into this homogenization scheme, what may be written as follows: Independent of the material properties of ideally adhering phases ( ) and the size of a RVE containing them, Equation (10) is fulfilled under certain boundary conditions (Bcs), whereas PDBcs (Equation (8b)) are a suitable choice for an accurate prediction of the homogeneous material behaviour. Detailed discussions on Bcs and convergence issues of effective properties within a homogenization may be found in [15][16][17][18][19]. Moreover, [7] reveals that Equation (5) is identical to the effective dissipation power density of a viscoelastic composite, which is homogenized as follows Here ( ) denotes the local microscopic volume-specific irreversible stress power of a linear viscoelastic phase ( ) with a constitutive law on the microscale according to Equation (1) including a relaxation function given by Equation (3) or (4). For the homogenized dissipation energy density it follows what needs to be equivalent to Equation (7) for a relaxation analysis of the composite under the same effective strain state .

Parameter identification of anisotropic linear viscoelastic materials
Within a parameter identification of anisotropic linear viscoelastic materials, the coefficients of the stiffness tensors ∞ and as well as the scalar relaxation times according to Equation (3) are determined. For that reason, the timedependent effective stress and strain (rate) tensors are calculated by a numerical homogenization of a microstructure to further calculate the relaxation function at several time points , which is then denoted by ( ). Since the effective relaxation function occurs in a convolution integral, special techniques for the parameter identification, including the viscoelastic correspondence principle [20], are necessary in general. The anisotropic relaxation function is then considered in the Laplace domain and the material parameters are identified for each direction according to the collocation method [1] or the multidata method [21]. Within a numerical homogenization arbitrary Bcs and therefore, according to Equation (8b), any time-dependent strain state can be applied. Thus, there is no reason to not apply a constant strain within an approximately infinitesimal time interval and conduct a relaxation analysis. This enables the calculation of the effective relaxation function's coefficients utilizing Equation (2) without additional transformations into the Laplace domain. In this manner up to 21 independent coefficients˜( ) of the relaxation function for general anisotropic behaviour can be calculated at = 1 … time points , by the application of independent strain states. The minimization of the scalar objective function then leads to the parameters and , whereas ∞ can be priorly found with˜( → ∞). The total number of relaxation times and stiffness coefficients depends on the observed time interval. The given procedure leads to accurate parameters for a wide range of anisotropic viscoelastic materials, whereas corresponding examples are given by [3][4][5]22].
On the other hand, the minimization of Equation (13) leads to a set of stiffness parameters for each coefficient of the relaxation function ( ), but the scalar relaxation times are redundantly involved within each optimization. Especially for composites with a high degree of anisotropy and a pronounced relaxation behavior, own studies revealed that different relaxation times occurred for each coefficient, what is in contradiction to Equation (3). For the considered composite, no convenient estimation of the relaxation times could be found to generate acceptable minimization results of Equation (13) for all directions. The same issue can be found, if a parameter identification for the general case is utilized. Here, the viscoelastic correspondence principle yields a similar minimization problem for each coefficient of the relaxation function in the Laplace domain (see [3]). For that reason, chapter 3 proposes an energy-based determination for a suitable distribution of relaxation times, what then results in accurate solutions for the optimization of Equation (13).

Procedure
For the proposed procedure of an energy-based identification of linear viscoelastic material parameters, an objective function based on Equation (7) is formulated as follows: The minimization of this objective function yields the parameter > 0 and > 0 within a non-trivial solution within a first step. Note that the coefficients of are not identified by this minimization. Their identification is done in a subsequent step described further below. If the applied strain tensor for the relaxation analysis contains only non-zero components, the entire material behaviour is included within the scalar energy response, what results in the determination of accurate relaxation times for the considered composite. For logarithmically spaced relaxation times, the following distribution with = − 1 = 0 … − 1 is introduced and represent the minimal and maximal relaxation times, respectively, to be regarded. For the case of one single Maxwell element ( = 1), 1 is set to .The parameter in equation (15) controls the distribution of relaxation times between and and depends on these interval boundaries. For that reason, the number of independent parameters within the minimization of the objective function (14) decreases from 2 to + 2 ( , and coefficients ). According to [3] the relative error allows the evaluation of the parameter identification. The nonlinear minimization of Equation (14) is then conducted multiple times. In each optimization, the number of Maxwell-branches is increased. In each analysis, the boundaries for and are set to the first time step and the maximum analysis time. Furthermore, Equation (16) is compared to a corresponding given maximum for each optimization with fixed . Once decreases below , the solution is taken. If this is not achieved between = 1 and a given the solution with the lowest is taken. Within a second step, the remaining stiffness parameters of the relaxation function given by Equation (3) are determined by the minimization of Equation (13), which is linear due to the already identified relaxation times. For this purpose, the considered composite is homogenized according to Section 2.2 within relaxation analyses under different strain states . The computation of the homogeneous stress tensor according to Equation (9a) for each analysis then allows the calculation of 21 time-dependent coefficients of the relaxation function utilizing Equation (2). Furthermore, the coefficients of ∞ are directly identified with ( ). As a restriction, ( ) and therefore are assumed to be orthotropic at most, what TA B L E 1 Elastic properties of T300 a carbon fibres, X -fibre direction, Y,Z -transverse directions (see Figure 1) is sufficient for a wide range of composites. Thus, only nine independent coefficients per tensor remain (see [9]). Moreover, the relaxation strengths have to be determined as positive definite, what requires the following coupled constraints for the optimization of the objective function (13) with respect to orthotropic tensors for each branch = 1 … of the GMM , , , , .

(17d)
The presented procedure is applicable for composites exhibiting an effective linear viscoelastic behaviour up to the orthotropic class of material symmetry. A lower class, e.g. isotropy or transverse isotropy, is automatically included and results in a corresponding appearance of the stiffness tensors within the relaxation function.

Example
In the following, the energy-based parameter identification according to Section 3.1 is demonstrated and evaluated for an unidirectional endless fibre-reinforced plastic by the use of the Finite Element Method (FEM). The employed software is ABAQUS/Standard [23], whereas the preprocessing (generation of the microstructure's geometry, meshing, analysis generation) and the postprocessing (homogenization, parameter identification) are realized with Python and the included packages NumPy [24] and SciPy [25]. Section 3.2.1 firstly illustrates the considered composite and the chosen RVE for it. Furthermore, the mechanical properties of the constituents are presented. Within Section 3.2.2, the results of the energybased parameter identification are shown, whereas Section 3.2.3 demonstrates that they are valid for the homogeneous description of the compound for arbitrary deformation processes.

Materials and representative volume element
The considered composite consists of carbon fibres (TORAYCA T300) and an epoxy-based matrix material. For the fibres an elastic transverse isotropic behaviour is assumed, whereas the material properties are taken from literature and illustrated in Table 1. The matrix material is a non-commercial polymer blend and was developed by the Fraunhofer Institute for Applied Polymer Research, Teltow, Germany. It combines stiffness characteristics of standard epoxy resins with an enhanced impact resistance as well as damping behavior at room temperature, what results in an increased relaxation behavior within a short time period. For its viscoelastic characterization a total number of three relaxation experiments at thirty degrees Celsius with a tensile testing machine from ZwickRoell including a climatic chamber were carried out. The constant longitudinal strain as well as the lateral strain were recorded by extensometers. The longitudinal stress was recorded with a load cell. Due to the special relaxation behavior of the polymer the experiments were conducted over a time period of 5 minutes. For the description of the matrix materials' mechanical behaviour the constitutive law given by

F I G U R E 1 SEM image, geometric model, FEM-RVE of the examined unidirectional FRP
Equation (1) including the isotropic relaxation function was assumed. Furthermore, the results of the relaxation experiments revealed that the viscoelastic behaviour of the volumetric part of Equation (4) was negligible. Thus, a constant bulk modulus was employed. The identification of the deviatoric relaxation times and stiffness coefficients was realized analogous to the nonlinear minimization of the objective function 13 for uni-axial material behaviour, whereas the remaining elastic properties (constant bulk modulus and deviatoric long-term stiffness) were directly determined with the corresponding stress-strain relations. Table 2 contains the identified material parameters, whereas no significant improvement of the parameter identification was observed with more than five relaxation times. Since all of the conducted relaxation analyses lead to similar results no averaging was taken into account and the parameters of the first experiment were chosen.
For the realistic modeling of the considered composite, a self-developed tool was utilized, which enables the conversion of a corresponding scanning electron microscope (SEM) image into a geometry file, what is further processed to a finite element RVE with desired features (size, with or without fibre contact, mesh type and mesh size). The considered endless carbon FRP was produced by resin-transfer molding (RTM), whereas a taken SEM image, a suitable large section as well as the geometry and the corresponding RVE of it are illustrated in Figure 1. To find a minimal necessary size of the RVE, a section of the entire SEM image (left in Figure 1) was taken from its origin and increased until the fibre volume fraction did not change anymore, what coincided with a convergence of the resulting homogeneous mechanical properties (see green dashed lines in Figures 3 and 4). The resulting fibre volume fraction amounted to 61.3%. Since the generated RVE represents a non-periodic microstructure, it exhibits a non-uniform finite element mesh. Thus, the application of PDBcs according to Equation (8b) was realized with an interpolation between the outer nodes similar to [27]. Furthermore, this was combined with the concept of pilot nodes according to [28] to improve the efficiency of the pre-and postprocesing within the numerical homogenization.

Results of the parameter identification
The following enumeration summarizes the parameter identification procedure according to Section 3.1 with respect to the current implementation and the example Variable Least-Squares algorithm [25] by taking the constraints (17) for identification of the relaxation strengths coefficients into account Due to the small values for the relaxation times of the matrix material (see Table 2), the relaxation analysis required relatively small time steps and only a total time of 100 seconds. The range for the number of Maxwell branches to be regarded had been chosen to be = 1, … , 12, whereas Figure 2 illustrates the obtained error within each parameter fitting for the dissipation energy function. For illustration reasons, the allowable error was set to zero to obtain optimization resulted until = = 12 is reached. As it can be seen from Figure 2, the best result is obtained for = 8 effective relaxation times, what exceeds the number of Maxwell branches for the pure matrix material ( = 5, see Table 2). The reason for this can be found in the reinforcement of the matrix by stiff endless fibres, which switch a relaxation process in the polymer to another level and time scale with respect to the fibre direction. Since the considered fibres are anisotropic, the entire direction-dependency of the composite's relaxation behaviour is more pronounced. For that reason, the stress relaxation of the compound would require different relaxation times for each direction, what is inconsistent to the given material law. However, within the evaluation of the dissipation energy density, the entire direction-dependency comes together in one scalar measurement and all necessary relaxation times are found. Their intensity corresponding to the direction is then scaled by the coefficients of the relaxation strengths .
In the current example, a number of = 8 relaxation times yielded the smallest optimization error. For higher numbers, the optimization resulted in a higher relative error, which might be attributed to the increasing number of sufficient solutions. However, more than eight relaxation times were not necessary for an excellent fitting result, as can be seen from Figure 3. It shows the numerically calculated and the adjusted dissipation energy density, F I G U R E 3 Numerically calculated and fitted density functions for energy dissipation of the examined unidirectional FRP for = 8 effective relaxation times, error: 0.05% (see Equation (16) which are nearly identical, what emphasizes the accuracy of the procedure. Table 3 contains the eight effective relaxation times of the compound. A complete presentation of the identified relaxation strengths and the longterm stiffness is omitted. Instead of this, Figure 4 shows four different comparisons between fitted and numerically calculated coefficients of the effective relaxation function. The error is calculated according to Equation (16) for each coefficient of ( ). Thereby, the maximum value for the error is 0.007% for the coefficient ( ).

F I G U R E 4
Numerically calculated and fitted functions for some coefficients of the relaxation function of the examined unidirectional FRP for = 8 effective relaxation times, error: ≤ 0.007% (analogous to Equation (16)) TA B L E 4 Identified relaxation times of the examined unidirectional FRP for = 5 Finally, the identification results for the considered FRP for the case of = 5 Maxwell branches are briefly presented as well. Table 4 summarizes the identified relaxation times. The relative error of the fitted dissipation function is 0.31% (see Figure 2). For the coefficients of the relaxation function, the worst result can be found for ( ) with an relative error of 0.02%. These errors still seem to be small. However, significant deviations can be found in simulated quasi-static analyses (see Section 3.2.3).

Validation of the parameter identification
Within this Section, the resulting sets of parameters according to Section 3.2.2 are validated. For that reason, the homogeneous material law given by Equations (1) and (3) was implemented in ABAQUS/Standard by the Fortran subroutine (UMAT) as proposed by [29]. In addition, the volume-specific dissipation power as well as energy was added within the UMAT in accordance with the procedure for the discretization of convolution integrals proposed by [29]. For the time interval from to +1 with the increment size Δ = +1 − and the initial conditions ( 0 ) = 0 as well as ( 0 ) = 0 it follows Herein, the index refers to a single Maxwell branch within the utilized GMM (see [7]). The numerical solution for the the evolution of the stress tensor ( ) can be found in [29]. For the validation, the implemented constitutive law, including the identified parameters, was applied to a single finite element (unit cell). Both the single finite element and the RVE of the microstructure were then loaded by a time-dependent homogeneous strain state prescribed by PDBcs. Finally, the comparison of the corresponding homogeneous responses was carried out. Figure 5 shows the time-dependent effective strain coefficients, which were separately applied within six different analyses. In an additional analysis, all strain coefficients were applied simultaneously. For the comparison of the results, the homogeneous stress coefficients, the effective dissipation power and the energy density were chosen. This is presented in Figures 6 and 7 for the analysis, where all six strain coefficients were prescribed simultaneously. In addition, these figures include simulation results with = 5 and = 8 relaxation times (see Section 3.2.2).
As it can be seen both the single finite element with anisotropic viscoelastic material behaviour and the RVE show nearly identical homogeneous responses for the solution with = 8 relaxation times. This result equally holds for the other analyses with separately applied strain coefficients. Altogether, only small deviations appear (see Figure 6 black versus green dashed), which are caused by the different discretization methods, which are implicit for the UMAT and driven by the evolution of viscous strains within the quasi-static analysis of the RVE (see [23]). Another source for deviations can be attributed to the interpolation of the PDBcs for the non-uniform mesh of the RVE (see [27]). For the case of = 5 relaxation times it becomes obvious that significant deviations occur for the homogeneous stress (see Figure 6 black versus red dashed), the dissipation power and the energy responses of the material model in comparison with the heterogeneous structure of the RVE. Here it is revealed that a homogenized material model consisting of = 5 Maxwell branches is not sufficient to simulate the RVE's behaviour with a matrix material that consists of the same number of relaxation times. In this case, especially the shear stress components and the dissipation energy are not reproduced accurately.
F I G U R E 7 Dissipation power (left) and energy (right) results within the test procedure for = 8 (green) and = 5 (red) effective relaxation times

CONCLUSION
The current work presents a numerical homogenization method for the microscopic dissipation energy density of anisotropic linear viscoelastic composites in order to identify effective material parameters of a corresponding GMM, which yields a homogeneous constitutive law by means of a convolution integral containing a tensorial relaxation function. For that reason, an expression for the GMM's dissipation energy is adjusted to the equivalent effective response of the microstructure, which results from a relaxation analysis under a complex deformation state. Within this first step, suitable homogeneous relaxation times for the composite are found. Secondly the homogeneous relaxation function is calculated by the evaluation of the composite's stress responses within several relaxation analyses. This enables the adjustment of the relaxation function with already known relaxation times. Thus, the remaining stiffness tensors can be found within linear optimizations. Especially the tensorial relaxation strengths scale the influence of the relaxation times in different directions. For the validation of this procedure, an endless carbon FRP with a polymer matrix exhibiting a significant viscoelastic behaviour on short time scales was chosen. At first, a successful identification of the effective relaxation times and stiffness tensors is illustrated. Secondly, the parameters are utilized within a homogeneous finite element analysis of a single finite element under a complex time-dependent strain, which is also applied to the RVE of the underlying composite. The comparison of the effective stress responses, the dissipation power and the energy density lead to nearly identical results, what highlights the applicability of the presented approach. As an extension, further studies should reveal an necessary maximum relative error for the dissipation energy fitting in order to reproduce the corresponding homogeneous material behaviour accurately. Within future work, the proposed approach could also serve as basis for different applications. This might include a realization of the approach in the frequency and Laplace-Carson domain.

A C K N O W L E D G E M E N T S
where is the Kronecker symbol. Within the infinitesimal strain theory, the material time derivativėis defined bẏ =̇⊗ , with fixed base vectors. Moreover, two isotropic fourth-order tensors are introduced as follows is the second-order identity tensor and 1 ( ) = ⋅⋅ the first main invariant of . Any part (elastic, inelastic) of the volume-specific stress power may be defined as = ⋅⋅̇. (A4) The corresponding volume-specific stress energy is then given by = ∫ dt. The expressions "density of power or energy" as well as "power-or energy density" always refer to this kind of volume-specific measures. Within the present work, the inelastic parts of these two expressions solely arises from viscous components and are also denoted with "dissipation power density" or "dissipation energy density", respectively. Furthermore, the current work is based on infinitesimal strains, what allows the additive split of the linear strain tensor in its elastic and viscoelastic parts. Finally, Table A1 gives an overview of further used symbols and indices.