Multiphysics modelling of the mechanical properties in polymers obtained via photo-induced polymerization

Photopolymerization is an advanced technology to trigger free radical polymerization in a liquid monomer solution through light-induced curing, during which mechanical properties of the material are significantly transformed. Widely used in additive manufacturing, parts fabricated with this technique display precisions up to the nanoscale; however, the performance of final components is not only affected by the raw material but also by the specific setup employed during the printing process. In this paper, we develop a multiphysics model to predict the mechanical properties of the photo-cured components, by taking into account the process parameters involved in the considered additive manufacturing technology. In the approach proposed, the main chemical, physical, and mechanical aspects of photopolymerization are modelled and implemented in a finite element framework. Specifically, the kinetics of light diffusion from a moving source and chain formation in the liquid monomer is coupled to a statistical approach to describe the mechanical properties as a function of the degree of cure. Several parametric examples are provided, in order to quantify the effects of the printing settings on the spatial distribution of the final properties in the component. The proposed approach provides a tool to predict the mechanical features of additively manufactured parts, which designers can adopt to optimize the desired characteristics of the products.


Nomenclature a
Thickness of the photopolymerized part A Attenuation coefficient or absorbance A a b s , A p o l , A mon Absorption due to photoabsorbers, monomer converted into polymer and un-polymerized monomers, respectively A,A e Depletion matrix of the discretized domain and the corresponding finite element one, respectively b Length of Kuhn's segments C I Concentration of photo-initiator molecules (β) Engineering strain tensor φ(r), φ 0 (r) Dimensionless distribution function of the chains' end-to-end vector and the corresponding one at the initial stress-free state, respectively μ, μ Shear modulus of the current photopolymerized material and that of the fully-cured material, respectively ψ Energy stored in a single polymer chain Ψ Energy per unit volume of polymer ρ(r) Distribution function of the chains' end-to-end vector ϱ Degree of cure (or degree of conversion, DoC) achieved during the photopolymerization σ Cauchy stress θ Initiator molar absorptivity [■] Concentration of the chemical species represented by ■ e ■ Nodal values of the generic quantity ■

Introduction
In recent years, additive manufacturing (AM) has revolutionized the modern industry by introducing a new concept to fabricate complex geometries by means of a threedimensional model data [1][2][3]. Differently from traditional methods based on material subtraction, in AM technologies, parts are built layer-by-layer by joining successive layers of material on top of each other. While in its early times AM was mainly employed for prototyping, nowadays, it is quickly expanding to cover a wide range of industrial sectors and daily life products [4], ranging from aerospace [5], robotics [6], automotive [7,8], construction [9], and healthcare sectors [10]. Additive manufacturing technologies based on photopolymerization are among the most widely used for polymeric materials, enabling the fabrication of parts with high resolution, and are nowadays exploited for several applications in advanced material science [10,11], including structured materials [12][13][14], stimuli-responsive materials [15], and bioprinting [16]. Photopolymerization is a chemical-physical process which converts a liquid monomer solution to a solid threedimensional polymeric material, by applying UV light in a spatially controlled way according to the component's CAD model [1]. Among the photopolymerization AM processes, the earliest technique is stereolithography (SLA) [17]. As illustrated in Fig. 1a, the resin is irradiated point-by-point from a moving laser (moving in the horizontal plane xy), which along the build direction (here denoted with z) acts from the top (SLA top-down system) or from the bottom (SLA bottomup system). More recent technologies, such as digital light processing (DLP), directly project the UV light over a whole portion of the resin surface, allowing layers to be cured in a single shot (Fig. 1b). Similarly, the continuous liquid interface production (CLIP) employs a bottom-up photopolymerization method to achieve continuous printing through the creation of a permeable window that mitigates the detrimental effects of oxygen inhibition [18].
Despite the wide and increasing use of such techniques, a critical issue in the AM field deserves a specific attention: the performance of printed components is determined not only by the raw material but also by the specific setup employed during the manufacturing process, i.e., the parameters used to print the component itself [19]. As a consequence, components sharing the same geometry and raw material can have totally different final properties by changing the process parameters. In order to ensure that 3D-printed parts can be used reliably and according to the desired features, designers should be able to predict the mechanical properties in relation to the process parameters used for printing. To pursue this aim, traditional methods are based on empirical relationships between a specific mechanical property and process parameters, derived by testing several components printed with different setups [20,21]. Although such an approach is well suited to visualize trends of parameters-properties correlation, thus being a valuable strategy for quality control in the manufacturing process [22], it does not provide a physical understanding of such trends. Furthermore, empirical correlations are limited to the specific scenario (in terms of geometry, raw material, environmental conditions during printing, etc.) and cannot be extended to more general cases.
Methods based on a computational implementation of the theoretical description of the phenomena involved in the printing process are recommendable. These approaches provide a description of the actual chemical-physical mechanisms of the printing process and accurately predict the properties of the final part. With particular reference to AM photopolymerization, the first theoretical works were based on an energetic approach, according to which a component is assumed to be cured when the energy supplied to the resin reaches a critical value, a property of the raw liquid resin itself [23,24]. These methods are typically limited at predicting the liquid-solid transition, whereas they do not consider the evolution of the polymer network in time. As a matter of fact, all the components are in a solid state after printing but can have developed radically different properties due to diverse evolutions of the chemical species. Kinetics models have been proposed to describe the whole photopolymerization process, by means of a set of first-order partial differential equations describing the evolution of the reactant variables involved in the curing process [25][26][27][28]. Therefore, these models do not simply assess the transition to the solid state but they can also provide a quantitative measure of the converted monomer, by means of the so-called degree of cure [27].
The conformation of the polymeric network evolves according to the photopolymerization process, in turn affecting the mechanical behavior of the components [29]. In this paper, we develop a multiphysics model to describe the main physical, chemical, and mechanical aspects involved in photopolymerization, by simulating the whole printing process starting from the liquid monomer resin and arriving to the final solidified component. Firstly, the light emitted by a moving source, as in a real stereolithography setup, is modelled by means of a generalized Beer-Lambert law describing the light diffusion in a semi-transparent medium. Secondly, the evolution of the polymeric network is described according to the rate equations of the free radical polymerization reactions [26]. Finally, we employ a micromechanical approach rooted in the network's chain statistics [30] to infer the mechanical behavior as a function of the degree of cure.
After developing the theoretical basis, the model is implemented in a finite element (FE) framework, capable of computing the spatial distribution and evolution in time of the chemical species, and therefore of the mechanical properties of the final material. The proposed framework can accurately model the chemical-physical process leading to the final product, irrespectively of the specific photopolymerization technology, in order to be used by designers to predict and optimize the desired mechanical features.
The present research provides a multiphysics approach, involving the light diffusion, chemical kinetics and mechanics coupling, aimed at quantitatively evaluating the physical and mechanical characteristics of photopolymerized materials. It enables to control the photopolymerization parameters to precisely tune the final properties of the solid material; this aspect is very important, for instance, in additive manufacturing processes based on such a polymerization technique since the desired quality and safety level of the final polymer, according to the application in turn, can be guaranteed and optimized.
The paper is organized as follows: Section 2 presents the multiphysics modelling of photopolymerization. Section 3 is devoted to illustrate the FE implementation of the coupled problem of light diffusion and kinetics evolution of the chemical species involved in the photopolymerization process within the medium. Section 4 shows the results of parametric analyses in order to capture the effects of the main process parameters on the final component, while Section 5 is devoted to concluding remarks and future perspectives on this topic.

Multiphysics modelling of photopolymerization
In this section, representing the core of our work, we provide a multiphysics approach in order to model the complex phenomena of photopolymerization, which transforms an initial liquid monomer into a solid component with its mechanical properties being dependent by the setup of the photopolymerization process itself. Specifically, In Section 2.1, we illustrate the adopted kinetic model for the prediction of the evolution of the chemical species involved in the photopolymerization process, triggered by a UV-light source, with the final aim to quantify the active chains (i.e., those fully connected to the network and taking part to the load-bearing mechanism of the material) concentration. In Section. 2.2, we explain the governing equations of the light propagation within the medium, which is responsible for the chemical species evolution, i.e., of the photo-induced polymerization. Finally, in Section. 2.3, we present the main features of the proposed mechanical model, emphasizing the role played by the active chains concentrationwhich represents the main output of the photopolymerization process linked to the mechanics of a final component.

Evolution of chemical species
Photopolymerization is a chemical-physical reaction where a UV light triggers free radical polymerization [26]. At the first step of the reaction, the liquid resin monomer M, whose concentration in space and time is indicated with C M (X, t), is irradiated by the UV light. The light irradiation allows photo-initiator molecules Ph I , initially placed inside the resin monomer with a concentration C I (X, t), to be converted into free radicals R • , whose concentration is herein indicated with C R (X, t). During the second step, free radicals react with monomer molecules, providing the activation of the functional groups P • . From a physical viewpoint, this functional group simply represents a polymer chain which behaves like a snake, growing into the medium by reacting with other monomer molecules. The propagation proceeds until a termination stage is achieved: this occurs when either the polymeric chain meets a free radical (P • þ R • → k t P dead ) or it binds to another chain encountered along its growth path P • þ P • → k t P dead À Á . A sketch of the photopolymerization process, along with the reaction scheme, is reported in Fig. 2, where β, k p , and k t are the rate constants characterizing this chemical process.
As the chemical reaction proceeds, the monomers in the liquid solution are gradually consumed by combining with radicals and forming polymer chains, whose density progressively grows providing a stiffness increase of the material. A parameter typically used in the kinetic description of photopolymerization is the degree of cure, defined as [27]: which quantifies the amount of the monomer molecules converted into polymer chains. In the initial stage of the process (we refer here to the reference configuration of the domain whose points are identified by the position vector X), where the resin is in a liquid state, the degree of cure is ϱ(X, t = 0) = 0, since the concentration of the monomer molecules is equal to the initial one. As the reaction proceeds, ϱ(X, t) increases in time because C M (X, t) < C M (X, t = 0), i.e., the number of monomer molecules reduces due to the polymer chains growth (Fig. 2), and tends to 1 when the chain formation process is almost complete. The degree of cure is of crucial importance since it can be related to the concentration of active chains c a (X, t) (number of chains per unit volume), i.e., chains capable of transmitting a force being connected at both ends to other neighboring chains. The elastic parameters of the polymer network, and consequently its stiffness, depend on such a concentration. Indeed, as shown in [27] and other previous works [31], the concentration of the active chains after the gelation point, i.e., when the material can be considered to become stiffer and stiffer as the photopolymerization proceeds, can be related to the degree of cure by means of an exponential relation, defined as [27]: where μ is the shear modulus, E c , E d , and s are fitting parameters, to be determined under the assumption of incompressible behavior. The concentration of the active chains within the medium is assumed to depend point-by-point by the degree of cure in the continuum, which is also point-wise dependent on the UV light intensity I(X, t) and its duration.
In the specific process of photopolymerization, light is irradiated over the resin surface and diffused within the layer thickness, allowing photo-initiators to be decomposed into active species. The evolution of the photo-initiators (which are consumed in order to provide activation of free radicals according to Ph I → β 2R • ) can be modelled by defining the time rate concentration: where■ stands for the time derivative of the quantity ■, i.e., ■ ¼ ∂■ ∂t , while β represents the photodecomposition rate, which is a property of the photo-initiators. After such an activation, free radicals evolves triggering the initiation, propagation and termination of polymer chains [27]. The evolution of the free radicals can be expressed as: where m is the number of radicals generated in the photodecomposition (for instance, it can be assumed m = 2). k t is the termination rate, expressed as a function of curing in the form k t (t) = k t0 · g(ϱ), where k t0 is the initial termination rate constant and g is the curing-dependent term related to the diffusion controlled termination. The resulting problem is nonlinear. As it can be noticed from Eq. (4), the concentration of radicals increases due to photo-initiator decomposition while it reduces because of the termination mechanisms of propagation. Additional mechanisms, such as oxygen inhibition and other related chemical problems [32][33][34] , are here neglected for the sake of simplicity, but could be taken into account by slightly modifying Eq.(4) (see section 4.4).
As the combination with radicals proceeds, the monomers in the solution are gradually consumed, allowing the polymer chains initiation and propagation to occur (R • þ M → k p P • and P • þ M → k p P • ). Such a mechanism can be described by the following differential equation: where k p is the propagation rate constant that depends, analogously to k t , on the curing evolution itself [27]. The system represented by Eqs. (3)-(5) can be solved by assigning the proper initial conditions known from the liquid resin properties, i.e., C I (X, t = 0), C R (X, t = 0) and C M (X, t = 0), together with the knowledge of the rate constants (both the initial values and evolution equations are needed to solve the problem), and with the light intensity I(X, t) within the continuum (see Sections. 2.2 and 3.1). As discussed above, the main aim of kinetic models is the evaluation of ϱ(X, t), in order to subsequently compute the concentration of active chains c a (X, t) which is responsible for the mechanical properties of the polymer at the mesoscale (see Section. 2.3).
Interestingly, at a given position in the continuum, the problem of defining the evolution of the mechanical properties of the material can be directly related to the increase of active chains' concentration caused by the curing process. By applying the time derivative to both members of Eq. (1), we geṫ by exploiting the monomer molecules evolution reported in Eq. (5), we obtain: The concentration of the active chains evolves according to an exponential function, see Eq. (2), which can be reformulated in this form: where it has been assumed that ϱ gel ≅ 0, while μ is the shear modulus of the fully cured material (i.e., for an ideal curing time t → ∞). We assume that a small number c a0 of active chains exists when the material is in the initial (nearly) liquid state, i.e., c a (X, t = 0) = c a0 . The knowledge of the current value of the degree of cure, obtained by solving the differential Eqs.
(3)-(6), suffices to evaluate the corresponding current chain concentration value by using Eq. (7). The concentration of active chains at the time t c (curing time) can be evaluated as 0ċ a X ; t ð Þ dt . Further, according to the previous assumption, the time derivative of the chain concentration can be expressed as c˙a , where α is a model parameter. The chain concentration at the end of curing can be finally obtained as follows: which emphasizes how the final concentration of the polymer chains can be controlled by properly driving the time rate of the degree of cure during the photopolymerization process.

Light diffusion in a partially solid material
As shown In Section 2.1, the evolution of the chemical species involved in the photopolymerization process is triggered by a light radiation provided by a UV-light moving source. The light radiation is absorbed by photoinitiators, photoabsorbers, and by other light reactive species present in the matrix. Light diffusion within the material is attenuated and the variation of its intensity is a function of the concentration of the absorbing species and the molar absorptivity. In order to quantitatively describe how the light intensity depends on the spatial position in the material, the well-known Beer-Lambert law can be used [29,35]. Let us consider the large deformation regime for which, in the reference undeformed configuration (having domain and boundary Ω 0 , ∂Ω 0 , respectively), the generic point is indicated with X, while in the current deformed configuration (having domain and boundary Ω, ∂Ω, respectively) the corresponding point is identified by x. In a one-dimensional setting, where the light beam is assumed to impact perpendicularly the surface of the material (along the vertical Z axis), the diffusion equation can be easily obtained. At a given time instant, it can be assumed that the radiant flux I(Z) of the light outgoing from a slice of material lying perpendicularly to the light beam is reduced-with respect to that of the incoming flux-by the amount dI( In this simple one-dimensional case characterized by the linear coordinate z, if a constant value for the term A is assumed, at the generic time t and with the boundary condition I(0, t) = I 0 (t), the solution of the Beer-Lambert equation reads: The above expression can be generalized to a threedimensional setting, leading to well-known Beer-Lambert, a first-order linear partial differential equation (PDE): where Ω 0 , ∂Ω 0 indicate the domain of the problem and its boundary in the undeformed configuration, respectively, l(X, t) is the unit vector of the incoming light beam, I(X, t) is the light intensity at the position X and time t, ∇ is the gradient operator, and A(X, t) is the local attenuation coefficient (or absorbance) quantifying the light attenuation due to photons that did not travel through the material because of scattering or absorption. The local attenuation coefficient A(X, t) can be written as follows [27]: where θ is the photo-initiator molar absorptivity, C I is the photo-initiator concentration, A abs (X, t), A pol , A mon are the absorption due to photoabsorbers, converted polymer, and unpolymerized monomers, respectively. It is worth noticing that Eq. (10) is highly nonlinear since the term A(X, t) depends on the light-intensity history I(X, 0 ≤ t ≤ τ) at the location X.

Micromechanical model for polymers
Polymeric materials, such as rubbers, gels, etc., have an amorphous microstructure constituted by an entangled network of long linear chains, coiled and reciprocally joined at discrete points termed as cross-links. The microstructural state of polymers, characterized by thermal fluctuations, and their mechanics, are dominated by the entropic energy of the network, at least for not too large deformation regime [36]. The consolidated and well-known freely jointed chain model [37] -assuming the chains to be made of N rigid segments of equal length b, randomly arranged and connected at their extremities with neighboring chains-is usually adopted to describe the network topology, while the chain's end-to-end vector r is used to define the deformed state of the chain itself [38]. Due to the disordered nature of the network chains, a statistical description of the endto-end vectors is suitable to define the chains arrangement within the network. To this end, a distribution function ρ(r) can be conveniently introduced [30], providing the number of network's chains with a given end-to-end vector r = (r x , r y , r z ) in a given state of the material, for instance the stress-free one. For our purpose, it is convenient to write the above-mentioned distribution in the form ρ(r, t) = c a · φ 0 (r, t), where c a corresponds to the number of reciprocally connected load bearing chains per unit volume (see Sects 2.1 and 2.2), while φ 0 (r) is a dimensionless distribution (usually obeying the standard Gaussian law with mean value r = 0 and standard deviation b ffiffiffiffiffiffiffiffiffi N =3 p [38]). A further assumption is the so-called affine deformation hypothesis, stating that the deformed chain's length is given by the chain's mean square length and λ the macroscopic stretch of the material. According to the above definitions, within the time interval 0; t ð Þ (note that this time interval follows the curing time interval (0,t c )), we have: where the operator 〈■〉 indicates the following integral evaluated over the chain configuration space, sinθdθdβ, with θ, β being the Euler angles of the end-to-end vector r in the 3D space [30,39]. The mechanical response of the polymer is assumed to take place after the completion of the photopolymerization process, i.e., the micromechanical model enters into play only after printing the material. It is worth noticing that the above definition of the chain concentration assumes that only the cross-linked chains enter into c a , i.e., only the chains linked at their extremities with other chains have to be considered for the mechanical properties of the material. The knowledge of the distribution function ρ(r, t), allows us to evaluate the elastic energy density stored in the material once the deformation energy of a single chain ψ(r) is expressed through its end-to-end vector r. According to the above statement, by adding up all the energies coming from the chains contained in a representative volume of material the energy density can be evaluated as: It is worth recalling here that the derivative with respect to r of the energy per single chain gives the force f existing in the chain, i.e., f(N) = ∇ ψ = ∂ψ(λ, N)/∂r, where the dependence of ψ on the chain length N has been emphasized. In general, being ψ a non-negative function of |r|, the energy per unit volume of Eq. (13) is nonzero also in the stress-free state of the material. It is thus much useful to define the deformation energy density of the material as is the dimensionless distribution function at t = 0 that is assumed to be in a macroscopically undeformed state (at the end of the photopolymerization process, i.e., when the material is assumed to be in the stress-free state), i.e., when F = 1 (being F = ∂x/∂X the deformation gradient tensor), while φ(r, t) is the corresponding distribution in a generic deformed state of the material induced by mechanical actions. The stress state in the material can be finally be determined as follows [30] : where the true Cauchy stress σ tensor has been used as the stress measure, while π is the hydrostatic pressure (playing the role of a Lagrange multiplier to enforce the proper volume change ratio J = det F accounting for the actual volume change of the material), L is the velocity gradient tensor. For an incompressible material J = det F = 1 and tr L = 1 : L = 0 so the above expression simplifies.
As it can be noticed, the above-described micromechanical model allows determining the polymer's mechanical properties and thus its response under a mechanical action; this can be done once the chain statistical distribution φ 0 (r, t) and the chain concentration are known at any point of the solid domain at the end of the polymerization time, t = t c . The end of the photopolymerization process is assumed to occur before the mechanical problem takes place. In our case, at a given point inside the initial liquid monomer material, the chain concentration is determined by the chains formation triggered by the light energy; the distribution of the chain end-to-end vectors at the formation time can be assumed to follow the Gaussian function, i.e., ρ 0 (r, t = t c ) = c a · φ 0 (r, t = t c ) being φ 0 the standard Gaussian, while the chain concentration is obtainable through the solutions of Eqs. (6) and (8), by using the kinetics of chains formation governed by Eqs. (3)-(5). Finally, according to Eq. (14), once the current chain distribution φ(r, t) determined according to the deformation process applied to the material has been determined, the stress state can be easily evaluated [30].
It is worth recalling that in this work, we apply the abovepresented multiphysics model for studying the tunability of the mechanical properties of photopolymerized materials (see Sections. 2.1 and 2.2 for the theory, Sections. 3.1 and 3.2 for the FE implementation, and Section. 4 for the numerical simulations); the study of the mechanical response related to the properties of the photopolymerized polymer, will be considered in future studies.
3 FE implementation of the coupled problem of light diffusion and chain growth evolution

Weak form of the light diffusion governing equations and its numerical implementation
In order to determine a form of the light diffusion governing equation suitable to be implemented in a computational code, the development of a weak form of the mathematical problem is usually adopted. Such a form can be obtained by introducing a proper functional or by transforming the strong form represented by Eq. (10) into an integral (weak) formulation. By multiplying the Beer-Lambert equations by the weight function w(X) and integrating over the domains where they are defined, yields to: By assuming that the vector field l(X, t) is constant in space (i.e., the incident light beam is made of parallel rays) and depends only on time t, i.e., l X; t ð Þ ¼ l (t), and by applying the divergence theorem to the first term of Eq. (10) we get: n being the outward normal to the boundary of the domain, while f(X, t) represents the light-intensity distribution over the irradiated boundary of the domain. The term on the right hand side of Eq. (16) can be rewritten as follows: where l n ¼ l À t ð Þ Á n is the projection of the light unit vector of the incoming laser beam on the normal to the boundary surface of the domain. By referring to Fig. 3, the incident light beam forms an angle γ with the normal to the outer top surface of the domain; by assuming γ = 0, we have l n = 1. The above integrals represent the weak form of the governing equation; they can be conveniently rewritten by introducing the domain discretization in which the field variables are interpolated through the corresponding nodal values, i.e.: where the subscript h indicates the values interpolated from the corresponding nodal quantities, whose values are indicated with e ■, being n n the number of nodes of the element; e I i ; e w i are the nodal light intensity and the nodal values of the weight function, respectively, while ∇ X represents the gradient operator in the reference (or undeformed) configuration.
In Eq. (10), the presence of only the convective (or transport) term, l · ∇ X I, leads to the so-called boundary layer, i.e., a region-generally arising close to the boundary-where the solution is characterized by large gradients or oscillations: this entails an unstable and meaningless numerical results [40]. These cases are typical of problems whose governing equations have a large Peclet number, defined as the ratio between the rate of advection of a physical quantity to the rate of diffusion of the same quantity. To overcome this problem, a stabilization term can be added to the original equation. This term is usually an artificial diffusion-like contribution term of the form Q hΔI(X, t), where Q is a coefficient and h a characteristic small geometrical size related to the discretized problem in turn [41]. The amended stabilized PDE governing the problem can thus be rewritten as follows: The introduction of the discretized field quantities of Eq. (18) in the weak form of Eq. (16), written for a generic finite element e having volume V e and boundary S e and thanks to the arbitrariness of the weight function w(X) (or, correspondingly, of its nodal values e w ), leads to the problem written in the following matrix form: where E e is the light gradient matrix, A e is the depletion matrix of the finite element, D e is the stabilization matrix (being Q > 1, h e a coefficient and a characteristic size of the finite element, respectively), while P e (t) is the vector of the nodal values of the distribution of the incoming light intensity on the boundary and f n (X, t) = l n f(X, t) represents the distribution of the of the beam light intensity normal to the boundary of the domain. Upon assembly of the elements' matrices over the whole discretized P e with A the assembly operator and ne the number of finite elements), we get: E t ð ÞþA t ð ÞþD ð Þ e I t ð Þ¼P t ð Þ. For the sake of simplicity, we do not consider the effect of light refraction since typically the light is irradiated by the photo-curing apparatus in a direction perpendicular to the surface of the monomer fluid; in cases where the light is not normal to the surface, the light propagation angle may change in space due to the difference of the refractive index, and the propagation direction must be determined accordingly.
In a 2D setting, the above discretized system of equations becomes: being a the element thickness, h e can be set as h e ¼ ffiffiffiffiffi A e p , being A e the area of the finite element e, while L e is the length of the finite element edge irradiated by the incoming laser light. The Fig. 3 (Color online) Scheme of the light irradiation of the top surface of the material to be cured, whose domain is discretized through finite elements. The distribution of the light intensity on the irradiated top surface is assumed to be distributed according to a Gaussian function which moves with a predefined speed v along the horizontal direction X. The degree of cure taking place in the material depends on the light intensity and its permanence in time and leads to a specific concentration c a of crosslinked chains scheme of the light diffusion in the medium, upon irradiation of its top surface by the laser beam source, is illustrated in Fig. 3.

Computational kinetics modelling of chains' growth
Once the light-intensity distribution has been evaluated by solving the system of Eq. (20) with the corresponding boundary conditions, such a quantity is used to quantify the degree of cure taking place at each Gauss point of the finite elements used to discretize the monomer vat domain. Firstly, the interpolated light-intensity value at the Gauss point is obtained through the standard expression I h;GP t ð Þ ¼ ∑ n n i¼1 N i;GP e I i t ð Þ, where N i, GP indicates the shape function related to node i and evaluated at the specific Gauss point GP, while 0 ≤ t ≤ t c is the generic time instant. The current value of the photoinitiator concentration is then evaluated by the time integration of the corresponding time rate, i.e.: where Δt = (t − t −1 ) is the time step interval adopted for the integration over the time domain. Once the photo-initiator concentration is known, the concentration of the free radicals can be determined by firstly evaluating its time rate: Due to the polymerization process, the concentration of the initial liquid monomer (referred, as done before, to the Gauss point positions GP) decreases in time at the rate: Finally, the degree of cure can be easily evaluated as follows: and the related chain concentration density can be determined as c a;GP t ð Þ ¼ c a;GP Á exp α ϱ GP t ð Þ−1 ð Þ ½ ¼c a;GP Á p GP t ð Þ ð26Þ being c a;GP ¼ c a;GP t ¼ ∞ ð Þ¼ μ k B T while 0 < p GP (t) ≤ 1 is the coefficient quantifying the degree of developed chains, with respect to the fully cured material obtainable for an exposition to the light source for a theoretically infinite time.

Numerical simulations
This section is devoted to the simulation of the photopolymerization process of a liquid monomer upon irradiation by a laser beam. Specifically, we mainly aim at illustrating the effects provided by: the laser beam translation velocity (Section 4.1), the peak light intensity of the laser beam (Section 4.2), and the material's light absorbance (Section 4.3), in order to understand their role in quantifying the mechanical characteristics of the obtained cross-linked polymer. In all the following examples, studied as 2D problems, the light intensity irradiated on the top boundary surface of the domain is assumed to be distributed according to the equation: where v is the linear constant velocity of the light source, while c defines the width of the assumed Gaussian distribution for the light intensity (Fig. 3).
In all the examples, we simulate the curing process of a liquid monomer plane vat with height h = 10 mm (usually known as layer thickness in the realm of AM) and width w = 30 mm. The laser moves on the top of the liquid monomer vat from the left to the right. Unless differently stated, in the following, we assume the following laser printing parameters: v = 10 −3 m/s, c = 10 −6 m, and I 0 = 150 W/m 2 . Moreover, the initial chemical properties of the liquid monomer being solidified are assumed as follows [27]: (i) initial photo-initiator concentration C I (X, 0) = 20 mol/m 3 ; (ii) initial liquid monomer concentration C M (X, t = 0) = 3000 mol/m 3 ; (iii) photo-decomposition rate β = 8 · 10 −4 s 2 /kg, (4) propagation and termination rates k P = k T = 0.21 m 3 /mol · s. For the sake of simplicity, we consider a constant value of the chemical rates k P and k T . Finally, unless explicitly stated, in the following we assume a constant material absorbance equal to A = 600 m −1 . It is worth mentioning that, since the aim is to perform a parametric study to physically quantify the effect of the curing process characteristics on the printed polymer, the parameters adopted are not related to a specific commercial liquid monomer used in 3D printing, but the order of magnitude of the adopted values are within the range of realistic cases as can be found for instance in [25,27]. To conclude this section, a real-case simulation is finally reported In Section 4.4.

Effect of laser beam translation velocity
In this section, we consider the effect of the laser beam moving speed on the degree of cure and, consequently, on the chain concentration related to the polymer's shear modulus. We consider two different cases, with laser speeds equal to v = 10 −3 m/s and v = 10 −4 m/s, respectively. It is assumed that the curing process ends when the laser covers the whole width of the layer being printed: the time required by the laser to cover such a distance, known in the AM literature as curing time, is defined as t c = w/v. A dimensionless curing time can conveniently be defined as t * = t/t c , with 0 ≤ t * ≤ 1. The geometrical dimension of the layer is kept fixed between the simulations (w = 30 mm), so that the curing times are equal to t c = 30 s and t c = 300 s for the two velocities considered. Furthermore, the material absorbance is assumed to be constant, so that the light-intensity field within the material is the same, irrespectively of the laser speed. In other words, the light intensity distributions shown in Fig. 4, are the same for both the laser speeds adopted.
For each laser speed, we are interested in determining the characteristics of the material being printed at two different time instants, namely t * = 0.5 and t * = 1, corresponding to two distinct positions of the laser with respect to the layer being printed. The first case corresponds to the laser beam placed at X = w/2, see Fig. 4a showing the light-intensity field within the material at that instant. The second case corresponds to the final position given by X = w, see Fig. 4b showing the light-intensity field within the material at the end of the curing process.
The corresponding distributions of the degree of cure are shown in Fig. 5. In particular, Fig. 5a illustrates the degree of cure within the material when the laser is placed at the middle of the layer being printed for the case characterized by the highest laser speed (or, correspondingly, for the lowest curing time), while Fig. 5b shows the degree of cure when the laser is at the same position but reached with the lowest laser speed (or, correspondingly, for the highest curing time). Since the degree of cure is the results of a kinetic equilibrium (see Section 2.1), a longer irradiation entails a more pronounced, homogenous and deeper degree of cure within the material. Analogous considerations can be made by looking at Fig. 5c and d, representing the distribution of the degree of cure at the end of the curing process. Figure 6 shows the dimensionless shear modulus, μ t ð Þ= μ À within the analyzed domain. Such a quantity can be obtained through the relationship with the chain concentration, Eq. (7), leading to the expression μ t ð Þ ¼ μ Á exp α ϱ x; t ð Þ−1 ð Þ ½ , where μ ¼ 267 MPa is the shear modulus for the fully cured polymer (this value is always assumed henceforth) and α = 3; in other words, the fields of the shear modulus distribution are a straightforward consequence of the degree of cure. As it can be appreciated from Fig. 6c, at the end of the curing process with the highest laser speed ( other parameters being fixed), roughly only one half of the layer is solidified. Moreover, the stiffness of such a cured part is not homogenous, showing a higher shear modulus close to the top boundary, while a nearly not-cured zone (μ(t) → 0) appears at the bottom. On the other hand, the curing process with the lowest laser speed allows to get a more pronounced, homogeneous and higher values of the shear modulus within the domain, Fig. 6d, even if there is still a zone with a low shear modulus value close to the bottom of resin vat. The previous considerations allow us to appreciate the importance of a layer-by-layer process in 3D-printing.
Assuming that the component to be printed had the dimension of the liquid vat investigated in this case (h = 10 mm and width w = 30 mm), it would require to be sliced into layers of thickness d < h in order to be optimally solidified. The curing process, here considered for a single layer for sake of simplicity and clearness, must be repeated for all the layers in which the component being printed has been divided. Firstly, the laser solidifies the top layer, irradiating the liquid monomer placed at Y = d (in the presented case one single layer is adopted, i.e., d = h) by moving rightwards in order to cover the whole part width; then, the part already printed is covered by a liquid layer of thickness d and the laser repeats the curing operation until the whole width is covered again. Such an operation must be repeated until the whole height h of the component is achieved, curing successive layers into which the total height has been subdivided. The number of layers (and therefore the proper layer thickness for the curing process) can be evaluated by using the present results with respect to the optimal distribution of the shear modulus. For instance, the optimum setup of printing parameters can be determined in order to achieve μ t c ð Þ=μ→1 in each material point of the component in the shortest possible curing time, if such a component is desired at the end of the printing.
In order to generalize the present case, it is worth considering what happen if the initial liquid monomer resin would be irradiated by two laser sources at the same time, one operating from above and the other from below. The light intensity at a point of the material can be evaluated through the superposition of the light intensities induced by both sources. However, since the absorbivity generally depends on the degree of conversion of the material (see Sections. 2.2 and 4.3), the lightintensity superposition from the two sources cannot be applied to get the whole degree of cure because the light diffusion process must account for their reciprocal interaction. For not too low laser speeds, this provides a final photopolymerized

Effect of light intensity of the laser beam
In the following examples, we investigate the effect of the laser beam light intensity, by keeping all the other parameters fixed. Two values of the peak light intensity are considered, namely I 0 = 25 W/m 2 and I 0 = 100 W/m 2 . The light-intensity distribution when the laser reaches the middle position of the layer being printed (X = w/2) is shown in Fig. 7a and b, for the two distinct light intensities considered, respectively. As expected, the case corresponding to the higher peak light hitting the vat surface, provides a more pronounced, diffused, and deeper light-intensity field within the medium.
In turn, the curing of the material is significantly affected: Fig. 8 illustrates the degree of cure obtained for the two adopted peak light intensities, corresponding to the laser position X = w/2. The curing process with the lowest peak light intensity (Fig. 8a, A1) provides a less pronounced and shallower degree of cure within the medium, in comparison with the curing process characterized by the highest peak light intensity (Fig. 8a, A2). The corresponding profiles of the resulting degree of cure obtained along the two straight vertical cross-sections 1 and 2 shown in Fig. 8a, are reported in Fig. 8b.
The dimensionless shear modulus distribution, corresponding to the degree of cure shown in Fig. 8, is reported in Fig. 9a for the two adopted peak light intensities. Close to the top boundary, the material that is being cured with the lowest I 0 presents values of the dimensionless shear modulus limited to 0.5, while it achieves greater values if cured with the highest light intensity  Fig. 9a, are reported in Fig. 9b. It should be noticed that, an increase of the peak light intensity provides a stiffer material at the same curing time (other parameters being fixed). However, it is important to highlight that a too large light intensity may induce a degradation of the material being printed [42,43]. This aspect is still an open issue in the AM research field and will not be considered in the present study.

Effect of the material's absorbance
In this paragraph, we investigate the effect of the material's absorbance on the curing process. Three different cases are considered. In the first case we model the curing process of a liquid monomer assumed to have a constant material absorbance, namely A = 600 m −1 , which has the magnitude order of a purely liquid monomer. In the second case we assumed the material absorbance still being constant in space, but with an higher value namely A = 2400 m −1 , which has instead the magnitude order of an already cured monomer (polymer in a solid state). However, real photopolymerization processes are characterized by a gradual transition from the liquid monomer to the solid polymer. Therefore, the absorbance ranges gradually between two extreme values according to the expression provided by Eq. (11). For this reason, a third case is considered, in which the curing process occurs with the absorbance varying in time according to the above-mentioned relationship. We assume the initiator molar absorptivity to be θ = 30 m 2 /mol, the absorbance of a fully cured polymer and of a liquid monomer (infinitely transparent) are, respectively, A pol = 2400 m −1 and A mon = 0, while no photoabsorbers are considered in the resin, i.e., A abs (X, t) = 0. From Eq. (11), by taking C I (X, 0) = 20 m 3 /mol, it can be noticed that when the material is in a liquid state (before curing, i.e., ϱ(X, 0) = 0) we get A = 600 m −1 , which is the value corresponding to the purely liquid material. On the other hand, when the material attains the fully cured condition, the photo-initiator concentration C I tends to zero while the degree of cure tends to one, leading to A = 2400 m −1 , which is the value corresponding to the fully solid polymer. In other words, in the first simulation, we assume to apply the curing process to a material that ideally maintains the absorbance of the original liquid monomer, while in the second case we ideally assume the absorbance to be constantly equal to that of the solid material; finally, in the third case, we adopt a more realistic hypothesis assuming the absorptivity to vary with the degree of cure. We firstly consider the same peak light intensity I 0 = 150 W/m 2 for all the examined cases and introduce a dimensionless intensity I 0 =Ī 0 that will be used in the following, being I -0 = 150 W/m 2 . The light-intensity distribution obtained by using a constant value of the absorbance are reported in Fig. 10a, b, for the lowest and the highest material's absorbance, respectively. The maps of Fig. 10 represent the results corresponding to the laser position X = w (i.e., at the end of the curing process, t * = 1.0); however, similar results can be obtained for other laser beam positions. As expected, the light penetrates deeper and deeper in the liquid monomer with the lowest absorbance, Fig. 10a, while it is more dimmed in the case with the highest absorbance, Fig. 10b.
The degree of cure distributions at the end of the curing process, characterized by the lowest and the highest value of absorbance, are reported in Fig. 11. Due to the higher intensity of the light in the material corresponding to the lowest A, a higher degree of cure is achieved with respect to that with the highest A (see Fig. 11a, A1 and A2, respectively). The corresponding profiles of the degree of cure evaluated along the three vertical sections 1, 2, 3 illustrated in Fig. 11a, A1, A2, are reported in Fig. 11b. It can be observed that ϱ presents a higher gradient in the Y direction for the liquid monomer with the highest A (black lines in Fig. 11b) while its gradient is less pronounced for the material with the lowest A (blue lines in Fig. 11b).
The corresponding dimensionless distributions of the shear modulus are reported in Fig. 12a, with the respective vertical profiles in Fig. 12b. As expected, the stiffness of the material decreases through the thickness according to the degree of cure. Figure 13a illustrates the degree of cure existing in the material at the end of the curing process, for a monomer with A varying in time according to the solid-liquid fractions in the material throughout the photopolymerization process. As expected, in this case the resulting field is in between those obtained for the lowest constant A (Fig. 11a, A1), and the highest constant A (Fig. 11a, A2). Finally, the time evolutions of the degree of cure and of the dimensionless shear modulus,  Fig. 13b, c and d, e, for all the cases investigated in this section. Taking the results corresponding to the varying A as reference, Fig. 13b, c shows that the simulation of the curing process with the highest A leads to an underestimation of ϱ, while the simulations performed by adopting the lowest constant A overestimates ϱ. It is worth mentioning that the more realistic simulation using a variable A requires a bigger computational effort due to the nonlinearity introduced by the relationship A = A(ϱ). Similarly, Fig. 13d, e reports the evolution of the dimensionless shear modulus during the whole curing process, leading to the same consideration made for the ϱ evolution. The achieved stiffness is underestimated by using the highest A while it is overestimated by using the lowest A, and it stays in the middle ground by means of A varying over time.
The three cases described above have been simulated also by using a lower peak light intensity, namely I 0 = 7.5 W/m 2 , i.e., for I 0 =I 0 ¼ 0:05. The evolutions of the degree of cure and of the dimensionless shear modulus during the curing process are reported in Fig. 14. As expected, the curves corresponding to a varying A fall in between the curves corresponding to the highest and the lowest constant absorbance values. However, compared to the analyses with the higher peak light intensity, the curve is now shifted towards the case with the lowest constant absorbance, Fig. 14.
In addition, the results highlight the potential of considering a constant material absorbance (liquid-like or solid-like) in the photopolymerization model, allowing to achieve a curing evolution which is close to the real one, i.e., with a variable absorbance. More specifically, the solid-like absorbance provides a good approximation for the curing process with higher light intensity, whereas the curve corresponding to the lower light intensity is better replicated by the liquid-like absorbance; however, it must be observed that the above conclusion depends also on the spatial position of the considered point within the cured material.
In other words, the use of the liquid monomer absorbance (the lowest one) rather than that of the fully cured polymer (the highest one) in predicting the curing evolution of a material, cannot be defined a priori but must be related to the photopolymerization velocity, and therefore depends on the peak light intensity and the travelling speed adopted in the process. The analyses performed allow comparing the quality of the photopolymerized material according to the simplifying hypotheses on the absorbance values mentioned above. In general, the assumption of the highest value for the absorbance (that of the solid material) underestimates the resulting mechanical properties of the printed elements, while the opposite occurs if the lowest value of the absorbance is used in the simulations. However, as discussed above, such under or overestimations have to be related to the light intensity, laser velocity and spatial-dependency; if the light beam is characterized by a sufficiently high intensity and the traveling speed is low enough, the adoption of the solid-like absorbance should provide reasonably results, especially for a small layer thickness .

A real-case simulation
In this last section, we investigate a real case of photopolymerization process as reported in the work by Wu et al. [27]. The monomer used is polyethylene (glycol)  Fig. 13 (Color online) Map of the degree of cure (a) at the end of the curing process for a material assumed to have an absorbance varying in time according to the degree of cure evolution, see Eq. (11) in the text. Evolution of the degree of cure (b, c) and the dimensionless shear modulus evolution (d, e), at the material points P (ξ = (h − Y)/h = 0.1) and M (ξ = (h − Y)/h = 0.5), respectively, during the curing process for all the cases investigated. In all the cases, the peak light intensity is kept fixed at I 0 = 150 W/m 2 , i.e., I 0 =I 0 ¼ 1 diacrylate (PEGDA) (with a molecular weight of 250 g/mol, Sigma-Aldrich, St. Louis, USA), and 0.3% weight percent 2,2-dimethoxy-2-phenylacetophenone (Sigma-Aldrich, St. Louis, USA) as the photo-initiator. We use the proposed FE model by assuming a stationary laser beam v → 0 m/s irradiating a single point lying on the resin surface, as made in [27].
The rate constants used for the simulations are given by k p ¼ k p0 1þ k p0 k p;D0 e cϱ and k t ¼ 1þ k p0 k p;D0 e cϱ for the propagation and termination rate constants, respectively. Such constants depend on the degree of cure, which is unknown; therefore, the step-by-step numerical procedure previously illustrated requires to update the rate constants at the end of each calculation step, before moving on to the next one. Several parameters are required to quantify the rate constants evolution. These depend on the resin property, environmental conditions etc, and need to be experimentally measured. Following Wu et al. [27], such parameters are: k p0 ¼ 1:86048 m 3 mol s , k p;D0 ¼ 8:994 Á 10 8 m 3 mol s , k t;SD ¼ 4:39 Á 10 3 m 3 mol s , k t;TD0 ¼ 10 024:43 m 3 mol s , β = 8.999 · 10 −4 s 2 /kg, c= 34.149 and C RD = 1.0146. Furthermore, the parameters adopted for the estimation of the chains concentration are: E d = 3.321 MPa, E c = 1.059 MPa, s = 5.248, and ϱ gel ≅ 0.15.
The results, expressed as usual in terms of the degree of cure and the polymer chain concentration, are reported in Fig. 15 for a monomer resin cured for t c = 100 s by using different light intensities. In the figure, the dashed lines correspond to the literature results [27] while the present results are represented by the continuous lines. The chain concentration shown in Fig. 15b has been made dimensionless with respect to that corresponding to the fully cured material, i.e., c a ¼ c a ϱ ¼ 1 ð Þ¼9:86 Á 10 27 m −3 . It can be noticed that the curing process as predicted by our model appears to be faster than that reported in the literature. This might occur because we are neglecting the role played by the oxygen concentration (indeed, the same assumption is often made in literature, see, for instance [28,44] ). In fact, the oxygen eventually present in the environment may inhibit the polymerization reaction by combining with radicals, thus hindering the propagation of the polymer chains. The specific issue is widely treated in [33]. As clearly demonstrated by our study, neglecting the oxygen role in a kinetic model provides an overestimation of the obtained degree of cure, resulting in a nonconservative prediction of the polymer's mechanical properties. Such a problem can also be overcome with specific strategies, by means of physical-methods (for instance, curing in an oxygen-free atmosphere thanks to physical barriers or inert gases) or using chemical methods (for instance, by inserting specific additives in the liquid monomer) as shown in [45] .
In conclusion, as already mentioned In Section 2.1, the chemical terms involved in a complex kinetic model like the one here presented, have to be slightly modified and adapted  (a, b) and of the dimensionless shear modulus (c, d), at the material point P (ξ = 0.1) and M (ξ = 0.5), respectively, for cases with constant and variable A investigated in this section by adopting a peak light intensity equal to I 0 = 7.5 W/m 2 , i.e., I 0 =I 0 ¼ 0:05 according to the actual condition of the problem under consideration which can vary from case to case, in order to correctly capture all the involved chemical-physics phenomena. For instance, the FE model presented herein can be easily updated in order to capture the role played by the oxygen inhibition: following [27], we modify Eq. (4) by adding the concentration of oxygen, namely C O (X, t), which affects the evolution of free radicals, i.e., The problem is now even more non-linear, since the production of oxygen depends on the evolution of free radicals themselves; in fact, according to [27] Þ which is the equation that completes the system. The curing process considered above has been re-simulated by taking into account the role played by the oxygen inhibition, for the exemplary case with I 0 = 2 W/ m 2 ; the results are plotted in orange solid squares in Fig. 15 (Model (*)). As can be clearly appreciated, by taking into account the oxygen inhibition, the trend of the degree of cure and of the chain concentration evolution, appear to be shifted (delayed in time) with respect to the same curves obtained in the case of oxygen inhibition neglected (continuous line with black squares). Such curves are in good agreement with those obtained in [27] for the same representative case (black dashed line with empty squares).
Finally, it is worth mentioning that in order to enhance the mechanical properties of the material after photopolymerization, additive manufactured components may be subjected to postcuring treatments; to this end several techniques are available, such as the use of an ultraviolet chamber, microwave, ultra-sound or conventional oven where the printed parts has to be placed for a while [46]. Post-curing has the main role of further improving the degree of cure of a material, by completing the polymerization process after printing [28] , which could have not reached the complete conversion during printing; this could be caused, for instance, because of the use of a too low light intensity [3].
The most common post cure techniques are the UV and thermal ones, both aimed at directly polymerize any uncured liquid trapped within the solid by radiation of the component with ultra-violet light or by promoting the development of polymer chains through higher temperatures, respectively.
In the present study, we have developed a model describing the photopolymerization printing process only, while all the possible post-treatments are beyond the scope of the research and have been intentionally neglected. A simplified phenomenological model for the post-curing treatment can be find in [47].

Conclusions
In the present paper we have proposed a theoretical approach, implemented in a finite element framework, for the simulation of the photopolymerization process used in additive manufacturing technologies, such as stereolithography (SLA), whose potentialities are nowadays of great interest because of its capability to fabricate objects at a very small scale.
The problem has a multiphysics nature, involving the phenomenon of light diffusion within a semi-transparent medium, the chemistry of the cross-links formation through the photopolymerization process and the description of the resulting mechanical properties of the final solid material. We have considered a moving light source, providing the energy for the polymerization to occur in the liquid monomer, and a micromechanical approachbased on the statistical distribution of the chains end-to-end vectorto quantify the mechanics of the material at the mesoscale level. The FE implementation of the multiphysics model has allowed to visualize the spatial time-dependent distribution of the degree of cure, and in turn of the mechanical properties, according to [27] Fig. 15 (Color online) a Evolution of the degree of cure and b of the dimensionless chain concentration, evaluated at the surface of a liquid monomer irradiated by a non-moving laser with different energy intensities. Continuous lines refer to model results obtained by intentionally neglecting the inhibition provided by the oxygen, differently by the results reported in [27] (dashed lines and cross symbols). Model (*) results obtained by accounting for the oxygen inhibition are reported for the representative case characterized by I 0 = 2 W/m 2 (orange squares) different simulated printing setups. It is worth noticing that this approach marks a remarkable improvement with respect to previous works (e.g., see [28]), where the spatial variation could not be described.
Several parametric analyses have been proposed in order to understand the role played by the main quantities governing the chains cross-linking phenomenon, such as peak light intensity, laser beam velocity (curing time) and material's absorbance. A real case simulation has also been considered, providing additional interesting considerations. For instance, we have verified the important role played by oxygen inhibition, which greatly delays the curing process. The approach proposed represents a comprehensive tool to simulate the complex photopolymerization process, providing a full control on the involved parameters, in order to optimize the design of AM components to meet the required mechanical performances.
Authors contributions RB conceived the research, designed numerical analyses, and wrote the paper. MPC designed and performed the numerical analyses and participated in the revision. LM, AS, and MT participated in the revision of the manuscript. All authors read and approved the final manuscript.
Funding Open access funding provided by Università degli Studi di Parma within the CRUI-CARE Agreement. The authors would like to thank the support from European Union's Horizon 2020 research and innovation programme (H2020-WIDESPREAD-2018, SIRAMM) under grant agreement No 857124.
Data availability There is no availability of data and/or materials.
Code availability The numerical analyses have been performed with an in house FE code that is made available by the Authors under request.

Declarations
Ethical approval The present paper contains an original research, it has not been published elsewhere, it has not been submitted simultaneously for publication elsewhere, and publication has been approved by all coauthors.

Consent to participate Not applicable.
Consent to publish The present paper does not require any consent to publish since all the figures, tables, and text are original.

Conflict of interest The authors declare no conflict of interest.
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/.