Thermo-mechano-chemical modeling and computation of thermosetting polymers used in post-installed fastening systems in concrete structures

As thermoset polymers find frequent implementation in engineering design, their application in structural engineering is rather limited. One key reason relies on the ongoing curing process in typical applications such as post-installed adhesive anchors, joints by structural elements or surface-mounted laminates glued by adhesive polymers. Mechanochemistry including curing and aging under thermal as well as mechanical loading causes a multiphysics problem to be discussed. For restricting the variety of material models based on empirical observations, we aim at a thermodynamically sound strategy for modeling thermosets. By providing a careful analysis and clearly identifying the assumptions and simplifications, we present the general framework for modeling and computational implementation of thermo-mechano-chemical processes by using open-source codes.

design-life of the structure. Different phenomena in many length scales [7][8][9][10][11] necessitate incorporation of multiscale approaches. Chemical reaction leading to hardening (hydration) needs to be considered by including dissipative effects [12] also caused by the internal friction during deformation [13,14]. Mechanical response change related to the loading rate [15,16] or modeling damage in concrete [17][18][19] has to be modeled and solved together with the hydration in order to reach accurate prediction of the design-life. Differing bonding mechanisms between structural parts like concrete and steel is an active research area. In the case of a mechanical interlocking, this so-called apparent adhesion [20] based on friction contact provides only a limited bond with respect to a chemical or dispersive adhesion. Amending adhesion is possible by gluing concrete and steel [21][22][23] that increases also the time to failure [24]. Such a chemical bond is achieved by a special type of polymer, specifically in this work, we concentrate on this adhesive material at the interface [25] between concrete and steel anchor.
Effected by their high stiffness properties, thermosetting polymers are proposed to be used at the interface. They are already available in the market; mostly an epoxy or vinyl ester is used. By mixing two fluid-like components, one is a multifunctional comonomer and the other one is the hardener, the chemical reaction called curing is ignited forming cross-links (chemical bonds) between polymers leading to an increased stiffness. The curing process never stops until it reaches the fully cured state. This process is of paramount importance for correctly characterizing thermosets. Until the so-called gel point, the process is fast and causes a significant shrinkage leading to residual stresses [26], after the gel point we consider the material as a solid-to be precise a fluid with an infinite viscosity-leading to increased stiffness with further cross-linking. It is possible to track the curing process experimentally [27] by introducing a degree of curing, experimental determination of mechanical response [28][29][30] is an important parameter to calculate the long term response [31]. Mechanical characterization [32] and thermomechanical considerations [33] especially below the gel point [34,35], or considering aging [36,37] alongside to the curing have been studied. The effect of the degree of cure on the mechanical response is inevitable. Various models have been used by exploiting numerical multiscale approaches [38][39][40] or also simple linear models. Multiphysics approaches with adequate numerical implementations are developed as well [41], not only for thermosets using small displacement assumption [42,43], but also for large displacements [44][45][46], has been applied in several applications [47][48][49], where the curing is established in a controlled environment (i.e. heat treatment of composite elements in the autoclave) such that the material characteristics yield the chemically best configuration. In other words, the polymer attains a fully cured state and the curing process stops. However, this is rarely possible in case of thermosetting polymers that find application in the construction industry and especially reinforced concrete structures. Owing to the nature of constructing a structure in an uncontrolled environment under definitely not optimal conditions, the thermoset polymer will typically fail to reach a fully cured state after the installation such that the curing phenomenon continues during the life-time of the structure. Considering the composition altering due to the curing and its effect on the design-life of the system is of paramount importance. In order to model the whole system response accurately, in this work, we start with a thermodynamically consistent modeling of mechanochemistry in thermosetting polymers. We present a computational implementation with the finite element method by using the open-source packages developed under the FEniCS project [50,51].

Modeling the curing process
The underlying thermosetting polymer is a resin being hardened due to the chemical reaction, called curing, leading to cross-links between existing polymer chains. Cross-linking affects (often decreases) the volume per mass, an increase in the overall stiffness, and at the same time altering the temperature due to an exothermal reaction. After mixing an agent to the resin, curing starts by using energy from the environment-mostly heat or radiation energy (UV light) is supplied to level off the temperature. Once the curing starts, from the viscous resin state until the so-called gel point, the curing process is fast. At the gel point the viscosity increases drastically so that we may call the material a solid and the curing process slows down; its rate asymptotically reaches zero at the fully cured state. Degree of cure is a difficult variable to measure. We start with a definition of this variable following [52] by using theory of mixtures in a simplified manner. In a unit volume, V , masses of resin, curing agent, and cured solid are introduced m R , m A , and m S , respectively. Their values vary in time; but the sum is constant, m = m R + m A + m S , throughout the curing process. We introduce mass fractions: leading to Now we introduce the degree of cure, ω = ω(t), as "measured" in time, giving the mass fraction of the solid under the assumption that the initial mass of solid is zero, Y S (t = 0) = 0, such that we have Consider that the ratio of masses, Y R (t = 0)/Y A (t = 0), remains the same during the reaction. In other words, same molar mass is used for resin and agent. Now by using z = Y R (t = 0), we obtain Hence, the whole formulation for calculating the masses of constituents has been subsumed to the degree of cure (conversion degree), ω. For obtaining its value, an evolution equation is developed dependent on the chemical reaction type. Epoxy resin reaction is based on autocatalysis mechanism [53], effected by hydroxy groups formed during catalysis. Hence, the following evolution equation [54] is found to be useful where amplitudes, A × , activation energies, E × , and powers, m, n, are constant parameters to be determined. Universal gas constant, R, is known. Temperature, T , is the key variable altering kinetic rates, k × , and thus curing rate, ω • , significantly. We refer to [55, Table 1] for parameter values of different materials. Obviously, ω = ω • dt is used for determining the current degree of cure in each material point. This model fails to incorporate vitrification such that different mechanisms below and above glass transition temperature, T g , are not modeled accurately. Especially curing at low temperatures (below 100 • C) leads to incomplete conversion that is of paramount interest in fastening systems, where the hardening occurs in environmental conditions. A possible approach as in [56] characterizes glassy and rubbery states occurring simultaneously. The idea was introduced in [57] based on phenomenological observation combined with the free volume reduction during cure leading to decrease in mobility. Mobility plays a dominant role in the rubbery state as molecules collide and form a network via cross-linking. Beyond a cross-linking density, material vitrifies and in this glassy state the chemical rate is more dominant. By following [58], we use only primary and secondary amines via autocatalysis and impurity catalysis (denoted by c in index), collective kinetic rates, K 1 , K 1c , involve diffusion controlled mechanism, K diff , as well as chemically steered mechanism, K 1,chem , K 1c,chem , as suggested in [59] as follows: The diffusion rate is modeled by the Williams-Landel-Ferry (WLF) equation based on [60,61] in the rubbery regime and greater than chemical rate in many orders in magnitude. In the glassy state, diffusion rate is nearly zero regarding chemical rate. This interplay delivers a realistic prediction of conversion degree, ω, at low temperatures, where a so-called partial freezing inhibits to attain ω = 1. This model incorporates glass transition temperature, which is often determined by a fit function [62] based on calorimetric measurements. For consistency, we follow [58] and use the relation from [63] with C = c ∞ / c 0 as the ratio of the heat capacity change at the glass transition temperature of the fully cured network with ω = 1, T g,∞ per the heat capacity change at the glass transition temperature of the monomer with ω = 0, T g,0 . The necessary coefficients for the evolution equations (6), (7) are determined by inverse analysis to data obtained by differential scanning calorimetry, we compile them from the literature as presented in Table 1. By knowing this mass fraction and the current mass density of the mixture, we can extract masses of each constituent by using Eq. (5) with z being the mass fraction of resin initially.

Thermodynamical formulation
Under the simplification that deformation is small such that no distinction is necessary between configurations, we use () • as the partial time derivative and ,i for the space derivative with respect to X i , where space coordinates, X, denote material particles composing a material system. Hence, the mass balance is satisfied and the mass density, ρ, is a given function in space-for homogeneous materials we set ρ = const in space, the presented formalism holds true for heterogeneous materials as well. We need to solve balances of momentum and internal energy, respectively, where u in m is the displacement to be calculated; stress, σ in Pa=N/m 2 , heat flux, q in W/m 2 , and specific internal energy (energy per mass), u in J/kg, need to be defined; g in N/kg is the gravitational specific force, and r in W/kg is the given specific supply term. For the sake of simplicity, we use a linear strain measure: However, we emphasize that the formulation is also suitable for nonlinear measures. Internal energy needs to be defined by exploiting thermodynamics. There are several methodologies; we refer to [64] for such examples based on thermodynamics of irreversible processes and non-equilibrium thermodynamics. Briefly, we introduce the so-called specific Helmholtz free energy: with the specific entropy, η, to be defined indicating the reversible part of heat flux and temperature, T , to be calculated. Free energy is modeled by assuming that it depends on temperature, T , strain, ε, and degree of cure, ω, We stress that we choose to exclude a dependency on the rate of these variables. Now by inserting the free energy into the balance of internal energy, dividing by T , and after rewriting, we obtain Based on the nonequilibrium thermodynamics, we decompose into equilibrium and nonequilibrium parts and then use the well-known relations for equilibrium parts For simplicity we assume that rate of temperature and rate of strain fail to alter mechanical response, leading to η neq = 0, σ neq = 0. Equation (15) 1 is often introduced as the definition of the entropy and Eq. (15) 2 is Cattaneo's theorem in its general form for small displacements, we refer to [65] for its derivation from a Lagrange function. We stress that these relations show that (equilibrium) entropy and (equilibrium) stress depend on the same arguments as the free energy, η = η T, ε, ω and σ = σ T, ε, ω , which is often called the principle of equipresence [66, §293.η]. By using these relations, we arrive at the balance of entropy: with the right-hand side called the entropy production that is zero for reversible, and it is positive for irreversible processes. By using this assertion-the second law of thermodynamics-we obtain Fourier's law and the restriction assuring that the entropy production is (zero or) positive as long as κ ≥ 0. The former is well known in thermodynamics. The latter is equivalent to the postulated Karush-Kuhn-Tucker relation originally introduced in plasticity, for a general discussion we refer to [67]. We emphasize that the chosen evolution equation for ω in Eq. (6) or (7) is positive such that ∂f /∂ω has to be negative. Heat is inserted into the system in a standard calorimetric measurement, δQ = c dT , with specific capacity, c. Now by assuming the relation T dη = δQ, as proposed in [68] and used in all thermodynamical formulations, we obtain All these relations will be fulfilled by defining the free energy adequately. After inserting Eq. (17) 1 into Eq. (13) 2 and using η = η T, ε, q as well as Eq. (15), we acquire This governing equation for the temperature is a general formulation under the assumption that the free energy depends "only" on temperature, strain, and degree of cure. By defining the free energy adequately, the governing equation will be closed and then calculated numerically.

Thermo-mechano-chemistry of thermosets
We aim at calculating the deformation by satisfying and the temperature by fulfilling for a hardening thermosetting polymer, whose state is characterized by the degree of cure, ω, modeled by the evolution equation (6) or (7). In order to close these governing equations, we need to model the system by defining a specific free energy, f . The whole formulation is reduced to define a correct model for the free energy.
In the following, we give a possible model and present its relation to the existing models in the literature. We emphasize that the free energy may depend on temperature, strain, and degree of cure. The real dependence is possibly determined by experiments; there are no thermodynamical restriction about its form. But we know, based on the numerical stability and also partly on our understanding of energy as well as intuition, the energy needs to be positive definite. In full accordance with Eqs. (18), (17) 2 , we propose to use the following free energy expression, where all material parameters, c, C, α, β, H ref , may depend on the same arguments as energy, namely temperature, strain, and degree of cure. This model is the simplest case for a thermoelastic material; the form is quadratic providing stability. The reference temperature, T ref , is often set as the room temperature under two assumptions: the material possesses no thermal stresses and the coefficient of thermal expansion (CTE), α, is measured regarding this temperature. Curing-related shrinkage is given by the material parameter, β. Heat release because of the exothermic reaction is given by H ref , the total heat (per mass) generated by a fully cured specimen. If the material response is such that c, C, α, β are constant in ε and T , in other words we neglect hyperelasticity as well as temperature-related changes in material coefficients, then we acquire the (generalized) Hooke's law with a linear shrinkage and well-known Duhamel--Neumann extension Analogously, the specific entropy and its derivative with respect to the strain read In the general case, we expect to have material parameters depending on the cross-linking. A rather obvious observation is made for stiffness, the material hardens as the cross-linking density increases. Accurate modeling of the relationship between material parameters and degree of cure needs to be achieved by experiments [69] by means of detecting ∂f /∂ω in the aforementioned general formulation. A possible way of introducing the degree of cure is simply to introduce two phases, solid (hardened) and fluid (non-hardened), creating the stiffness in the following simple rule of mixture: where C s , α s and C f , α f mean the maximum stiffness, maximum CTE reached after full cross-linking; and, initial stiffness, initial CTE of the non-hardened fluid alike material, respectively. Even the specific heat capacity, c, may depend on the degree of cure, ω, we direct to [70] for such a study. The simple relation between material parameters and degree of cure-motivated by the mixture theory-has a formal benefit of constructing thermodynamically sound material properties [71, Sect. 4]. Since the parameter for shrinkage, β, is inherently related to the curing process, the explicit (linear) dependence is established in the free energy definition instead of introducing phase depending coefficients. In this manner, we have such that we reach from Eq. (21) leading to Especially for modeling fastening systems in concrete, all assumptions are easily justified. Furthermore, for the sake of a better comprehension, we introduce more simplifications as follows: 1. In the material state before cross-linking, stiffness of the fluid-like or gel-type material vanishes, C f = 0, and partly as a consequence, a reversible thermal expansion is negligible, α f = 0. Then Eq. (27) for the temperature is obtained as follows: 2. A further simplification, especially beyond the gel point, relies on the fact that the shrinkage and thermal expansion are small. By neglecting these effects, we obtain Especially by neglecting deformation, ε = 0, we end up with the formulation used in the literature often for thermal analysis during hardening.
3. Further assumptions emerge in calorimetric measurements, without a supply term, r = 0. We write the global form of the latter, after using Gauss--Ostrogradskiy theorem and Fourier's law, The surface normal, n i , is outward the body such that that the second term is the heat released from the body. In an experiment, the heat supplied to the body is measured, In the case of an experimental setup, where the specimen is almost free on all sides, we may assume that no energy is stored. Furthermore, as the specimen size is small, we assume constant temperature within the body as well as constant rate of curing degree. Heat flow, Q, and temperature, T , are controlled and measured in the experiment for the specimen of mass, m = B ρ dV . For a known specific heat capacityeven if the heat capacity or H depends on the temperature, we stress that the temperature is constant within the body, B-we calculate The mass, m, and mass density, ρ, are constant, and we obtain H ref and parameters in the evolution equation for ω • . Two special cases arise in the measurement. The first case is the isothermal calorimetric measurement,

Generating the weak form for the finite element method
We consider a general case and set the objective to compute the primitive variables, ω, u, T , in space, x, and time, t. By starting with the initial conditions of a partly cured body given by ω ref as follows: with the evolution equation for ω • . We use the quadratic free energy as in Eq. (22), leading to linear material equations; yet emphasize that the implementation is designed in such a way that more advanced models are possible to be employed as well. Therefore, we leave the free energy as f in the following. All continuous functions are approximated by using their discrete representations in space and time. In order to discretize in time, we use constant time steps and compute in a time series with (·) 0 , (·) 00 denoting the calculated values of one, two time steps before, respectively. Time derivative is approximated by the finite difference method (Euler backwards scheme) The governing equations become We remark the linearized strain measure in order to find out the value one time step before, For space discretization, we use the finite element method with a triangulation of the continuum space, B, and exchanging continuous primitive functions with their counterparts in a Hilbertian Sobolev space [72] by using nth polynomial degree form functions for each discrete element, For the sake of notational easiness, we simply skip to distinguish between continuous functions and their discrete representations. We use the standard procedure to construct the integral forms by building residuals from governing equations and multiplying (weighting) them by corresponding test functions, δu, δT , where they are expressed in the same space as the primitive variables, called the Galerkin procedure, We emphasize that stress already possesses space derivative of primitive variables, because of the additional divergence on it, we need to choose at least n = 2 in the space discretization. This continuity condition is weakened by integrating by parts for all second derivative terms, where we have insertedt i = n j σ ji = n j ∂f /∂ε i j andq = n i q i = −n i κ T ,i on boundaries. For the displacement with a given traction,t in Pa, two sets of boundaries appear: Dirichlet boundaries, ∂B D , and Neumann boundaries, ∂B N . We restrict the formulation such that these two sets fail to intersect, ∂B D ∩ ∂B N = {}, so the whole boundary is the union of them, ∂B = ∂B D ∪ ∂B N . At the Dirichlet boundary, the displacement itself is given, hence, no test functions are necessary, δu = 0 ∀x ∈ ∂B D . For the rest, traction is given, t = given ∀x ∈ ∂B N . Analogously, the temperature may be modeled by given temperature or heat flux; however, we use a more natural approach by using a mixed boundary condition called Robin boundaries, herein applied to the whole boundary,q with the convection parameter, h, and the ambient (environment) temperature we use as T a = T ref . Two weak forms are in different units, the former is in the unit of energy and the latter is in the unit of power. After multiplying by the constant time step t, we can sum them up and use for computing thermomechanics of a curing thermoset.

Computational implementation and numerical examples
In order to examine the strength and accuracy in predicting the behavior of realistic systems, we implement the weak form in Eq. (44)  First, the weak form is nonlinear such that we use the conventional Newton--Raphson linearization approach for solving a linear problem and update the solution by starting from the last converged solution. Even if the system is governed by stiff differential equations, using small time steps leads to a convergence in each iteration. The linearization is handled by exploiting symbolic differentiation that allows to implement also highly nonlinear forms-for example a hyperelastic material model can be easily implemented with the approach developed herein, we refer to [64] as well as to the supply code of [65] for examples of such an implementation.
Second, the specific free energy in Eq. (22) has been implemented also with the aid of symbolic differentiation. Technically, the whole material formulation is reduced to the free energy definition such that the implementation is very general and can be easily adapted for other systems simply by changing the free energy. We stress that the material modeling is limited by the assumed free energy formulation; however, implementation is novel in the sense that any more advanced formulations are possible to utilize. We make the code publicly available in [73] for encouraging its use under GNU Public licensing [74] for an efficient and transparent scientific exchange.

Curing models comparison
Two evolution equations (6), (7) have been introduced and implemented. Both equations consist of parameters depending on the current temperature such that the weak form F T in Eq. (43) is nonlinear. Material parameters belong to typical epoxy type of materials with coefficients compiled in Tables 1, 3. We compute temperature and degree of cure for a simple geometry of 50 mm×10 mm×10 mm, made of homogeneous material. First, we solve the temperature distribution by using model A in Eq. (6), and secondly, we obtain the temperature by applying model B in Eq. (7). Three-dimensional simulation results are recorded, and both results at the same location (geometric midpoint of the domain) are compared qualitatively. We emphasize that the material parameters in Table 1 are realistic; but from different sources such that models A and B fail to be compared quantitatively.
We begin with the well-known effect of the ambient temperature on the curing phenomenon. By setting the convection parameter high, h = 1 kW/(m 2 K), we enforce nearly a constant surface temperature equivalent to the ambient temperature, T amb . Applying model A and model B leads to results in Fig. 1. As expected, higher temperatures increase the conversion rate such that the possible upper limit is attained in less than 90 minutes at 90 • C than in 20 hours at 50 • C. These simulations model an isotherm calorimetry measurement. Model A is capable of quantifying the process only in the case of complete curing. We point out the case of an incomplete curing taking place at low temperatures. The model B given by the evolution equation (7) is indeed capable of simulating this freezing behavior effected by the vitrification. Such a vitrification, 10-20 • C above the curing temperature, is typically observed in isotherm calorimetry measurements. This glass transition temperature, T g , is also detected and measured with the aid of this observation. During cross-linking, the exotherm reaction produces heat which is taken out of the system in order to hold the temperature constant. For simplicity, we use a constant value for the specific heat capacity in simulations. In reality, the change of specific heat capacity between glassy and rubbery states causes an overshoot in this heat energy, which is used to determine the numerical value of T g in a non-fully cured structure. Material showing a distinct T g justifies the use of model B. A typical example are the epoxy systems used in post-installed fastening applications subject to low (ambient) temperature curing.
For a better understanding of the role of the evolution equation on the released heat, we simulate a temperature ramp test, where structures with different degrees of cure are simulated as being in a temperature chamber where the ambient temperature increases by 10 K per minute. Figure 2 demonstrates the specific heat flow, Q/m, obtained by Eq. (32) from the temperature and degree of cure computation. Obviously, model A is symmetric leading to the sigmoid relation, whereas model B shows a kink at the glass transition temperature such that the symmetric character is lost.
In general, the computational implementation results in reliable simulations with realistic parameters. Such a tool is of importance to study and comprehend the behavior of different curing models. We have selected two models without assessment of the quantitative accuracy to a specific material.  . 3 For a comparative analysis, curing of polymer adhesive material as used in Brussels (Belgium) on October 1, 2019. Left: Using the temperature change over the day by using 5 picked values as given in Table 2. Right: Mean temperature value applied on the boundary and shown in the core of the sample

Curing under alternating thermal loading
For a realistic simulation, we use model B and consider that the epoxy with parameters in Table 1 has been used for adhering components in an outdoor facility. Complete curing of the material is desired for an adequate stiffness leading to high adhering. We obtain hourly temperature, as an example on October 1, 2019 in Brussels, 1 and use with a rough approximation as given in Table 2. By using a convection parameter h = 100 W/(m 2 K), we model a mild windy ambient for heat exchange over the surface. Again for a simple geometry of 50 mm×10 mm×10 mm, we compute transiently the temperature distribution in tree-dimensional continuum coupled to the degree of cure given by the model B. The simulation starts directly after casting, ω = 0, computes for 24 hours temperature and degree of cure. The temperature values in the core of the structure are recorded, which are nearly the same as the surface temperature since the structure is small. We emphasize that ambient temperature and structure's temperature are not identical. Convection over the surface as well as curing produced heat are responsible for this difference. Important observations are attained from results demonstrated in Fig. 3. After one day of curing in environmental temperatures, only ω = 0.35 has been reached in Fig. 3(left). With high mobility, rate of conversion is accelerating in the beginning, especially after the early morning when the temperature starts rising. As a coincidence, on this day the temperature drops before noon such that the material presumably vitrifies and afternoon temperature increase is only partly amending the conversion rate. If the material has a gel point of ω = 0.5, then 24 hours of curing in a mild autumn day reaches not even beyond the gel point. Obviously, the expected stiffness is not delivered in this condition, at least for this material. Therefore, different commercial products use modified compositions, allowing an increased reaction kinetics in these temperatures. By utilizing the implementation herein, we enable investigating the degree of cure of a material with known parameters. Moreover, we stress that the temperature history is of importance as well. When we use the mean value for the ambient temperature, instead of the alternating temperature, we observe that a significant difference of approximately 10% is captured in Fig. 3(right). Even a linear relation between stiffness and the degree of cure is assumed, this 10% deviation is of importance to consider in design. This analysis is for a small sample, where the temperature conductivity is high enough to generate a homogeneous temperature distribution. However, the exothermal reaction is still altering the temperature in the middle of this sample as visible in Fig. 3(right). But this change is insignificant. Herein, we demonstrate how valuable such a study becomes especially whenever applied to more realistic geometries.

Mechanical loading superposed by hardening
Under a mechanical loading, the amount of deformation depends on the stiffness that changes during hardening. Two possible implementations can be developed. One approach is to use a stiffness, C s , for the fluid phase leading to a soft matter even before the gel point. Then we start with ω = 0 and with the aid of this initial stiffness, the numerical implementation converges. Another approach is to set the fluid stiffness to zero, C s = 0, α s = 0, and start with an initial curing, ω ref , greater than zero such that the simulation begins with an initial stiffness again leading to a convergence. We choose the second approach, for example the gel point, ω ref = 0.5, where the material response is more solid type than fluid. We use the latter approach in the following simulations. Materials response is chosen to be elastic and the initial geometry is used as the reference configuration to which the continuum body recovers after unloading. This modeling is based on the fact that the hardening beyond ω ref fails to affect the reference configuration. We leave an experimental verification of this assumption to further research.
A simple, thick beam of 50 mm×10 mm×10 mm is modeled. Material is assumed to be isotropic such that the solid properties are given by Lame's constants, λ, μ, coefficient of thermal expansion, α, and shrinkage parameter, β, The Lame parameters are related to the (measurable) engineering constants, E, ν, as follows: We use fictitious but realistic material coefficients as compiled in Table 3. Nearly all coefficients can be determined by experiments on the corresponding material. But the thermal convection parameter modeling the rate of heat exchange between the surface of the material and environment depends on the material, surface roughness, and ambient state [75]. For example, moving air is known to extract more heat from the material than still air. Analogously, humid air possesses a higher heat transfer coefficient than dry air. We consider a uniaxial tensile test, the beam is clamped on one end and pulled on the other end with a controlled force, first, relatively quickly until it reaches its maximum and then, secondly, it is held at this force throughout the whole simulation of 10 minutes. In the case of an elastic material, the response is instantaneous, the beam gets stretched and remains at the same length under constant force. This case occurs in Fig. 4(top), where the ambient temperature is less than T g such that the material remains at the same degree of cure, ω ref . In the case of an ambient temperature higher than the glass transition temperature, with the additional assumption that Young's modulus remains the same, response to the applied force is the same; however, the material undergoes further curing. Hence, the material hardens, stiffness increases, and the beam gets shorter under the same force held constantly, see Fig. 4 (bottom). Since we have modeled stiffness to degree of cure linearly, the length change is linear as well. In reality, Young's modulus depends on the temperature, especially around the glass transition. The resulting decrease of the modulus is significant. Herein we neglect the viscous character and model this behavior by using inverse tangents function and constant moduli for rubbery and glassy states, E r and E g , respectively, as follows: For the simulation, we use again fictitious but realistic parameters, E g = 2 GPa and E r = 0.5 GPa, leading to a stark decrease at the T g , for example as in Fig. 4. We emphasize that the glass transition temperature, T g , increases with evolving degree of cure given in Eq. (8). Therefore, we repeat the last analysis and apply a given boundary condition, hold this displacement, increase the ambient temperature, measure the necessary force. In the case of constant Young's modulus, as we model (linear) elastic material, the response remains the same as long as no postcuring occurs. Since Young's modulus decreases at the time crossing the current glass transition temperature-we stress that the simulation starts off with ω ref -stiffness drops because of decreasing Young's modulus. Furthermore, curing begins and material hardens as well. Such a response is  experimentally observed, we refer to [76] for a discussion in this sense. Herein, we observe analogous response in Fig. 6. The curing rate depends on the local temperature that converges to the ambient temperature steered by the convection parameter, h. Even for the same ambient temperature T ref = 300 K, a variation of the convection parameter shows that the increased heat exchange, larger h in the simulation, increases the time necessary for attaining the same degree of cure. This phenomenon is demonstrated in Table 4 quantitatively, where the time is determined, at which, in the middle of the same structure, ω = 0.95 has been reached. The justification of this phenomenon is as follows: Increased heat exchange leads to the increasing energy loss to the environment in form of heat energy that is generated as a consequence of the exothermal curing reaction. Obviously, smaller heat transport parameter represents a better isolation from the environment, practically sealing the structure. The generated heat during cross-linking is captured within the material, as seen in the temperature evolution at the middle point of the beam demonstrated in Fig. 7. In other words, the heat energy, which is released free at the moment of curing, increases locally the temperature and accelerates the curing process. The evolution equation for curing depends on the temperature.

Conclusion
A general framework has been developed for thermomechanics of thermosetting polymers under curing governed by an evolution equation. The numerical implementation of such a thermo-mechano-chemical process has been established by using open-source packages. Simple geometries are utilized for promoting a better understanding of underlying multiphysics as well as capability of the computational implementation. The following assumptions have been undertaken: • Small strain formulation is used; there are several applications especially in composite materials, where this assumption is accurate. • Dissipative parts in the thermal and mechanical terms are neglected. Depending on the chosen material, a viscoelastic term might be necessary. • The evolution equation for curing depends solely on the temperature. The temperature dependency is intuitive; however, it is not known, if a mechanical deformation would increase the rate of degree of cure. • The reference configuration is chosen at the gel point, which has been the initial condition as well. Although this assumption is very tempting, we leave an experimental validation of this assumption to further research.
Thermodynamically sound modeling of thermo-mechano-chemical processes allows us to accurately model the coupling effects between variables like temperature, deformation, and degree of cure. An important contribution is that we have subsumed all modeling by using the free energy such that any possible restrictions of the Fig. 7 Local temperature evolution for different cases by varying heat transfer coefficient for the same ambient temperature simulation is recovered by revisiting the free energy definition. The computational implementation has been utilized by using an automatized generation of relations with the aid of symbolic differentiation of the assumed free energy. Therefore, any more advanced formulation is easily implemented into the framework developed herein. For demonstrating the strength of the numerical implementation as well as comprehending the material model's capabilities, three different simulations have been utilized: 1. Two evolution equations, model A and B, are compared by means of the released heat in order to present the kink in the evolution denoting the glass transition temperature. 2. A realistic scenario has been used for demonstrating the vitrification in the case of an uncontrolled ambient temperature like the daily alternating temperature in Brussels. 3. Superposed hardening and mechanical loading is simulated during the cross-linking related glass transition change leading to a softening. This experimentally observed interplay between hardening and softening because of thermal and chemical contribution is possible by the proposed multiphysics simulation by adding glass transition dependence to the Young's modulus. Additionally, the role of the heat convection across the boundaries has been discussed.
We stress that the code is publicly available in [73] for encouraging its use under GNU Public licensing [74].
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/.