Experimental analysis and modelling of temperature- and humidity-controlled curing

There has been much discussion about modelling the reaction kinetics of a curing polymer. Typically, curing is described by the development of a variable called degree of curing as a function of temperature and time. The material considered in this paper exhibits two different curing mechanisms, namely temperature-activated and diffusion-based. To be able to describe the complex hardening process, the material is extensively analysed experimentally, and a thermodynamically consistent coupled reaction kinetics model is formulated based on experimental observations. This model enables the implementation of the thermal, caloric, and mechanical properties of the material into a finite element (FE) framework.


Introduction
In general, curing mechanisms of adhesives can occur on completely different time-scales. While almost complete curing by temperature activation is possible within a few seconds to minutes, it may take hours to days (depending on the layer thickness) to cure an adhesive by humidity. In technical applications, moisture-curing diffusion-controlled polyurethane adhesives are only utilised where curing is possible over a longer period within the production process as stated by Cognard (2006). In other cases, the acceleration of the humid curing reaction in the polyurethane adhesives via so-called booster systems is used, e.g., Cognard (2006), Koch et al. (2005). To further accelerate temperaturedependent curing processes, many adhesives can be fastcured by conductive or inductive heating, e.g., Mahdi et al. (2003), Frauenhofer et al. (2012), Ratsch et al. (2020).
Based on the current investigation, the polyurethane adhesive system analysed can react in two ways: by temperature and by a diffusion-controlled reaction. Curing by temperature is recommended by the manufacturer for industrial applications. The combination of temperature and subsequent moisture curing is possible, as well as the other way around. At the same time, however, storage and curing under less than optimal conditions can have an unpredictable effect on the material if the diffusion of water into the adhesive is not taken into account. Therefore, it is necessary to investigate and model diffusive curing to make a reliable prediction of the material behaviour.

State of the art
The experimental investigation of reaction kinetics using differential scanning calorimetry (DSC) is common knowledge and practice. The approach using dynamic measurements (constant heating or cooling rates) has already been reported in Kamal and Sourour (1973) and Sourour and Kamal (1976), who developed the widely used the so-called Dr Michael Johlitz, Guest Editor of this Special Issue [New Challenges and Methods in Experimental Investigations and Modelling of Elastomers] confirms that where he was co-author of a research article in this Special Issue, he was not involved in either the peer review or the decision-making process of that particular article.
Kamal-Sourour reaction kinetics model. DSC experiments also allow the identification of the specific heat capacity; for instance, a comprehensive work on the identification and interpretation using DSC and temperature-modulated DSC (TMDSC) can be found in Lion and Yagimli (2008) and Lion and Yagimli (2009). For the identification of the diffusion coefficient for water in thin layers, a gravimetric measurement setup is utilised. Such gravimetric measurements of the vapour intake are primarily used in the food and pharmaceutical industry, e.g., Surface Measurement Systems Ltd. (2020). Density measurements using the buoyancy of the sample can be employed to investigate the volumetric shrinkage through curing, as is described in Kolmeder et al. (2011). In Leistner et al. (2018), the authors even define the evolution of the degree of curing based on the curing shrinkage. Furthermore, thermography measurements during the curing process have also been applied to validate the reaction kinetics model, see Leistner et al. (2019).
In this study, to incorporate the strain-rate-dependent elastic properties, a hypoelastic approach is considered. The most distinguishing feature of this approach is that the change in stress depends on a change in deformation, so that the progression of curing alone does not cause an increase in stress. For a more detailed description of the hypoelastic theory, refer to Hossain et al. (2009). Any desired hyperelastic model, see Ogden (1997) for an overview, could be considered based on this formulation. In Hossain and Steinmann (2014), the authors directly included the degree of curing in the their hypoelastic formulation in contrast to Hossain et al. (2009). Herein, a similar methodology is also employed. The investigation of the change in viscoelastic properties during the curing process using a shear rheometer was elaborated in Yagimli and Lion (2011) and Yagimli (2013). In Mahnken (2013), the influence of temperature and degree of curing on the viscoelasticity as well as the coupling between the two were treated.
The aforementioned Kamal-Sourour reaction kinetics model has been used in various applications, e.g., in the medical application of bone cement, refer to Kolmeder and Lion (2013). Landgraf (2016) gives a comprehensive overview of the literature concerning material modelling as well as an application of a Kamal-Sourour-type reaction kinetics model in the scope of thermally coupled finite viscoelasticity. In Landgraf and Ihlemann (2017) and Landgraf and Ihlemann (2018), the authors describe the locally defined curing of a PU adhesive that is modelled using a Kamal-Sourourtype model, as well. By introducing a so-called diffusion factor, the reaction kinetics model can also map the transition of a reaction from the temperature-controlled to the diffusion-controlled region, e.g., Rabearison et al. (2011). In view of these studies, a thermodynamically consistent model describing the temperature-controlled reaction between the resin and hardener together with the diffusion of water as a reactive agent constitutes a novelty. As an important remark, there is a parallelism between the modelling of curing and the thermo-oxidative ageing, e.g., Johlitz (2015) and Musil (2020), where the network reformation process is treated similar to the hypoelasticity and the reaction-diffusion system is considered for the diffusive ageing. Here, the mentioned similarities are exploited to construct a sophisticated continuum mechanical model. In Jennrich et al. (2020), the authors published a preliminary one-dimensional approach to the modelling for the same material used in this work.
In the framework of FE implementation, there are several points that may need to be addressed. First, the ordinary differential equations (ODEs) of Kamal-Sourour-type reaction kinetics models are inherently unstable and prone to influences of the time-step size. In Leistner et al. (2020), the authors introduced high-order, time-adaptive time integration schemes to counter this effect. Second, for the efficient solution of the inelastic strains in a finite strain viscoelastic Maxwell-type model, the robust single-step algorithm proposed in Shutov et al. (2013) can be applied.
The current study is organised as follows: the next section gives an overview of the literature regarding the experimental analysis and modelling of elastomers and adhesives. The following section elaborates on the chemistry of the polyurethane adhesive. The next section describes the experimental studies in detail. The following section introduces the Kamal-Sourour-type reaction kinetics model. The next section outlines the basic approach to the theory of mixtures. The following section presents the thermodynamically consistent mathematical model. In the next section, the staggered parameter identification process is outlined. The last section illustrates a representative finite-element example based on a well-known rheometer experiment.

General description of the investigated polyurethane adhesive
Sikaflex®-360 HC (Sika Deutschland GmbH, Stuttgart, Germany) is a stable one-component polyurethane adhesive that cures with heat to an elastomer in a few minutes for industrial applications, e.g., Cognard (2006). The adhesive can be used for dynamically highly stressed structural bonds and has been investigated, among others, for laser communication optics by LEICA Heerbrugg AG (1992) or applications in shipbuilding with regard to its curing behaviour by Wang et al. (2003). Typical curing conditions of this adhesive are approximately 10 min at 150 • C and longer than 20 min at temperatures above 110 • C . According to the technical data sheet, Sika Deutschland GmbH (2011), the adhesive should not be cured with air humidity alone as the final properties and curing time may vary. Due to its chemical structure, the adhesive is capable of building up a polymeric network with highly elastic properties through a number of cross-linking reactions. Even if the exact formulation is not known, some information on ingredients and possible resulting reactions can be acquired from the safety data sheet (SDS) Sika Deutschland GmbH (2016).
Prepolymers based on difunctional long-chain polyether polyols are often employed for such highly elastic, moisturereactive one-component adhesives and sealants. The monomer isocyanate used is composed almost exclusively of 4,4' methylene diphenyl diisocyanate (MDI), Fig. 1. The resulting highly flexible properties result from strictly linear isocyanate groups in the prepolymer. Isophorone diisocyanate (3-isocyanate methyl-3,5,5-trimethylcyclohexylisocyanate, IPDI, Fig. 2) is often used as cycloaliphatic diisocyanate in the production of moisture-curing polyurethane prepolymers with good photostability. The advantage of this isocyanate is the different reactivity of the primary aliphatic and the secondary cycloaliphatic isocyanate group. As a result, the products have a low viscosity, a narrow molecular weight distribution, and a low content of free isocyanate. Studies have shown that due to the different reactivity and cis/trans isomer contents in commercial IPDI products, a wide variety of reaction products can develop depending on the temperature and possible catalysts. A study of this is described in Lomölder et al. (1997). In principle, IPDI is less reactive than aromatic diisocyanates. It will therefore normally be catalysed with, e.g., dibutyltin dilaurate (DBTL), which is not mentioned in the SDS, Sika Deutschland GmbH (2016). However, at elevated temperatures, it also reacts without a catalyst. It decomposes slowly in water, producing nitrous gases and hydrogen cyanide. In the presence of water, the described prepolymers react by decomposing the unstable carbamic acid, forming polyureas, Fig. 3. The reaction proceeds relatively slowly and is diffusion-controlled. The SDS indicates the thermolabile (unstable at elevated temperatures) 4,4' methylenebis(aniline) (MDA, Fig. 4) complex with sodium chloride as thermal curing hardener. The relatively simple production of crystalline complexes of MDA with, e.g., alkaline earth chlorides was described by DuPont in VanGulick (1972). The crystalline structure of such complexes is well known since Jarvis and Owston (1971). An advantage of these complexes, which can be interpreted as blocked MDA, is that they only reversibly release the amine at a certain temperature and react with isocyanates to form polyureas (so-called thermolatency). For this reason, one-component, at room temperature, storage-stable adhesive formulations can be obtained. In the absence of water, the curing starts only when the required temperature is exceeded. In general, the possible reaction products in Sikaflex®-360 HC can result from the interaction of moisture and temperature curing as well as the fast reaction of the MDI-based prepolymer, the slower reaction of the IPDI prepolymer, and the latent hardener system.

Experimental investigations
The investigation of the reaction kinetics using DSC measurements is not sufficient to characterise the dual curing mechanism, since diffusion-based curing does not take place under the inert gas atmosphere of the DSC. Hence, the material is first pre-cured under controlled humidity and temperature for a defined time and then thermally post-cured in the DSC under an inert gas atmosphere. The  1 3 specific enthalpy of the thermal post-curing can be observed in the DSC, but the part of the diffusive reaction cannot be measured directly. By assuming that the specific reaction enthalpy is constant overall and is distributed between the thermal and diffusive parts, the diffusion-based part of the specific reaction enthalpy can be determined indirectly. Calorimetric investigations using DSC are usually carried out either at a constant temperature (isothermally) or at a constant temperature rate (dynamically). In the present case, only a dynamic measurement is appropriate as the reaction at elevated temperatures progresses rapidly. In the case of an isothermal measurement, the reaction would start as soon as the sample is heated up, but before reaching the isothermal measuring temperature, some of the specific enthalpy of the reaction would be missing in the measured data. In measurements with a constant temperature rate, after correcting the data by the baseline, one can extract the reaction peaks and determine the total specific enthalpy of the reaction. If the material has previously been pre-cured under a moist atmosphere, a smaller reaction peak can be observed. The signal measured by the DSC is the specific enthalpy rate ̇h . The area under the curve yields the specific reaction enthalpy. The reaction kinetics model is adapted to these two quantities, see Sects. "The reaction kinetics model" and "Parameters of the reaction kinetics model". The duration of pre-curing varies between 2, and 6 h, with 30% , 50% , and 70% relative humidity. The sample material is filled into aluminium DSC pans and smoothed at the top to ensure consistent diffusion boundary conditions. After pre-curing, the pans with the partially cured material are heated in the DSC from 20 • C up to 200 • C at a heating rate of 20 K min −1 , thus ensuring complete curing. The reaction mainly occurs between 80 and 160 • C . The fresh, uncured material is also dynamically tested with heating rates between 5 and 50 K min −1 , and the reaction peaks can be seen in Fig. 5.
The negative sign of the specific enthalpy rate signifies an exothermal reaction. The residual specific enthalpy can be observed to decrease with increasing pre-curing duration and pre-curing humidity. These are the parameters that drive the curing process, and the model must be able to reproduce this effect. Figure 6 shows the residual specific reaction enthalpy that can be measured in the DSC after pre-curing. The Fick's diffusion coefficient D is identified by sorption experiments in a gravimetric sorption analyser. The cured material is exposed to different relative humidities and temperatures, and the change in mass is observed until a state of saturation is reached. The diffusion coefficient at 25 • C and a relative humidity of 40% is chosen as a reference for the FE implementation, since this corresponds to real operating conditions.

The reaction kinetics model
As a measure for the progress of the curing reaction, the degree of curing q is formulated. It can be computed from the corrected, extracted measurement data as the integral of the specific enthalpy rate ̇h in relation to the total specific enthalpy H 0 : The degree of curing ranges between 0 and 1, where q = 0 describes the uncured state and q = 1 represents the completely cured state. Since in the present case, a dual curing mechanism is involved, two partial degrees of curing are introduced accordingly: the thermal degree of curing q and the diffusion-driven degree of curing q D . Since the total specific enthalpy is assumed to be constant in this modelling approach and is distributed between the two curing mechanisms, the sum of the partial degrees of curing also yields the total degree of curing. Based on the well-known reaction kinetics model of Kamal and Sourour (1973) and with reference to Yagimli (2013), we formulate the following ODE for the thermal curing: R is the universal gas constant. The empirical parameters E a1 and E a2 are usually referred to as the activation energy. is the temperature and R is the reference temperature, chosen as 20 • C . A 1 and A 2 are pre-exponential factors. The autocatalytic behaviour that can be observed in the experiment is represented by the term ⟨q + q D ⟩ m with the parameter m. Furthermore, the Macaulay brackets ⟨x⟩ = x+�x� 2 prevent the very small negative values that may occur at the beginning of the simulation. Similarly, the self-limiting behaviour is ensured in a later phase of the reaction by the term (1 − q − q D ) n with the parameter n in Eq. (2). This term further ensures that the sum of the partial degrees of curing can never become greater than 1. The term f step is a temperature-dependent step function that disables the development of q below the starting temperature of start = 248.55 K , see Eq. (5). The step function increases in value over a temperature range of = 3 K for numerical reasons. The PU adhesive investigated in Landgraf andIhlemann (2017, 2018) showed a sharply defined starting temperature for the curing reaction, as well. This behaviour could be observed in our experiments and is explained by the thermolabile chemical bond of the amine: For the diffusion-based curing model, the following ODE is formulated with reference to chemical ageing, e.g., Johlitz (2015) and Musil (2020): In Eq. (6), A D is another pre-exponential factor and E aD is another activation energy. In parallel with Eq. (2), it also contains a term for slowing down the evolution of q D when the reaction is near to completion. The physical interpretation here is that the structure of the cured material makes it difficult for reaction partners to interact. There is no autocatalytic term for diffusion-based curing in Eq. (6). The term c is the relative concentration: The superposition of different effects, for example, the chemical reaction of the water and the reduction in mass caused by the outgassing, makes it impossible to define the water concentration in the material during curing via the change in mass. It is therefore necessary to find an alternative definition that suits the experimental facilities. For identification purposes, the relative concentration is defined as the absolute humidity F abs relative to its expected maximum value under real conditions (at 40 • C and 100% relative humidity). The specific enthalpy of the diffusion-driven reaction H D cannot be measured directly. Instead, it is indirectly measured in post-curing DSC experiments under the assumption that the total specific enthalpy H 0 is constant and is split between diffusion-driven and thermal curing depending on the conditions. As long as a fully cured state is reached, the following applies: H D is the diffusion-driven specific enthalpy and H is the thermal specific enthalpy. The diffusion-driven degree of curing indicates which part of the total specific enthalpy is related to the diffusion-driven curing reaction: Keeping these assumptions in mind, the diffusion-driven specific enthalpy can be extracted from the DSC post-curing experiments: The index ' meas ' refers to measured quantities as opposed to ' sim ' referring to simulated quantities. The identification of the parameters is done according to the following relation, see also section "Parameter identification":

Theory of mixtures and the diffusion model
Based on the amount of fluid that can diffuse into the solid material under realistic circumstances, it is acceptable to treat the fluid molecules as a tracer inside the solid. The experiments have shown that at the maximum diffusion boundary conditions of 40 • C and 100% relative humidity, a water sorption of 1.4% by mass can be expected. Thus, it can be considered as the so-called problem type I as described in Hutter and Jöhnk (2004). For the treatment of different types of mixture problems, the reader is referred to Hutter and Jöhnk (2004) and Greve (2003). The model is formulated in the context of large deformations. Therefore, it is developed with respect to the actual configuration. The index ' f ' will mark variables and equations concerning the fluid (water), and ' s ' denotes the ones concerning the solid. If no index is specified, the variables refer to the mixture. The fluid concentration c f can be defined as the ratio of the fluid density f to the density of the mixture : The fluid concentration c f in Eq. (12) will be used throughout the modelling and the simulation. The relative concentration c in Eq. (7) is necessary solely for the identification of the diffusion-driven curing parameters. Then, the balance equation of mass for the mixture in the current configuration is given as: with the density and the velocity vector v . As the fluid phase is very small compared to the rest of the material (i.e., = s ), there is no mass balance of the fluid phase. As a consequence, the mixture mass balance equation [Eq. (13)] is considered in terms of the concentration: with the diffusion flux j f and the reaction term ĉ f , which will be explained further in Eq. (53). The diffusion flux in Eq. (14) is defined as: where f is the fluid density, v f is the fluid velocity, and d f is the barycentric velocity of the fluid as a tracer. Equation (15) will be characterised in the section "Diffusion and flux" by a constitutive law as a consequence of the Coleman-Noll procedure, refer to Coleman and Noll (1963), using the form of the entropy flux in Eq. (34), as formulated by Johlitz and Lion (2013).

Kinematics
The deformation gradient is multiplicatively split into two parts: an isochoric ̂ = J −1∕3 and a volumetric ̄ = J 1∕3 part, such that: The isochoric deformation gradient has the following property: Furthermore, the isochoric deformation gradient ̂ is decomposed into elastic and inelastic parts to incorporate viscoelasticity: e is the isochoric elastic deformation gradient and ̂ i is the isochoric inelastic deformation gradient. In consequence, the right Cauchy-Green deformation tensor is: Then, the isochoric part of the right Cauchy-Green deformation tensor, based on Eq. (18), reads as: Similarly, the elastic and inelastic right Cauchy-Green deformation tensors are defined as: In the present study, isotropic tensor formulations are considered. Therefore, for ̂ , the following relations hold: where ̂ is the isochoric spatial velocity gradient. The derivation of the specific hybrid free energy will be shown in section "Hybrid free energies". Finally, the time rate of ̂ e is given by the following relation: with the isochoric inelastic spatial velocity gradient:

Liu-Müller entropy principle
For models based on the theory of mixture, a different approach has to be considered to ensure thermodynamical consistency as explained in Liu (2002), Hutter and Jöhnk (2004), and Greve (2003). There are two main reasons: first, the existence of the absolute temperature is not valid in the mixtures; second, considering the entropy flux and the entropy supply s as arbitrary thermomechanical quantities is presumptuous. However, instead of deriving the equations anew using the Liu-Müller approach, the entropy flux defined in Johlitz and Lion (2013) is employed. Thus, the standard Coleman-Noll procedure can be applied. To begin this tedious and meticulous process, the independent variables P are chosen as: where ̂ , c f , and are the isochoric right Cauchy-Green deformation tensor, the concentration, and the temperature, respectively. Furthermore, the process variables Π , which have to be defined by constitutive laws, are given as: which are the second Piola-Kirchhoff ( 2 nd PK) stress tensor , the diffusion flux j f , the specific entropy s, the entropy flux , the entropy supply s , and the heat flux q . Following the approach of Johlitz and Lion (2013), by the principle of equipresence, see also Haupt (2002), all process variables Π depend on the same set of variables including the hydrostatic pressure p: The specific Helmholtz free energy is related to the specific hybrid free energy by a Legendre transformation as shown in Lion et al. (2014): where the volumetric strain is given as vol = J − 1 and R is the density in the reference configuration. It should be further noted that there are alternative definitions for vol , for instance in Simo and Taylor (1982). Similarly, the specific internal energy e is related to the specific hybrid free energy by: The rate of temperature ̇ is arbitrary; therefore, the specific entropy must be defined as: Since the pressure rate ṗ is also arbitrary, the volumetric strain vol is related to the pressure by: In accordance with Lion et al. (2014), we define the total stress power per unit volume in the actual configuration and then transform it to the reference configuration as necessary for the formulation of hypoelasticity (see section "Equilibrium and non-equilibrium parts") using the pull-back operation on the Cauchy stress tensor : The decomposition of the total stress power per unit volume into a volumetric and an isochoric part then reads, refer to Lion et al. (2014): with the isochoric 2 nd PK stress tensor ̂ . As explained in Johlitz and Lion (2013), the following ansatz for the entropy flux is chosen: Next, a general form of the Clausius-Duhem inequality suitable for the constitutive modelling in the context of this study will be introduced. First, Eq. (34) is inserted into the balance of entropy ( ≥ 0 ), where ̂ is the specific entropy production: Second, the time derivative of Eq. (28) is inserted intothe time derivative of Eq. (29). After the organisation of the terms and multiplying both sides with , the following relation is obtained: Next, with Eq. (36), the balance of heat has the form: Finally, modified balance of heat in Eq. (37) is inserted into the balance of entropy in Eq. (35), and after some mathematical manipulations, the Clausius-Duhem inequality reads:

Hybrid free energies
The hybrid free energy concept was first introduced by Lion et al. (2014) to be consistent with differential scanning calorimetry (DSC) measurements, where the isobaric specific heat c p is measured rather than the isochoric specific heat c v . In DSC, the hydrostatic pressure p is the main controlling variable. Based on this concept, the specific Helmholtz free energy is converted to a Gibbs-type specific free energy by making use of the Legendre transformation (see Eq. (28)). In particular, the Gibbs-type specific free energy is a function of the hydrostatic pressure, the isochoric right Cauchy-Green deformation tensor, and other variables. In essence, the use of the specific hybrid free energy only affects the volumetric part, and it is more accurate, since the measured isobaric specific heat is directly used in the model. For this reason, the hybrid free energy approach will be used in this work. For different interesting aspects about the hybrid free energy formulation, the reader is referred to Lejeunes and Eyheramendy (2020), such as its relation to mixed u-p models and thermomechanical splitting.
To embark upon the constitutive modelling, the decomposition of into a concentration part c , a volumetric part vol , and an isochoric part, which itself consists of an equilibrium part iso eq and a non-equilibrium part iso neq , is postulated: The specific hybrid free energy depends on the complete set of variables: but these dependencies can be distributed, as motivated by physical observations: Considering Eq. (41), the total time derivative of the specific hybrid free energy yields: with ̇i so eq and ̇i so neq as placeholders for the hypoelastic formulation of the equilibrium stress and the evolution of the which results in the derivative: a 1 is a parameter to be identified. The volumetric part of the specific hybrid free energy is defined as: c c p0 and c u p0 are the specific heat capacity at constant pressure for the cured and uncured material, respectively. c c p and u c p define the temperature dependence of the heat capacities. K is a constant, chosen as 1 K . c vol and u vol are the volumetric coefficients of thermal expansion for the cured and uncured materials, respectively. is the volume change caused by the curing process and is the bulk modulus. Taking Eqs. (31) and (41) into account, the partial derivative of vol with respect to the pressure p yields the volumetric strain: The second partial derivative of vol with respect to yields the equation for the specific heat capacity: where the linear dependence of the heat capacity on both temperature and degree of curing is obvious. The product of Eq. (47) with the temperature rate ̇ corresponds to the specific enthalpy rate ̇h , which is the outputquantity of a DSC measurement: The partial derivatives of vol with respect to q and q D are identical: Now, the only terms left undefined in Eq. (42) are ̇i so eq and ̇i so neq . They will be defined in Eqs. (55) and (64).

Diffusion and flux
For the diffusion flux, similar to Johlitz and Lion (2013), we define the following ansatz: where c is a parameter directly related to Fick's diffusion coefficient D by: Based on this assumption, the heat flux q can be determined from the following relation: where is the coefficient of heat conduction. From the mass balance of the fluid, we gain the reaction-diffusion Eq. (14), which shows that the change in concentration is composed of the change through diffusion flux and the chemical reaction, denoted ĉ f . Knowing the diffusion flux j f and the heat flux q , the reaction term in Eq. (14) can be defined as: where k r is a reaction parameter.

Equilibrium and non-equilibrium parts
The hypoelastic formulation of the equilibrium part is motivated by the dependence of the stiffness parameters on the degree of curing. The volumetric-isochoric split dictates the stress-free evolution of curing in the isochoric part, meaning that the isochoric stress may not increase based on the evolution of curing alone. For a change in the isochoric stress, there must always be a change in the isochoric deformation. Furthermore, Hossain et al. (2009) describe the simultaneous deformation and curing of the material, such that existing cross-links are deformed, while new cross-links develop in accordance with the already deformed state of the material. Thus, new cross-links cannot contribute to the previously developed stress, but will only contribute to the stress after the deformation is changed again. Based on that, the time derivative of the equilibrium part of the isochoric specific hybrid free energy is formulated as: where ̂ is the isochoric Green-Lagrange strain tensor. The fourth-order material tensor ℂ is further defined as: with the strain energy density function w and the correction factor eq . The strain energy density function is defined from a series of terms commonly referred to as the polynomial series. For incompressible materials, it takes the form: with the first and second invariant of the right Cauchy-Green deformation tensor I and II . In this case, using the isochoric invariants, the strain energy density function reads: with the stiffness parameters c 10 , c 01 , and c 20 . Its partial derivative with respect to ̂ follows as: with being the second-order identity tensor. Its second partial derivative with respect to ̂ reads: with the fourth-order tensors: The correction factor eq is necessary to guarantee that the isochoric stress tensor is indeed a deviator. It is defined as: For the time derivative of the non-equilibrium part of the isochoric specific hybrid free energy, we postulate: containing a stress contribution and the evolution equation of the Maxwell branches. Considering Eqs. (55) and (64), the Clausius-Duhem inequality now yields as follows: The isochoric 2 nd PK stress tensor also consists of an equilibrium and a non-equilibrium part: The correction factor for the non-equilibrium stress is defined as: With ̂ being the isochoric part of , the 2 nd PK stress tensor reads as: where ̂ is related to the deviatoric Cauchy stress tensor D through the following relation: The pull-back operation for is shown in Eq. (32). From the equilibrium part of the isochoric specific hybrid free energy, we obtain the time derivative of the isochoric second Piola-Kirchhoff stress tensor's equilibrium part:

3
A simple Neo-Hookean approach is used for each Maxwell spring with the following strain energy density function: where is a stiffness parameter. Inserting Eq. (71) into the first term of Eq. (64) yields the non-equilibrium stress for a number of n Maxwell elements: with the stiffnesses of the Maxwell springs k . Finally, the only remaining term in Eq. (65) will be used to determine the evolution equation for the viscoelasticity: Because of the isotropy, the following relation holds: Therefore, Eq. (73) takes the following form: which motivates the ansatz for ̂ i as in Johlitz and Lion (2013): symbolises the dashpot viscosity of a Maxwell branch and is a correction factor necessary to ensure that the evolution equation satisfies the deviatoric property ̂ i ∶ = 0 . It is defined as:

Considering the relation:
we can obtain the evolution equation for each Maxwell branch from Eqs. (76) and (77): where the relaxation time is defined as the ratio of the viscosity to the stiffness for every Maxwell element: To account for changes in the viscoelastic properties caused by curing and temperature changes, the basic relaxation times 0 are modified by a shift factor a sh in the following way: The shift factor is defined as: where the parameters a and a q are related to temperature and curing, respectively.

Parameter identification
The following paragraphs give a brief overview of the experiments used to identify the parameters, as summarised in Tables 2 and 3. Therefore, the advantages of the model formulation in terms of independent parameter identification of the different parts should become apparent.

Parameters of the reaction kinetics model
The parameters of q D are determined by comparing the diffusion-driven specific reaction enthalpy H D with the product of the diffusion-driven degree of curing and the total specific enthalpy, as shown in Eq. (9). The diffusion-driven specific reaction enthalpy H D results from the difference between the total specific enthalpy and the area under the reaction peak for thermal post-curing, see Eq. (10). The identification of the diffusion-driven curing parameters can be done independently of q . Figure 7 illustrates the comparison of the diffusion-driven specific reaction enthalpy derived from the measurements and the simulation. This can be considered as a good fit. The quality of the identification results varies (81) = 0 10 a sh .
(82) a sh = a ( − 0 ) + a q (q + q D − 1), depending on the chosen measurements. The selection of data points for the identification of the diffusion-driven curing includes a variation of both temperature and concentration ( c = 0.178...0.5 ). To provide a recommendation on the type and scope of the measurements (number of specimens, temperature, humidity, and pre-curing time), more investigations would first have to be carried out. For the parameter identification of q , the simulated specific enthalpy rate ̇h sim is generated from the time change of the degree of curing by multiplication with the total specific enthalpy, as shown in Eq. (11). Then, ̇h sim is compared to the DSC signal. The previously identified parameters of q D are used for the solution. This is done for the DSC measurements of fresh, uncured samples and for the residual curing of pre-cured samples in the DSC to determine the parameters for q . A comparison between the DSC measurement and the simulation for previously uncured material is shown in Fig. 8. The relative concentration c resulting from the relative humidity and temperature used in the measurements is summarised in Table 1. To have the relative concentration c better reflect the fluid concentration c f in terms of range, it could be normalised to a maximum value of 1.4% by mass fluid concentration, a value that was measured in equilibrium for the cured material. Such considerations are of a purely practical nature and have absolutely no bearing on the theoretical definition of the fluid concentration c f . This has not been done in the present work; the identified parameters were used 'as is' to showcase the model in the section "Modelling results".

Parameters of the volumetric part
In addition to the specific enthalpy of the curing reaction, the heat capacity as a function of temperature can also be measured using DSC. The measurements were carried out for both the uncured and the fully cured material in the same temperature range. Thus, the specific heat capacity parameters c u p0 , c c p0 , u cp , and c cp are identified as displayed in Fig. 9. The density of the uncured material R and the volume shrinkage through curing of approximately 5% are given by the manufacturer in the technical data sheet, Sika    (2011). They can be validated by density measurements using the Archimedes principle. The linear coefficient of thermal expansion for the cured material c lin can be identified using thermomechanical analysis (TMA) and subsequently converted to the volumetric coefficient of thermal expansion for isotropic material c vol . The corresponding identification results are shown in Fig. 10. At the intersection of the two simulated linear curves lies the glass transition temperature T g . For the volumetric coefficient of thermal expansion of the uncured material u vol , density measurements are applied.

Parameters of heat flux and diffusion
The coefficient of heat conduction for the cured material has been identified using the HotDisk transient plane source (TPS) procedure. The evaluation of the temperature drift indicates that the measurement is valid. For the diffusion flux, Fick's diffusion coefficient D is identified using gravimetric measurements under a constant temperature of 25 • C and a constant relative humidity of 40% . The parameter a 1 , cf. Eq. (51), is a constant of proportionality, making D compatible with the coefficient of heat conduction in terms of its unit.

Parameters of the equilibrium part
In Mezger (2019), the author summarises different approaches to the experimental analysis of equilibrium elasticity, and compares them in terms of the time expenditure and the quality of the results. One of the experiments with high quality of results is the multi-step relaxation test. It is used to exclude rate-dependent behaviour while identifying the parameters of the equilibrium part c 10 , c 20 , and c 01 . The maximum desired strain is approached in several equal steps. At each step, a certain amount of time is allowed for relaxation to occur, though the equilibrium state is never reached. From the maximum strain, unloading is performed in the same manner, so that there are two stress responses for each strain state, one on the loading curve and one on the unloading curve. Between these two points, the mean values are calculated which together form the so-called equilibrium curve. These measurements were performed on the cured material both in uniaxial tension and in uniaxial compression with very good results, as shown in Figs. 11 and 12. The tensile specimen geometry was S2 according to DIN Deutsches Institut für Normung e.V. (2009), and the compression specimen was a cylinder 10 mm in height and 20 mm in diameter.

Parameters of the non-equilibrium part
As a rule of thumb, the number of Maxwell elements chosen for a viscoelastic model should be two per decade. Nevertheless, seven elements have been found to be sufficient for this material and use case. The viscoelastic and temperature shift parameters are identified using cyclic shear rheometer measurements with a plate-plate geometry, constant frequency, constant small shear amplitude, constant gap, and constant temperature rate. The change in stiffness can be observed over the whole curing process and for a wide range of temperatures. Figure 13 shows the measurement and simulation for every fifth data point at the maximum shear strain.

Weak forms
The model introduced in the section "Mathematical model" is implemented in the commercial finite-element  Fig. 13 Identification of the viscoelastic parameters

20
• C framework COMSOL Multiphysics, which is advantageous when applying sophisticated continuum mechanical models like the one presented in this contribution. The partial differential equations are transformed into the weak form and applied to the software. In this way, the error-prone computation of the numerically consistent tangent matrices can be avoided. Another significant aspect is that during model development, new modelling ideas can be tested quickly and efficiently. Hence, the model can be changed and improved. The equations are solved with a second-order backward difference time integration method (BDF2), which is unconditionally stable, using linear elements. To solve the coupled system of partial differential equations (PDEs), a staggered approach is used. For the solution of the hydrostatic pressure, the mixed u-p approach is necessary to resolve the numerical problems associated with incompressibility, see also Chang et al. (1991). In this sense, a perturbed Lagrangian formulation with linear displacements and constant pressure elements is utilised. To begin with, the weak form of the heat flux Eq. (52) is calculated as: The term (⋅) denotes the test function in Eqs. (83), (84), (85), and (86). Then, for the fluid diffusion, the weak form reads as: Similarly, the weak form for the balance of linear momentum in the reference configuration can be computed as: which is inserted into the 'solid mechanics' physics environment in COMSOL. Also, the hydrostatic pressure p can be computed via the following weak form in a Lagrangian multiplier sense: Equation (86) is a direct result of the specific hybrid free energy method, which simply provides a mixed u-p form without requiring a variational formulation.

Modelling results
Because of the complex nature of the mathematical and FE model containing multiple couplings between temperature, concentration, reaction, and mechanics, smaller benchmarks are designed to showcase the different features. Another notable point is that the diffusion and thermal curing have a large difference in time-scale; therefore, it is necessary to investigate them separately.

Benchmark 1: diffusion
To observe the effects of diffusion, the first benchmark includes a very thin piece of material with a Dirichlet boundary condition for the concentration on the upper surface and zero flux boundary conditions on the other surfaces as well as a constant initial concentration in the bulk. The geometry, boundary conditions, and loading are shown in Fig. 15. This emulates the real conditions in a thin film diffusion experiment in a simplified manner. The temperature is constant and there is no mechanical loading. Figure 16 shows the diffusion of fluid into the material. The evolution of the concentration over time that can be seen in Fig. 16 looks physically plausible.

Benchmark 2: diffusion-reaction
The second benchmark visualises the progression of diffusion-driven curing and the decrease in concentration that results from it. Here, the diffusion flux is zero over all surfaces. Instead, there is a relatively high initial concentration in the bulk as the distribution of q D is homogeneous in this case. Figure 17 depicts the geometry and boundary conditions. The temperature is constant. Figure 18 shows the evolution of the diffusion-driven curing q D at an arbitrary point in the bulk, which leads to a decrease in concentration c f via the reaction parameter ĉ f .

Benchmark 3: thermal reaction and heat conduction
The third benchmark emulates the process in a rheometer with a plate-plate geometry but without mechanical (86) loading. The fluid concentration is 0 to show pure heat conduction. The temperature increases with a constant rate starting at room temperature up until 200 • C . Heat transfers on the upper and lower surfaces into the specimen. The geometry and boundary conditions are shown in Fig. 19. The partial degrees of curing q and q D start at 0 and q D remains 0 due to the concentration being 0. The progression of thermal curing becomes visible. Figure 20 shows the temperature and thermal degree of curing in the core of the specimen. Note that the evolution only starts after a certain temperature threshold is reached, which is due to Eq. (2). Moreover, to counter the inherent instability of the curing model's ODEs, Macaulay brackets have been introduced. A study on the influence of the time-step has shown no instability and no influence in the case at hand. The comparison for thermal curing based on benchmark 3 is displayed in Fig. 14. The evaluated history point is in the centre of the specimen.

Benchmark 4: curing under cyclic loading in simple shear
The final benchmark includes all the effects seen in the curing rheometer experiment. The geometry includes the aluminium plates, as displayed in Fig. 21. The temperature is increased from room temperature to 200 • C using a constant heating rate on the aluminium plates. On the outer surface, a zero flux boundary condition mimics the isolation through There is no reaction the surrounding air. There is no diffusion of water. The cyclic loading with a constant frequency and small constant amplitude is introduced at the upper surface, i.e., the loading curve with a frequency of f = 0.1 Hz depicted in Fig. 21.
The lower surface of the lower plate is completely fixed. Hence, the time-step is taken as t = 2.5 s . Displacement in y-direction is free for the upper plate. Figure 23 shows the evolution of the thermal curing ( q D stays 0) as well as the temperature in the centre of the specimen. The volumetric strain caused by both thermal expansion and curing shrinkage is shown in Fig. 24. In the beginning, there is only thermal expansion and no shrinkage ( q is still 0), and hence, the gap increases. Then, the thermal expansion and the curing shrinkage overlap, causing a decrease in the gap when the shrinkage dominates between t ≈ 400 s and t ≈ 800 s , and an increase when thermal expansion dominates above t ≈ 800 s . The evolution of q causes an increase in the stiffness of the elastic part. The viscoelastic properties change both through curing and temperature. Finally, the symmetric torque response of the upper face of the specimen is demonstrated in Fig. 25. The torque response is calculated by the wellk n o w n r e l a t i o n , s e e S o k o l n i k o f f ( 1 9 5 6 ) , T R = ∫ S R P 23 X 1 − P 12 X 3 dA = ∫ S R T A dA , where T A can be denoted as torque per area. The terms P 23 and P 12 are the components of the 1 st PK stress tensor. Furthermore, Fig. 22 displays a symmetric distribution of torque per unit area T A for two different times, namely at the beginning of the simulation t = 107.5 s and near the end of the simulation t = 1107.5 s.

Summary and outlook
In the present work, a comprehensive experimental analysis and a multiphasic modelling approach have been introduced for the thermal and the diffusion-based curing of a polyurethane adhesive. The chemical components of polyurethane adhesives, and their expected effects on the macroscale and on the finished product were summarised. For the modelling of the material behaviour, a physically motivated thermodynamically consistent continuum mechanical model was introduced. The resulting model is applicable in the framework of an FE simulation. In consequence, the weak form of the mathematical model was implemented in the commercial FE software COMSOL Multiphysics. Another interesting point is that the modular nature of the model makes the identification of material parameters straightforward. The experiments necessary for the identification of the reaction kinetics, which include dynamic DSC and gravimetric measurements, were described in detail as well as the corresponding identification results. The application of a hypoelastic formulation for the equilibrium part enabled the stress-free evolution of curing considering its profound influence on the material stiffness. The viscoelastic formulation included the influence of temperature and curing on the relaxation behaviour through a shift of relaxation times in the Maxwell model. Several smaller benchmarks were used to demonstrate the implemented complex FE model successfully capturing the effects of diffusion, heat conduction, curing, and mechanical loading. For future studies, the complex model will be validated in comparison with a component test. Furthermore, a parameter study on the curing in a high-temperature rate environment created by inductive heating will give an indication of the