A viscoelastic Mooney–Rivlin model for adhesive curing and first steps toward its calibration based on photoelasticity measurements

The transition of polymer adhesives from an initially liquid to a fully cured viscoelastic state is accompanied by three phenomenological effects, namely an increase in stiffness and viscosity in conjunction with a decrease in volume (curing shrinkage). Under consideration of these phenomena, some of us (Hossain et al. in Computational Mechanics 46:363-375, 2010) have devised a generic, viscoelastic finite strain framework for the simulation of the curing process of adhesives, which renders a thermodynamically consistent model regardless of the selected free energy density. In the present work, this generic curing framework is modified by means of more precise integration schemes and is applied to a hyperelastic Mooney–Rivlin material based on an additive volumetric-isochoric split of the strain energy density. The benefit of this decomposition is directly related to the distinct material responses of various polymers to volumetric and isochoric deformations [4]. The resulting Mooney–Rivlin curing model provides the foundation for implementing a user-defined material subroutine (UMAT) in Abaqus requiring the Cauchy stress and a non-standard formulation of the tangent operator. To this end, the corresponding transformations are presented. Additionally, a first attempt to determine the evolution of the curing-dependent material parameters through optimization with respect to a photoelasticity measurement is presented. A subset of the material properties, which reflect the emergence of shrinkage stresses inside a ceramic-epoxy composite after its fabrication, is determined via inverse parameter identification. However, due to a lack of experimental data and some rather strong assumptions made on the physics involved, this demonstration can currently be considered only as a proof-of-concept.


Introduction
In today's industry, polymer-based adhesive bonds are an indispensable way of coupling different components and materials, which are frequently utilized in various branches, e.g., in the automotive, electronics, aerospace or even medical sector. With the practical application of polymer adhesives, the demand for meaningful constitutive equations goes hand in hand with predicting the time or degree-of-cure dependent mechanical material behavior. At the molecular level, the curing process of polymer adhesives is characterized by chain polymerization and conceivably cross-linking causing substantial changes of the macroscopically observed mechanical properties. That is an increase in stiffness, relaxation time and molecular weight.

Phenomenological foundations of adhesive curing
In the first contribution of Hossain et al. [1], a one-dimensional small strain curing model is developed by successively taking into account the aforementioned effects. Indisputably, these considerations are absolutely crucial, as they provide the foundation for the corresponding three-dimensional finite strain models presented in [2] and [3]. Thus, the most important aspects are briefly summarized as follows. 1. Increase in stiffness: A macroscopic deformation of a curing solid is accompanied by an elongation of the polymer chains inbetween existing crosslinks. Hence, as polymerization progresses, new crosslinks will fit into the deformed configuration resulting in a gain in stiffness, but these new, not yet deformed crosslinks will only contribute to the stress answer if and only if a new deformation is applied. These considerations at the molecular level can now be simply transferred to a rheological model, which initially consists of a single Hookean spring sustaining a strain increment. By continuously adding linear springs in parallel to the already strained springs and consecutively applying new strain increments to the resulting spring network, ultimately a time-continuous constitutive relation for the rate of stress can be derived, cf. [1]:σ (t) = c(t)ε(t), with a stiffness evolution c(t) and the strain rate historyε(t). This constitutive equation is in line with the preliminary considerations at molecular level, as it holds thatσ (t) = 0 if ε = constant. 2. Increase in viscosity: Besides the gain in stiffness, the successive cross-linking will cause a restriction of the molecular movements leading to a macroscopic gain in viscosity. To this end, the elastically curing spring is extended by a Maxwell element in parallel, consisting of a regular spring and a damping device with time-dependent viscosity in series. The assumption of a time-invariant spring does not cause a qualitative change in the stress response due to the decay of the non-equilibrium stress response over time, as discussed in [3]. 3. Volume shrinkage: Moreover, cross-linking narrows the distance between molecules during the polymerization, which causes under the assumption of mass conservation an increase in density. In a one-dimensional setting, the curing shrinkage can be taken into account by extending the elastic-viscous configuration by a new rheological device in series, as can be seen in [1].

Previous works on modeling adhesive curing
A detailed review of phenomenologically motivated models capturing the most important mechanical properties of curing polymers can be found in our previous contributions [1][2][3]. It is noteworthy that these contributions also provided the foundation for the extension toward a magneto-mechanical curing model as devised in a series of publications by Hossain et al. [5][6][7]. Recently in [8], this purely mechanical framework has also been extended to particle-filled electro-active polymers, which is based on the considerations in [9]. In Liebl et al. [10], a 3D viscoelastic-viscoplastic material model for curing adhesives including chemical shrinkage is proposed, which is based on the one-dimensional model developed in [11].
Among others, the rapid evolution in the field of additive manufacturing is driving the development of new curing models capturing the various underlying physical effects of the related manufacturing process. In view of stereolithography and digital light processing, both belonging to the type of light polymerization techniques, a curing model, which takes into account the light intensity as well as the temperature of the manufacturing process, is developed by Rehbein et al. [12]. For the highly dynamic process of extrusion-based additive manufacturing, a curing model for soft polymers capable of reproducing the initial fluid-like behavior, which requires a mesh-free solution, is devised by Hartmann et al. [13]. Note, this excerpt is by no means complete but is rather intended to picture the new upcoming application fields requiring reliable and innovative polymer curing models.

Photoelasticity and shrinkage stresses
Probably, the most prominent phenomenon arising in polymer adhesives is the curing shrinkage effect, which may lead to significant stresses and distortions, even without external loading, especially in the case of composite components. A widespread whole-field optical experimental method with a long history in engineering science suitable for visualizing and determining such shrinkage stresses is a photoelasticity measurement. Frequently, this methodology is utilized to experimentally determine the shrinkage stresses, as done, e.g., in the field of dental implantology in [14,15].
After its discovery, photoelasticity became an indispensable tool for determining stress distributions in complex or large objects, by analyzing transparent small-scale replicas. Although temporarily overshadowed by numerical methods like FEM, the importance of experimental stress analysis by means of photoelasticity picked up again in recent years [16]. Paradoxically, among other things this can be attributed to the significant increase in computing power, breaking the ground for digital photoelasticity [17], which promoted the disappearance of photoelasticity in the first place. Additionally, the advent of stereolithography has sparked renewed interest, which enables to prototype complex transparent models rather in hours than days.
Nowadays, stereolithography in conjunction with three-dimensional photoelasticity via stress freezing and slicing provides an accessible experimental stress analysis method for complex geometries, whereby the experimental stress data can be used to validate corresponding FEM models, as, e.g., done for a stressed bearing in [18].

Aims and outline
The present contribution aims at a fusion of the above-mentioned aspects, i.e., (1) to measure shrinkage stresses with photoelasticity, (2) to simulate curing with a new, viscoelastic and isothermal finite strain model, (3) to calibrate this model with data from (1), i.e., to determine the evolutions of the material parameters during curing from an inverse optimization procedure.
The resulting approach is beneficial since a large number of parameters and their temporal evolutions are required to simulate curing adhesives. Note that our approach still presupposes a subset of the material parameters to be known from other experiments. The details of the optimization approach are discussed thoroughly in Section 5. Besides that, this paper is organized as follows: Section 2 summarizes the generic finite curing model as proposed in [2,3]. Based on this framework, the corresponding decoupled Mooney-Rivlin curing model is derived in Section 3. In view of the optimization procedure, the theory of photoelasticity as well as its application to a ceramic-epoxy composite are presented in Section 4.

Generic viscoelastic finite strain curing framework
The true mechanical energy density m accumulated from deformations and stiffness increase during the curing process adopts the frequently employed approach of finite viscoelasticity modeling, in which the traditional mechanical strain energy density is additively composed of an (elastic) equilibrium and a (viscous) nonequilibrium part 1 : where eq m denotes the equilibrium energy related to the increase in stiffness and neq m the non-equilibrium energy associated with the increase in viscosity. The particular format of both energy contributions will be clarified in detail in Sections 2.1 and 2.2, respectively. A superposition like in Eq. (1) is motivated by kinematical considerations, mostly from a multiplicative decomposition of the mechanical deformation gradient F m and was firstly proposed by Lion [19] for the curing of polymers including shrinkage. Contrary to the approach in Lion [19], we do not provide an additional chemical potential for Eq. (1) to drive the curing shrinkage, but incorporate the shrinkage directly by scaling the mechanical deformation F m according to a prescribed shrinkage coordinate s := s(t) ≤ 0. This is motivated by the fact, that, unlike a chemical potential itself, the effect of volume shrinkage, and thus, the shrinkage coordinate s(t) can be straightforwardly determined experimentally.
To this end, the deformation gradient is decomposed according to Fig. 1 into 1 Remark: Hereafter, we distinguish two types of potentials or strain energy densities, namely traditional and true, denoted by and , respectively. A traditional potential could be any hyperelastic law, e.g., the Neo-Hookean model, that is equipped with time-dependent Lamé parameters to cover the increase in stiffness during curing. The stress-strain relation following from such an enhanced traditional potential does not satisfy the curing specific properties as discussed in Subsection 1.1, i.e., the stress would increase permanently, even when the strain rate is zero. To resolve this problem, a true accumulated potential requires a particular format like given in Eq. (17), but is still based on a traditional with evolving material parameters via the stiffness operator included. This notation coincides with that used in our previous contribution [2], in which a more in-depth explanation can be found.  1 Multiplicative decomposition of the deformation gradient F, which maps the tangent space of material configuration 0 to its spatial counterpart t , i.e., F : T 0 → T t , into a shrinkage F s and a mechanical part F m , which in turn is divided into a viscous F v and an elastic contribution F e , by adding two stress-free intermediate configurations: a shrinkage¯ s and a viscous v one, adapted from Belytschko et al. [32] where the shrinkage deformation gradient F s depends on the aforementioned shrinkage coordinate s. The evolution of s is described by an exponential saturation function, which can be expressed in a generic format as follows: with the initial value χ 0 at the gel point, the final value χ ∞ and the curvature parameter β χ . By substituting s for χ in Eq. (3) and noting that the volume shrinkage is negligible before the gel-point, i.e., s 0 := 0, the volume shrinkage evolution simplifies to . The evolutions of all other curing-dependent parameters are analogously derived from Eq. (3) and are introduced in the remainder of this section. The choice of a multiplicative decomposition like Eq. (2), cf. e.g., [20][21][22], usually provokes the question for an alternative additive decomposition, cf. e.g., [23][24][25]. This is a long-standing debate that has led to controversies and comparative studies, e.g., [26][27][28][29], in the material modeling community. There are compelling arguments in favor of both options, e.g., the physical appeal of the multiplicative decomposition in single crystal plasticity or the numerical ease of additive decompositions, cf. e.g., [30] and the references therein. However, there is no consensus on the single-one best option. In the context of isotropic curing, nevertheless, the multiplicative decomposition separating the isotropic shrinkage induced contribution to the deformation gradient is following the well-accepted route as also pursued in thermo-elasticity. This leads to no particular numerical obstacles, while, at the same time, fitting nicely into a sound geometry-inspired approach. Thus, the multiplicative decomposition has been favored here for the finite strain curing framework.
Unlike the shrinkage and mechanical deformation gradients F s and F m , neither the viscous F v nor the elastic part F e are prescribed or computed directly, as will be discussed in detail in Section 2.2. Nonetheless, the multiplicative decomposition of the deformation gradient from Eq. (2) enables the definition of mechanical elastic and viscous right Cauchy-Green tensors analogous to the right Cauchy-Green tensor C = F T · F as follows: which in turn are utilized to apply the same functional relationship for the respective shares of the total accumulated energy density in Eq. (1), but in dependence of different strain measures, as is usual with phenomenological models: To derive constitutive equations compatible with the second law of thermodynamics the curing model has to satisfy the Clausius-Duhem inequality requiring the dissipation power density D to be non-negative for all considered processes. By restricting ourselves to the isothermal case, the dissipation inequality boils down to: with the stress power density P, which equals the double contraction of the Piola-Kirchhoff stress tensor S with the material rate (•) of a work-conjugate strain measure, either Green-Lagrange strain E = 1 2 C − I or right Cauchy-Green tensor C. In order to evaluate dissipation inequality Eq. (6), the material time derivatives of the free energy contributions from Eq. (1) are computed using the chain rule, which yields: with the reduced viscous dissipation term Inserting Eq. (7 -8) into Eq. (6) renders the Clausius-Duhem inequality for the generic finite strain curing model: To ensure that the dissipation inequality (10) is fulfilled for all processes, we follow the standard Coleman-Noll procedure [31], which yields the following conditional equation for the Piola-Kirchhoff stress: as well as a non-negativity constraint for the reduced viscous term (9): which is fulfilled by means of a compatible evolution equation for C v as will be demonstrated in Section 2.2. For the successful implementation of a constitutive model into an FEM package, the so-called tangent operator is required to set up the system stiffness matrix in every iteration, which emerges from Eq. (11) as follows:

Elastic equilibrium contribution
As discussed in Section 1, the increase in stiffness during curing can be taken into account in the onedimensional case by the following rate-dependent constitutive equation valid for small strains: whereby a prescribed stiffness evolution c(t) relates the strain rate historyε(t) to the stress rateσ (t). A naive but expedient approach for the extension toward large strains and three dimensions is to simply transfer Eq. (14) into the following hypoelastic constitutive equation for the mechanical Piola-Kirchhoff stress tensor in the material configuration 0 :Ṡ eq In contrast to the one-dimensional setting, the stiffness operator T eq m (t) is no longer prescribed directly but is rather derived from an arbitrary mechanical equilibrium strain energy density eq m (C m , ϑ eq ), which is dependent on C m and governed by a set of r evolving, i.e., curing-dependent material parameters ϑ eq := {ϑ eq 1 , . . . , ϑ eq r }. In this framework, each ϑ eq i is either explicitly or implicitly characterized by means of an evolution equation like Eq. (3). Therefore, eq m is not only implicitly time dependent due to the external loading C m , but also explicitly due to the curing progress reflected by the evolving material parameters ϑ eq . For ease of notation, only this explicit time dependency of eq m is indicated for further elaborations. Hence, the derivation of the stiffness operator can be written as follows: Equation (15) has not yet been proven to be thermodynamically consistent. By applying the standard Coleman-Noll procedure to the following convolution integral: characterizing the true accumulated mechanical energy density, it can be shown that constitutive Eq. (15) is recovered, i.e., that eq m (t) is the potential required to make Eq. (15) a thermodynamically consistent relation. In Eq. (17), T eq m = dT eq m (v)/dv denotes the total differential of stiffness operator (16) with respect to the integration variable v. The thermodynamical consistency of this stiffness operator is intrinsically given since it is derived from a traditional potential . For a detailed derivation and discussion as well as a physical interpretation of the convolution integral, we refer to [2].
Due to its beneficial properties, which are unconditional stability and second order accuracy, the midpoint rule is chosen for the numerical integration of tensor-valued differential equation (15), yielding the Piola-Kirchhoff stress at time t n+1 : 1 2 , 1} abbreviates the evaluation at time t n+k . In accordance with this numerical integration of the stress, the corresponding consistent algorithmic tangent operator is required to set up the system stiffness matrix such that a quadratic convergence of the underlying local Newton-Raphson scheme is preserved [33]. Note that we need to carefully distinct between the elastic stiffness operator T eq m , which follows from the underlying potential eq m , and the consistent algorithmic tangent operator E eq m , which relates the changes in stress to the changes in strain at time t n+1 and can be derived from Eq. (18) as follows: Therein, the sixth-order tensor A eq m = ∂T eq m /∂C m appears, which is the derivative of the current, curingdependent material stiffness with respect to the current strain. Note that E eq m is utilized for Eq. (13) due to the aforementioned reasons.

Viscoelastic non-equilibrium contribution
Contrary to the equilibrium part, the true accumulated mechanical non-equilibrium energy density is simply identified with a traditional strain energy density, i.e., in which only the relaxation time } is assumed to be curing-dependent and specified by means of Eq. (3), i.e., τ = τ (τ 0 , τ ∞ , β τ , t). All other material parameters in neq m are assumed to be independent of the curing progress, following the same argumentation as for the time-independent spring in the Maxwell element in the one-dimensional small strain setting, cf. Section 1.1. Due to this simplification, the non-equilibrium stress is straightforwardly obtained by substituting Eq. (20) into Eq. (11), followed by an evaluation at time t n+1 : The same applies for the non-equilibrium tangent operator, which follows from inserting S neq m as in Eq. (21) into Eq. (13), again evaluated at time t n+1 : Note that C e , which presumably remains after the differentiation in Eq. (21), cannot be computed explicitly by means of Eq. (4) due to the absence of an equation for F v . However, this can be obviated by rewriting S neq m in terms of the mechanical C m and viscous right Cauchy-Green tensor C v , whereby the latter represents the strainlike internal variable, associated with the damping device in a Maxwell element in a one-dimensional setting. Irrespectively whether the viscous right Cauchy-Green tensor C v is known, the corresponding deformation gradient cannot be recovered due to the lack of rotation information in the strain tensor. In order to obtain a sufficient evolution equation for C v not violating the second law of thermodynamics, C v has to satisfy D red v ≥ 0. Inserting the chosen ansatz of Eq. (20) into Eq. (12) yields: After the application of two consecutive linearization steps to Eq. (23), a simple ordinary differential equation These two steps are as follows: (1) Linearization around the thermodynamic equilibrium state renders an evolution equation still valid for large deformations, but restricted to only small deviations around the equilibrium state [34]. (2) Consecutive linearization of C v obtained via (1) around the material configuration then results in Eq. (24), with an additional limitation to small deformations [35].
To preserve the benefits of the midpoint-rule used for integrating the equilibrium stress in Section 2.1, an integration scheme with the same properties is chosen for Eq. (23), that is the trapezoidal rule, which yields: The above derived generic equations for stress tensors and tangent operators of a curing material are now specified for a particular constitutive model, which is frequently used for polymers.

Viscoelastic Mooney-Rivlin curing model
In Hossain et al. [3], the previously presented generic finite strain curing model has successfully been applied to the compressible Neo-Hookean solid, based on an indistinguishable coupled ansatz. This approach yields a curing model captivating by its compact formulation and low computational costs, but lacks in view of determining the material parameters through purely distortional and dilatational experiments [4] and in accuracy at large strains [36]. Another practical drawback is attributed to the fact that commercial FE packages predominantly rely on decoupled hyperelastic solids, cf. Steinmann et al. [37] for an overview of selected representatives. Hence, the equivalence to the built-in hyperelastic material models in the limits of the curing model vanishes. These disadvantages with acceptance of higher computational costs may be overcome by considering a compressible Mooney-Rivlin solid based on an additive volumetric-isochoric energy decomposition as given by: in whichĪ 1 = trC andĪ 2 = 1 2 [tr 2C − trC 2 ] denote the first and second principal invariant of the isochoric right Cauchy-Green tensorC =F T ·F, which in turn arises from the multiplicative decomposition of the deformation gradient into F = J 1 3F . Contrary to isochoric part iso with parameters c 1 and c 2 , volumetric contribution vol depends on the whole deformation gradient in terms of its Jacobian J = det(F), with corresponding volumetric parameter d 1 . Due to this decomposition, a volumetric deformation strictly belongs to vol , whereas an isochoric deformation solely invokes iso . This can be easily verified by showing that the isochoric deformation gradient does not cause a change in volume, i.e., det(F) = 1.
As a direct consequence of additive energy decomposition (26), the stress tensor and the tangent operator can likewise be divided into isochoric and volumetric parts:

Elastic Mooney-Rivlin equilibrium part
Following the discussion in Section 2.1, the true accumulated mechanical equilibrium energy density depends on a traditional energy density with time-dependent material parameters ϑ eq := ϑ eq (t), i.e., eq m := eq m eq m (C m , ϑ eq ) . To this end, Mooney-Rivlin energy function (26) is equipped with the curingdependent parameters ϑ eq = {c 1 (t), c 2 (t), d 1 (t)} and rewritten in terms of C m : in which (•) t := •(t) and (•) = (•) m are utilized to highlight the curing dependency of the material parameters and the dependency on C m , respectively. The isochoric invariants have been replaced by the traditional principal invariants given in Appendix A, which also includes some relations from tensor calculus required in the following. By taking advantage of Eq. (28), the elastic stiffness operator, which is required to update equilibrium stress (18) and tangent operator (19), can be derived as follows: With the components of the Mooney-Rivlin energy density (28), the corresponding volumetric and isochoric parts of the elastic stiffness operator then follow as: whereby some abbreviations have been introduced for the sake of clarity: The different non-uniform dyadic products used therein are given by the following definitions for two arbitrary second-order tensors A and B: The consistent algorithmic tangent operator n+1 E eq m from Eq. (19) additionally requires the computation of the sixth-order tensor n+ 1 2 A, i.e., the partial derivative of (30) with respect to C m , which has been outsourced to Appendix B for the sake of readability.

Viscoelastic Mooney-Rivlin non-equilibrium part
The application of ansatz (20) to Mooney-Rivlin model (26) provides the following true accumulated mechanical non-equilibrium energy density: wherein ϑ neq = {c 1 , c 2 , d 1 , τ (t)} summarizes the non-equilibrium material parameters and (•) e = (•)(C e ) highlights the dependency of the invariants on only the elastic right Cauchy-Green tensor C e . The elastic Jacobian J e = J m J −1 v follows as the quotient of the overall mechanical Jacobian and its viscous counterpart. The only curing-dependent material parameter, the relaxation time τ (t), is not immediately apparent from the energy density (33), but appears implicitly in the non-equilibrium stress via the evolution equation (24) for the internal variable C v , the viscous right Cauchy-Green tensor. Due to the separate isochoric and volumetric contributions in Eq. (33), the non-equilibrium stress (21) can be subdivided analogously: and some rigorous derivations then provide the individual stress contributions as: with the abbreviation K : eq m , the non-equilibrium Mooney-Rivlin tangent operator required in Eq. (13) has been outsourced to Appendix B for the sake of readability.

Photoelasticity measurements
Photoelasticity is a non-destructive optical testing method to determine stress distributions inside transparent specimens. In a plane-stress setting, the initial optical isotropy of the specimen is disturbed by the applied load, which results in two different refraction indices n 1 and n 2 at each point of the specimen. This essential property, a load-induced optical anisotropy, is denoted as stress birefringence. An experimental setup of a twodimensional plane-stress photoelasticity measurement in a dark-field configuration using linearly polarized light is depicted in Fig. 2. The notion dark-field stems from an angle of π 2 between polarizer and analyzer [38]. After leaving the polarizer, the linearly polarized light represented by amplitude vector A is split into two orthogonal parts A 1 ⊥ A 2 according to the specimen's inherent principal stress directions σ 1 ⊥ σ 2 , cf. Fig. 2. Due to the unequal refraction indices, the two waves propagate at different velocities through the specimen, causing a phase shift , which emerges from the stress-optic law as follows [39,40]: Therein, t denotes the thickness of the specimen and C is referred to as stress-optical constant, which has to be determined from experiments for each material. Subsequently, the phase shifted vectors A 1 and A 2 impinge the analyzer, which only allows the respective horizontal components H 1 and H 2 to pass.This does not alter the phase shift , but results in identical amplitudes |H 1 | = |H 2 | = |A| sin α cos α, which follows from trivial geometric considerations. Here, α denotes the angle between the horizontal axis and the axis of the first principal stress σ 1 , cf. Fig. 2. By taking into account , the amplitude of the resulting light vector emerges from the superposition of the two shifted horizontal light vectors H 1 and H 2 as follows: where δ = λ with observation wave length λ denotes the relative phase shift or fringe order. A more detailed derivation can be found, for example, in [38][39][40][41]. According to Eq. (37), complete extinction, i.e., a dark point, occurs if either sin 2α or sin πδ equals zero, which also applies to the intensity I , since I ∝ A 2 . The former holds if the principal stress direction coincides with the polarization direction, i.e., α ∈ {0, π 2 , π, 3π 2 , . . . } and thus, the principal stress trajectory at any point of the specimen can be obtained by rotating the specimen in Fig. 2. The resulting cancellation fringes are referred to as isoclinics. The latter holds if δ ∈ N 0 . Conversely, maximum brightness is obtained in between the fringe orders, i.e., δ + 1 2 . Both yield so-called isochromatics, which are lines of equal principal stress differences remaining unchanged through rotation. Frequently, only the magnitude of the stresses is of particular interest. To this end, the isoclinic fringe pattern can be suppressed by inserting quarter wave-plates on both sides of the specimen, which renders a circularly polarized light. As a consequence of that, the amplitude from Eq. (37) simplifies to: A typical example for an isochromatic pattern obtained by a circular dark-field polariscope is depicted in Fig. 3 (left), which represents the stress distribution inside a ceramic-epoxy composite [42]. The composite consists of regularly distributed ceramics cubes (black), which are connected by an epoxy resin (gray). To arrive at the Experimental setup for a two-dimensional plane-stress photo-elasticity measurement in dark-field configuration. Due to the stress-induced birefringence, linearly polarized light (amplitude vector A) passing through the specimen is decomposed into two orthogonal parts A 1 ⊥ A 2 according to the specimen's inherent principal stress directions σ 1 ⊥ σ 2 , finally leading to a phase shift . quantified pattern in Fig. 3 (right), δ is evaluated point-wise with the experimental setup in Fig. 2, whereby a compensator was additionally utilized. In connection with the stress-optical constant of the epoxy, which was determined experimentally as C = [26.8 ± 0.9] · 10 −12 mm 2 N , σ P E can be calculated from the stress-optic law (36). A detailed description of the presented experiments and results can be found in [42]. Due to the absence of external loadings, the clearly visible stresses in the cured epoxy layers connecting the ceramic blocks solely arise from alterations of the physical state of the composite's constituents as caused by the manufacturing process, which consists of four stages: (1) Infiltration of the ceramics cubes grid with liquid epoxy.
Due to the temperature changes involved and the mismatch of the thermal expansion coefficients of epoxy and ceramic blocks, • ζ epo = 76 · 10 −6 K −1 , • ζ cer = 6.8 · 10 −6 K −1 , both the curing shrinkage of the epoxy and thermal strains contribute to the residual stresses visible in Fig. 3. To analyze this interplay, we start at the initial phase of stage (3). As a consequence of the temperature increase from stage (2) to stage (3), the epoxy experiences compression, since it expands much more than the surrounding ceramic blocks. At the same time, however, it also experiences tension due to the curing shrinkage. Both deformations cannot evolve freely but are hindered by the much stiffer ceramic blocks, which leads to internal stresses in the epoxy. These stresses are in turn continuously relaxing, at least partly, according to the current curing-dependent viscosity of the cross-linking epoxy. During the cooling in stage (4), the resultant strain state from stage (3) is again overlaid by further volumetric contractions. The final stress distribution depicted in Fig. 3 hence corresponds to a superposition of tensile and compressive strains from curing shrinkage and thermal expansions and contractions, potentially followed by further relaxations after room temperature has been reached everywhere in the sample. Unfortunately, the exact time of the photoelasticity measurement after stage (4) and the particular strain states after stages (2) and (3) are beyond the authors' knowledge. Consequently, the respective contributions of curing shrinkage and thermal strains on the overall stress distribution cannot be separated retrospectively, but would only be accessible from two more photoelasticity measurements after stages (2) and (3). Due to this lack of information, we assume hereafter for simplicity that the stresses in Fig. 3 are exclusively caused by the curing shrinkage, which corresponds to a kind of worst-case scenario, since the reducing effect of the thermal strains is neglected. The identification of the curing parameters described below is therefore not quantitatively exact, but nevertheless gives an idea of the potential of the method proposed here.

Optimization of curing model parameters
The above derivation of the Mooney-Rivlin curing model revealed that quite a large set of material parameters ϑ is required to fully describe the curing process of polymeric adhesives. The experimental identification of the whole set ϑ may pose a difficult task, but is, of course, indispensable to achieve quantitative agreement between numerical simulations and reality. Still, one can assume that a certain subset ϑ f ix ⊂ ϑ can easily be determined from experiments, like the moduli of the uncured and fully cured epoxy. More difficult to measure are, for example, the curvature coefficients of the exponential saturation functions for the curing evolutions of the material parameters or the overall curing shrinkage. The optimization approach subsequently discussed hence concentrates on the determination of this remaining set of curing parameters ϑ opt with ϑ f ix ∪ ϑ opt = ϑ.
The application of a FEM-based inverse method is a popular choice for the characterization of Mooney-Rivlin material parameters in a non-curing setting, which has been successfully applied using for instance biaxial [43] or nanoindentation tests [44,45]. The optimization procedure used in the following is similar and relies on the stress distribution inside the cured epoxy around the ceramic blocks in Fig. 3, as measured with photoelasticity. Due to the already discussed lack of experimental information and for the sake of simplicity, we assume that this stress distribution has emerged solely from the curing shrinkage after a total curing time of 60 seconds, like in a UV-curing epoxy. The shrinkage stresses thus coincide with the photoelastic stresses, which serve as target values for the downstream parameter optimization. Discrete target values for the goal function are determined from sampling the relevant area of the stress distribution in Fig. 3 (right), like in Fig. 4. We now mimic the curing process of the epoxy around the ceramic blocks and the subsequent photoelasticity measurement by means of FE-simulations, which apply the viscoelastic Mooney-Rivlin curing model and its corresponding parameters ϑ. The unknown parameters ϑ opt are thereby iteratively adapted by a gradient descent least squares optimization, until the simulated photoelastic stress distribution matches the measured one from Fig. 4 within a certain tolerance. Taking advantage of the periodic structure of the composite, it is sufficient to simulate only a segment of the overall structure, as illustrated in Fig. 5. All dimensions are chosen such that the real geometry, cf. Fig. 3 (left), is reproduced. Furthermore, the analysis area ( ), which is then compared to Fig. 4, is sufficiently far away from the two front surfaces. The two back surfaces are fixed such that symmetry is realized, i.e., they cannot move in their respective normal directions. At the bottom surface, the ceramic blocks are fixed in vertical direction, but lateral movements are allowed. To avoid delamination or slippage between ceramic blocks ( ) and epoxy ( ), their surfaces are tied to each other.
Concerning the constitutive laws, the ceramic blocks are modeled as linear elastic, with the Abaqus internal implementation, a Young's modulus of 390 GPa and a Poisson's ratio of 0.3. For the epoxy, the above-derived viscoelastic Mooney-Rivlin curing model is used, which has been implemented in Abaqus in terms of a Fortran user material subroutine (UMAT ). Thereby, particular attention must be paid to the element type used for the discretization. The curing model is formulated in terms of the conventional deformation gradient, which is passed to the UMAT by Abaqus if and only if second-order elements are used. When using linear elements, Abaqus takes advantage of the B-bar method and accordingly modifies the deformation gradient passed to the UMAT, which leads to unacceptable results. A general guideline for the implementation of a UMAT can be found in [47], and Appendix C provides a detailed discussion about the correct formulations for the required stress tensor and tangent operator.
The 2D photoelastic stress distribution in Fig. 4 is a result of the phase shift accumulated by the light during its way through the specimen. To enable the comparison of this measured photoelastic stress and the FE-simulated 3D stress state caused by the curing shrinkage, which is required by the optimization procedure, the latter has to be post-processed to arrive at σ P E at each grid-point, as e.g., shown in [48]. To simplify the post-processing, the discretization of the analysis area ( ) is chosen such that the projection of the nodes onto the top x y-plane matches the grid-points. Thus, σ P E of each grid-point can be straightforwardly computed by means of the corresponding nodal path, which coincides to the light path through the specimen (with thickness t) in the undeformed configuration, as depicted in Fig. 6. Under the assumption that the analysis area ( ) is straightly transilluminated in z-direction, only the stress components in the x y-planes have to be considered [41]. Photoelastic stress σ P E at nodal position z k can hence be rewritten in terms of the in-plane stress values provided by Abaqus, by applying the following principal stress transformation: In general, the principal stresses and their directions α i vary slightly throughout the nodal path, cf. Fig. 6. By assuming the orientation to be constant, the corresponding photoelastic stress can be derived by means of an integral over the thickness of the specimen, according to [48]: whereby n is denoting the number of elements passed by the light in z-direction, cf. Fig. 6. The integrals in (40) can further be approximated using the trapezoidal rule with the nodal values σ P E z k of the photoelastic stress: wherein t k+1 denotes the thickness of the related element, cf. Fig. 6. The application of the trapezoidal rule is not restricted to the second-order elements, but is also compatible with linear elements, which does not hold for the midpoint rule, that has been used for the integration of the curing stress in Eq. (18). The combination of Eqs. (39)(40)(41) finally yields the numerical counterpart to the experimental photoelastic stress at each grid-point. The practical realization of the identification of the curing material parameters through optimization works as follows. In each iteration loop, the curing process is simulated with the Abaqus FE epoxy-composite model (Fig. 5), based on the finite Mooney-Rivlin curing UMAT and current material parameter set ϑ. A Python script then extracts the resulting photoelastic stress distribution according to (39)(40)(41) and passes it to an optimization routine available in MATLAB (lsqcurvefit, a gradient descent method optimizer with finite differences for the computation of the sensitivities), which iteratively updates the parameter subset ϑ opt by means of the least squares difference between experimental and current photoelasticity stresses. These steps are repeated until a prescribed tolerance between simulated and experimental σ P E is underrun.
Fundamental to the proposed procedure is that a subset of material parameters ϑ f ix has been determined experimentally prior to the optimization. Material testing generally provides traditional parameters like Young's modulus or Poisson's ratio, rather than the separate Mooney-Rivlin parameters for theĪ 1 -andĪ 2 -terms in (26). We therefore begin with total values ϑ eq f ix and ϑ neq f ix for the overall equilibrium and non-equilibrium contributions, each in terms of traditional material parameters, which are derived from the values reported in [42] for the (un)cured epoxy and summarized in Table 1. Thereby, μ(t) and ν(t) denote the curing-dependent shear modulus and Poisson's ratio of the equilibrium part, both known at gel-point t 0 and at fully-cured state Table 1 Fixed material parameters ϑ f ix = ϑ eq f ix ∪ ϑ neq f ix , based on [42] 0.1 7200 s The same applies for the relaxation time τ (t) of the non-equilibrium part. The chosen initial and final values of ν(t) reflect the transition of the epoxy from a nearly incompressible fluid to a compressible solid, whereas the non-equilibrium Poisson's ratio ν e is set to a constant and nearly incompressible value due to the already discussed invariance of the Maxwell element stiffness. Contrary to these equilibrium and non-equilibrium parameters, all shrinkage parameters summarized in ϑ s are to be optimized, i.e., ϑ s ⊂ ϑ opt . In order to preserve consistency with linear elasticity for small strains, the following relations between traditional elasticity and Mooney-Rivlin parameters have to hold ∀t: The volumetric parameter d 1 (t) hence is computed by means of (42) 2 , rather than prescribing an independent saturation function like (3). Since the distribution of the shear modulus μ(t) to the isochoric parameters c 1 (t) and c 2 (t) is not known, a linear dependency is a priori assumed for simplicity: i.e., the quotient r c is introduced as another variable to be optimized. Equations (42) and (43)  Since no information about the stiffness of the non-equilibrium part is available in terms of experimental data, another linear relation is introduced: i.e., r m is to be optimized and relates the first isochoric parameter of the equilibrium with that of the nonequilibrium part. The second isochoric non-equilibrium parameter again follows from (43) as c e 2 = r c c e 1 . Of particular interest are the curvature parameters β i governing the velocity of the exponential saturations (3) of the individual parameters, which are difficult and expensive to identify experimentally. However, approximately identical values for the individual curvature parameters can be expected for simplicity, which motivates the definition of a generalized curvature parameter: which, together with the overall shrinkage, implies an optimization set with only four parameters: Initial values and definition ranges chosen for each parameter to ensure physically reasonable behavior are summarized in Tab. 2, which also states the final values after optimization.  Table 3 and Eq. (3) Table 3 Optimal Mooney-Rivlin curing parameters according to Table 2 with β i = 0.03617 s −1 The overall curing shrinkage s(t ∞ ) ≈ 0.6% meets the expectations, as values of up to 2% are reported in the literature. Moreover, the optimization provides a generalized curvature parameter β ≈ 0.036 s −1 , which means that all curing-dependent parameters following exponential saturation attain approximately 88.6% of their maximum values after the total curing duration of 60s. The relaxation time τ (t) rapidly exceeds the total curing time (after ≈ 0.24s), i.e., the non-equilibrium stiffness contributes completely to the overall response beyond this point. The optimized values are comparatively insensitive to the initial values since only four variables need to be determined, cf. Eq. (46), due to the assumptions and dependent parametrizations made before. Therefore, the optimization manifold is low dimensional, and the gradient descent procedure is well posed. The complete set of the optimized Mooney-Rivlin parameters following from r c and r m is summarized in Table 3. Figure 7 visualizes the exponential saturations which are followed during curing by the stiffness parameters, the relaxation time and the shrinkage parameter. A comparison of optimized (Fig. 8) and goal (Fig. 4) distributions of the photoelastic stress shows a convincing agreement. The stress decay toward the center of the analysis area is qualitatively and quantitatively reproduced. Likewise, the large stresses at the corners of the analysis area are also evident. However, the stress state at the left interface is not yet perfectly recovered. This may primarily be attributed to the manual quantification procedure  of the experimental stress distribution without access to the measured values. Although the discretized stresses reflect the experimental picture very well, at least visually, it is rather unlikely that the gradients between adjacent grid points fully coincide with the true values. Thus, the reconstructed stress distribution represents a goal function that can be reproduced by the FE-Model only up to a certain degree. This can also be seen in Fig. 9, which compares one-dimensional versions of the discretized and optimized nodal σ P E values. The optimized curve is much smoother in between neighboring grid points than the goal curve and provides a more physical profile. Significantly smaller errors are hence expected if more accurate experimental data were available.
Finally, Fig. 10 shows how solely the curing shrinkage and the increasing material stiffness lead to the emergence of a von Mises stress distribution inside the epoxy. These shrinkage stresses obviously can reach significant magnitudes, especially at the interfaces with the ceramics blocks and in the corners, even without any external loadings present. This highlights the importance of a decent curing model, particularly for the design and dimensioning of adhesive bonds, which obviously can be subject to significant preloads, simply by virtue of their manufacturing process.

Summary and conclusion
The presented extension of the coupled Neo-Hookean finite curing model of Hossain et al. [1][2][3] towards a decoupled Mooney-Rivlin finite curing model is accompanied by an augmentation of the curing-dependent material parameters to be identified. Although the decomposition into strictly volumetric and isochoric con-tributions facilitates or at all enables the experimental determination of related material parameters, indeterminate parameters may still remain. To bridge this gap, the presented Mooney-Rivlin curing model is applied to photoelasticity measurements within the scope of an optimization procedure to identify remaining curing parameters. A first methodological investigation revealed that the proposed method appears to be technically feasible as the numerically identified parameters meet the expectations. However, the derived material parameters have limited significance since the optimization goal is detached from measurement raw data. Based on the presented proof-of-concept, the determination of more quantitative photoelasticity data and its comprehensive application for the identification of material parameter evolutions during curing are planned for the near future.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding Open Access funding enabled and organized by Projekt DEAL.

Declaration
Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest.

Appendix A: Some basic tensor relations required
Principal invariants and their partial derivatives: For second-order tensors A: For symmetric second-order tensors B, i.e., B = B T :

Appendix B: Auxiliary tensors for the Mooney-Rivlin curing model
Sixth-order tensor A to update equilibrium tangent operator (19):   For a successful implementation of a UMAT, two tensors have to be computed at time step t n+1 and must be returned to Abaqus: the tangent operator n+1 t abq , which is crucial for the assembly of the tangent stiffness matrix and the Cauchy stress n+1 σ , cf. [49]. The tangent stiffness matrix results from the linearization of the nonlinear finite element equations required for the solution within an iterative scheme such as the Newton-Raphson method. This section is intended to provide the fundamental theoretical background to implement the Mooney-Rivlin curing model in Abaqus, i.e., to provide the transformations necessary for stress and tangent, which have been derived in the material configuration so far. We do not aim at a comprehensive discussion of all aspects on practical implementation and theoretical background of an UMAT, and more details can be found in [46,47,[49][50][51].
Appendix C.1: Cauchy stress The updated Cauchy stress is simply obtained from a push-forward operation of the Piola-Kirchhoff stress in Eq. (11) to the spatial configuration: by means of the complete deformation gradient n+1 F and its corresponding Jacobian n+1 J = det( n+1 F). It thereby is important to follow the discussion on element types in Sec. 5 to ensure that Abaqus provides the intended deformation gradient.

Appendix C.2: Tangent operator
The particular format of the tangent operator required by Abaqus is given by the following relation: [46,47,49,50]: i.e., the tangent J e abq relates the Jaumann rate (•) ∇ J of the Kirchhoff stress τ with the rate of deformation d, which is defined as the symmetric part of the spatial velocity gradient d = l sym . The Jaumann rate belongs to the class of objective stress rates such as the Truesdell rate and the Green-Naghdi rate, which are frequently employed in different contexts of commercial finite element software [32]. To end up with a closed relation for e abq , we begin with the definitions of the Truesdell rate (•) ∇T R and the Jaumann rate (•) ∇ J of τ : in which w is the skew-symmetric contribution of the spatial velocity gradient, i.e., w = l skw . By additively decomposing l into its symmetric and skew-symmetric parts, the following relation between the two objective stress rates can be derived: It is apparent from (50) that both rates coincide for rigid body rotations, i.e., for d = 0. Since both stress rates depend onτ , we compute the material time derivative of the Kirchhoff stress to arrive at the following relation: in which E is specified by means of (13). By comparing (49a) and (51), the Truesdell rate of the Kirchhoff stress is identified. Thus, by replacing τ ∇T R in (51) and taking into account that d is the push-forward ofĖ, we obtain: in which the [ ]-bracketed term is just the tangent relating the rate of deformation to the Truesdell rate of the Kirchhoff stress. Inserting (52) into (50) results in: In view of the definition (48) of e abq , the rate of deformation (53) must be factored out. However, since a double contraction with a symmetric tensor like d is not unique, there exists an infinite number of feasible fourth-order tensors. One possibility, which is usually employed, stems from the requirement for the tangent to preserve major and minor symmetry [51]. This approach provides the required formulation for the Abaqus tangent, evaluated at time t n+1 :