Rate-dependent energy dissipation of graded viscoelastic structures fabricated by grayscale vat photopolymerization

A 3D printer utilizing the gMSLA method primarily comprises a UV light source, an LCD chip, a vat of liquid photopolymer resin, and a print platform, as depicted in figure 2. During the printing process, UV light is emitted from the source onto the LCD chip, where grayscale masks are loaded to control the light intensity and thus the curing process on a pixel-by-pixel basis. Every mask comprises grayscale pixels, spanning from white with RGB $= (255,255,255)$ or G = 1, signifying full light intensity, to black with RGB $= (0,0,0)$ or G = 0, indicating zero exposure, i.e. keeping the photopolymer at the material point uncured. These masks illustrate the pattern to be solidified in a printed layer [31]. The UV light is filtered and directed by the mask onto the optic lenses, which control the transmitted area and focus on the tank of liquid resin.

Zoom In Zoom Out Reset image size

The photopolymer resin, composed mainly of free monomers and photoinitiators, undergoes a transformation under the UV light's influence, see figure 3. The photoinitiators decompose into active radicals, initiating a reaction with the monomers, converting them into active monomers. Subsequently, these active monomers react with other monomers, forming polymeric chains. As the process progresses, the activation and propagation of active monomers lead to an increase in the crosslink density, resulting in the formation of networks of molecules. Consequently, a solid-state polymer materializes in the designated domain, with the degree of cure depending on the incident light energy [35–37]. For additional insights into photopolymerization and cross-linking during printing, refer to [36, 38, 39] among others.

Zoom In Zoom Out Reset image size

Moreover, the solidification process and crosslink density play a crucial role in determining the mechanical properties of the polymer [40]. In figure 3, the fully crosslinked polymer network, depicted in green, represents the material's elastic behavior. Conversely, the orange regions in figure 3 represent the remaining free monomer chains, which are essential for the material's viscoelastic behavior. Consequently, the viscoelastic properties of the material are contingent on both the degree of monomer conversion and the density of the crosslinked network [33, 41].

The mechanical properties of viscoelastic materials can be effectively represented by rheological models [42–44]. The one-dimensional spring-dashpot presentation of a generalized Maxwell-type viscoelastic model that demonstrates relaxation characteristics is shown in figure 4. A set of springs, which respond linear elastic, is assumed to model the solid behavior. The stiffnesses of the free spring on one end and the spring for the so-called Maxwell element are determined by Young's moduli, denoted as $E_{\infty}\gt0$ and $E_i\gt0$, respectively, where i takes values from 1 to n and represents the number of Maxwell elements, which are connected in parallel to the elastic spring. The flow behavior is represented by a Newtonian viscous fluid that acts like a dashpot. The viscosity of the fluid in the ith-Maxwell element is adequately specified by a material coefficients denoted as $\eta_i\gt0$.

Zoom In Zoom Out Reset image size

For a generalized Maxwell rheological model with n viscous branches, the starting point for deriving the relationship between the Cauchy stress σ and the strain tensor ε is the assumption of a Helmholtz free energy potential in the form

$$\begin{equation} \begin{aligned} &\Psi\left(\boldsymbol{\varepsilon},\gamma_1,\ldots,\gamma_n\right)\, , \end{aligned} \end{equation} \tag{ 1 }$$

which is defined in terms of the total strain ε and internal variables denoting the viscous strains γ_i . The energy function can be divided into the equilibrium (elastic) part $\Psi^e$ and the non-equilibrium (viscoelastic) parts $\Psi^v$ as

$$\begin{equation} \begin{aligned} &\Psi\left(\boldsymbol{\varepsilon},\gamma_1,\ldots,\gamma_n\right) = \Psi^e\left(\boldsymbol{\varepsilon}\right) +\sum_{i = 1}^n \Psi^v_i\left(\boldsymbol{\varepsilon},\gamma_i\right)\, . \end{aligned} \end{equation} \tag{ 2 }$$

The Clausius–Planck inequality for purely mechanical problem can be rewritten for small strains as

$$\begin{equation} \mathcal{D}_{\text{int}} = \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}} - \dot{\Psi} \unicode{x2A7E} 0 \, , \end{equation} \tag{ 3 }$$

where $\dot{\boldsymbol{\varepsilon}}$ denotes the rate of the strain tensor, and $\dot{\Psi}$ is the material time derivative of the energy function. This inequality asserts that the scalar product of the stress tensor and the rate of strain tensor (stress power) must be greater than (for irreversible process) or equal to (for reversible process) the ratio of the rate of internal energy. It captures the irreversible processes and energy dissipation occurring within a material system. Calculating the time derivative of energy function equation (1) using the chain rule and substituting in equation (3) results in

$$\begin{equation} \begin{aligned} \mathcal{D}_{\text{int}} = & \left( \boldsymbol{\sigma} - \frac{\partial \Psi\left(\boldsymbol{\varepsilon},\boldsymbol{\gamma}_1,\ldots,\boldsymbol{\gamma}_n\right)}{\partial \boldsymbol{\varepsilon}} \right) : \dot{\boldsymbol{\varepsilon}} &\\ &-\sum_{i = 1}^n \frac{\partial \Psi\left(\boldsymbol{\varepsilon},\boldsymbol{\gamma}_1,\ldots,\boldsymbol{\gamma}_n\right)}{\boldsymbol{\gamma_i}}:\dot{\boldsymbol{\gamma}_i} &\unicode{x2A7E} 0 \, . \end{aligned} \end{equation} \tag{ 4 }$$

The Coleman–Noll procedure is applied to fulfill the requirement of $\mathcal{D}_{\text{int}} \unicode{x2A7E} 0$. By considering arbitrary choices of the tensor variable ε , a physical expression for the Cauchy stress σ is deduced, along with an inequality that governs the non-negativeness of the internal dissipation $\mathcal{D}_{\text{int}}$, as it is required by second law of thermodynamics. Substituting the right side of equation (2) into equation (4) results in

$$\begin{equation} \begin{aligned} \boldsymbol{\sigma} & = \boldsymbol{\sigma}^{\infty}+\sum_{i = 1}^n\boldsymbol{\sigma}^v, \quad \text{with}\\ \boldsymbol{\sigma}^{\infty}& = \frac{\partial \Psi^e\left(\boldsymbol{\varepsilon}\right)}{\partial \boldsymbol{\varepsilon}},\quad \boldsymbol{\sigma}^v = \frac{\partial \Psi^v_i\left(\boldsymbol{\varepsilon},\boldsymbol{\gamma}_i\right)}{\partial \boldsymbol{\varepsilon}} \, . \end{aligned} \end{equation} \tag{ 5 }$$

In the rheological model illustrated in figure 4, the total strain in the ith Maxwell element is the sum of the elastic spring and dashpot strains for that element. This relationship results in

$$\begin{align} \begin{aligned} \boldsymbol{\varepsilon}_i = \boldsymbol{\varepsilon}-\boldsymbol{\gamma}_i\, . \end{aligned} \end{align} \tag{ 6 }$$

In the context of viscoelasticity theory, the volumetric component does not have any impact on the viscous behavior and can be considered negligible [44, 45]. Hence, the Cauchy stress tensor σ can be expressed for incompressible material in reduced form as

$$\begin{equation} \begin{aligned} \boldsymbol{\sigma} = p\mathbf{I} +\mathbb{C}^e:\boldsymbol{\varepsilon}^{\text{iso}}+ \sum_{i = 1}^n \mathbb{C}^v_i:\left(\boldsymbol{\varepsilon}^{\text{iso}}-\gamma_i\right), \end{aligned} \end{equation} \tag{ 7 }$$

where p is the hydrostatic pressure and $\boldsymbol{\varepsilon}^{\text{iso}} = \boldsymbol{\varepsilon}-\frac{1}{3}\, \text{trace}({\boldsymbol{\varepsilon}}) \mathbf{I}$. In consideration of the observation that only the isochoric component of an incompressible material contributes to viscoelastic deformation, the material's behavior can be described using (isochoric) fourth-order stiffness tensors as

$$\begin{equation} \mathbb{C}^e = 2\mu^e\,\mathbb{I} , \qquad \mathbb{C}^v_i = 2\mu^v_i\,\mathbb{I} = \alpha_i\, \mathbb{C}^e, \end{equation} \tag{ 8 }$$

where $\mu^e\gt0$ represents the elastic shear modulus, $\mu_i^v = \alpha_i\mu^e\gt0$ corresponds to the viscoelastic shear moduli, and $\mathbb{I}$ denotes the fourth-order identity tensor. The total shear modulus, denoted as µ, can be defined as

$$\begin{equation} \mu = \mu^e+\sum^N_{i = 1} \mu^v_i\, . \end{equation} \tag{ 9 }$$

The linear equations that govern the complementary evolution of viscous strains $\boldsymbol{\gamma}^v_i$ are expressed by the following equation

$$\begin{equation} \begin{aligned} \dot{\boldsymbol{\gamma}}_i = \frac{1}{\tau_i}\left(\boldsymbol{\varepsilon} - \boldsymbol{\gamma}_i\right), \end{aligned} \end{equation} \tag{ 10 }$$

where $\tau_i\gt0$ represent the relaxation times.

In this study, n = 1 viscous branches were deemed sufficient to accurately model the viscoelastic behavior of a gMSLA 3D-printed photopolymer. The simulations were conducted using Ansys software (version 2021 R2, ANSYS Inc. Canonsburg, PA, USA). For finite element calculations, the Prony shear viscoelastic constitutive model in Ansys was selected to apply the described theory. In the ANSYS Mechanical simulation, a mesh comprising second-order tetrahedral elements was employed, totaling 220 000 elements in the domain. Ansys Parametric Design Language (APDL) was utilized for grading the structures. The corresponding code is detailed in appendix A.

All specimens considered in this study are manufactured using the commercial MSLA printer, 'Original Prusa SL1S', developed by Prusa research a.s. This printer incorporates a UV LED with a 405 nm wavelength and a monochrome LCD module for uploading masks with dimensions of $120\,\text{mm} \times 68\,\text{mm}$ and a resolution of $2560 \times 1620$ pixels. The printing process employs the UV-sensitive 'Prusament Resin Tough Prusa Orange' by Prusa Research a.s. The material composition includes epoxy resin (20%–40%), color pigment (2%–5%), and photoinitiators (3%–5%).

The geometries of the test samples are created using CAD software, and they are exported in STL file format. The masks (slices) are generated from the STL files using the 'PrusaSlicer' software with a desired layer thickness. In configuring the print parameters, the exposure time t and layer thickness h are set within the slicer software and exported with the printable file. Notably, the slicer software defaults to providing only black and white masks. Hence, for printing graded structures, 'MATLAB' (version R2022a, The MathWorks Inc.) is employed for image processing, adjusting the white pixels to achieve the desired grayscale values G. A MATLAB code for grading structures is provided in B.

The light intensity $I_e(G)$ is measured on the LCD surface using a slim photodiode power sensor (PM100D with S120VC sensor, Thorlabs Inc. Newton, NJ) with a ±5% measurement accuracy. The maximum light intensity of the printer under consideration is $I_e^\text{max} = 3.6$ mW cm⁻². The relation between grayscale and light intensity is well-approximated as $I_e(G) = G^2 \cdot I_e^\text{max}$, see also [40].

For mechanical characterization, uniaxial tensile tests are conducted following the ASTM D638 standard, employing sample type IV. A T500-1200-5kN machine (MFC Sensortechnik GmbH) is used for testing, adhering to EN ISO 7500-1, and providing a 1 µm travel resolution accuracy. To obtain the elastic and viscoelastic material responses, a minimum of three different test samples are utilized in tensile and relaxation tests. The samples have the same form as the uniaxial tension test samples in the ASTM D638 standard.

To formulate a constitutive model for the materials, four distinct grayscale values, namely G = 0.6, 0.75, 0.85, and 1.0, were selected for calibration purposes, as explained in [31]. The exposure time for the calibration process was fixed at t = 5 s, while maintaining a layer thickness of h = 0.05 mm. Note that we use a minimum value of G = 0.6 for this particular choice of exposure time, layer thickness, MSLA printer, and photopolymer material in order to ensure solidification of the resin. The viscoelastic material model was tailored to fit experimental curves obtained at three different strain rates. Initially, the elastic component of the model in terms of the shear modulus µ^e was calibrated using a uniaxial tensile test performed at a strain rate of $\dot{\varepsilon} = 2.1\cdot 10^{-5}$ 1/s. Subsequently, the viscoelastic constitutive model with n = 1 viscous branches was fitted to higher strain rates of $\dot{\varepsilon} = 2.1\cdot 10^{-3}$ 1/s and $\dot{\varepsilon} = 4.2\cdot 10^{-3}$ 1/s to determine the relaxation time τ₁ and relative modulus α₁. The resulting constitutive parameters for each grayscale value are illustrated in figure 5.

Zoom In Zoom Out Reset image size

Importantly, as can be seen in figure 5, the behavior of the constitutive parameters exhibits nearly linear characteristics across the range of grayscale values. Consequently, a linear function can be utilized to approximate the constitutive parameters $\Gamma\in\{\mu^e,\alpha_1,\tau_1\}$ for arbitrary grayscale values between G = 0.6 and G = 1.0 as

$$\begin{equation} \begin{aligned} \Gamma \left(G\right) = a\, G+b\, , \end{aligned} \end{equation} \tag{ 11 }$$

where the corresponding coefficients are provided in table 1.

Here, a linear interpolation approach is suitable for predicting mechanical properties based on grayscale variations within the specified range. This approach aligns with similar characterizations of grayscale-graded materials [31].

To avoid the need for repetitive physical experiments, simulations offer a valuable alternative for investigating diverse physical conditions and grading strategies. Thus, the primary aim now is to validate the framework using a graded 3D printed complex structure.

In this case, a shell lattice structure comprising $5\times 5\times 3$ 'Schwarz primitive' cells, with dimensions of $50\,\text{mm} \times 50\,\text{mm} \times 30\,\text {mm}$, was selected for this purpose, as illustrated in figure 6. To grade the lattice, grayscale variation was employed, ranging from G = 0.6 to G = 1.0. The gradation initiates from the bottom, with the first layer of unit cells printed using G = 1.0, followed by the second layer using G = 0.9, and concluding with G = 0.6 in the top layer. The gradation is visually represented in figure 6(c).

Zoom In Zoom Out Reset image size

To validate the grayscale-dependent viscoelastic constitutive model and finite element simulations, two test case scenarios were examined. The first scenario involves cyclic loading, where the graded lattice structure is compressed up to 2% strain and then immediately unloaded. The strain rate in loading and unloading is $\dot{\varepsilon} = 1.56\cdot 10^{-3}$ 1/s. The corresponding force-strain and force-time curves are illustrated in figure 7. The graph clearly shows a high level of agreement between the experimental results and the numerical simulations conducted using Ansys. In the second scenario, also the relaxation behavior of the material and structure are examined. In this case, the graded lattice structure is loaded to 2% compressive strain and then held at that that strain for a duration of 30 s before being unloaded. Once again, a comparison between the numerical and experimental data demonstrates a good level of agreement in figure 8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

After validating the numerical framework and simulations, we focus on exploring the influence of material grading strategies on the viscoelastic mechanical behavior of 3D lattice structures.

Like in the previous investigation, 'Schwarz primitive' shell lattices with $5\times 5\times 3$ cells such as the one shown in figure 6 are considered. Now, different grading strategies are explored, where the grayscale value varies locally, i.e. $G = G(\boldsymbol{x})$, but ensuring that the average grayscale value $\bar{G}$ in each structure is consistent as

$$\begin{equation} \bar{G} = \frac{1}{\text{vol}\left(\Omega\right)} \int_{\Omega} G\left(\boldsymbol{x}\right)\, \text{d}V \,\stackrel{!}{ = }\, 0.8\,. \end{equation} \tag{ 12 }$$

Figure 9 shows 5 different grading strategies, including (a) the initial scenario of the lattice structure presented in figure 6, which exhibits grading from G = 1.0 at the bottom to G = 0.6 at the top by assigning distinct grayscale values to individual cell layers, resulting in a discontinuous shift between layers. The subsequent case (b) involves a continuous linear grading of the lattice structure, progressing smoothly from G = 1.0 at the bottom to G = 0.6 at the top. In case (c), the structure undergoes grading from the left side with G = 1.0 to the right side with G = 0.6, where each vertical layer of cells maintains a uniform and constant grayscale. A disordered vertical cell grading strategy is demonstrated in (d) and in (e) the cells are inhomogeneously graded in a random fashion. It should be highlighted that the average grayscale value remains consistent across all grading cases with $\bar{G} = 0.8$.

Zoom In Zoom Out Reset image size

Similar to section 3.2, the 5 different shell lattice structures are subjected to cyclic compression tests with loading up to an applied strain of 2% and subsequent unloading, both at a strain rate of $\dot{\varepsilon} = 1.56\cdot 10^{-3}$ 1/s. The resulting force-strain curves corresponding to different gradings, which were all obtained numerically using finite element simulation, are depicted in figure 10. It is evident that the mechanical behavior remains largely consistent when the average grayscale values are identical, indicating minimal impact of the grading strategy on the mechanical behavior within the linear viscoelastic region.

Zoom In Zoom Out Reset image size

Finally, the rate-dependence of the energy dissipation of the graded shell lattice structure illustrated in figure 9(a) is examined.

The structure undergoes cyclic loading up to 2% strain across various strain rates, spanning from quasi-static conditions at $\dot{\varepsilon} = 1.5\cdot 10^{-5}$ 1/s to higher rates at $\dot{\varepsilon} = 0.3$ 1/s. The dots in figure 11 represent the normalized energy, calculated as the ratio of dissipated energy to absorbed energy under various strain rates. Notably, the dissipated energy is relatively low at lower strain rates. However, with an increase in strain rate, the dissipated energy attains peak values before subsequently decreasing and converging to very small magnitudes for high strain rates. Comprehending this observed behavior and the dissipation of energy across different strain rates holds significant importance in the design of structures intended for applications requiring specific or high energy dissipation.

Zoom In Zoom Out Reset image size

In this context, a rational function can be utilized to model the energy dissipation for different strain rates, adhering to the specified form

$$\begin{equation} \begin{aligned} \Pi_d\left(\dot{\varepsilon}\right) = \frac{{p_1 \cdot \dot{\varepsilon} + p_2}}{{\dot{\varepsilon}^2 + q_1 \cdot \dot{\varepsilon} + q_2}} \, , \end{aligned} \end{equation} \tag{ 13 }$$

where the coefficients are found to best approximate the values from figure 11 as

$$\begin{equation} \begin{aligned} p_1 & = 0.04848\, , \\ p_2 & = -3.843 \cdot 10^{-5}\, , \\ q_1 & = 0.07236\, , \\ q_2 & = 0.002545\, . \end{aligned} \end{equation} \tag{ 14 }$$

The solid line in figure 11 represents the fitted dissipated energy. The strain rate at which maximum energy dissipation occurs can be identified from this fit by obtaining the stationary point with $\mathrm{d}\Pi_d(\dot{\varepsilon})/\mathrm{d}\dot{\varepsilon} = 0$ at $\dot{\varepsilon} = 8.6\cdot 10^{-4}$ 1/s.

This behavior can be explained by the fact that, at lower strain rates, the viscoelastic element in figure 4 moves freely and is unable to endure stress. Consequently, the rheological system behaves like an elastic element with an elastic modulus of $E^{\infty}$. Conversely, at higher strain rates, the dashpot in the viscoelastic element behaves like a rigid element, resulting in an effective modulus of $E = E^{\infty} + \sum_i E_i$.

It is crucial to note that plastic deformation of the material is not considered in this investigation. Therefore, when designing structures, it is important to take into account that at larger deformations and higher strains, the material may undergo plastic deformation.