Computer Simulations of Continuous 3-D Printing - Macromolecules

Article pubs.acs.org/Macromolecules

Computer Simulations of Continuous 3‑D Printing Zilu Wang, Heyi Liang, and Andrey V. Dobrynin* Department of Polymer Science, University of Akron, Akron, Ohio 44325, United States S Supporting Information *

ABSTRACT: 3-D printing is a revolutionary manufacturing technique which makes it possible to fabricate objects of any shape and size that are hard to reproduce by traditional methods. We develop a coarse-grained molecular dynamics simulation approach to model the continuous liquid interface production (CLIP) 3-D printing technique. This technique utilizes a continuous polymerization and cross-linking of the liquid monomeric precursor by the UV light within a thin layer while pulling the cross-linked polymeric object out of a pool of monomers. Simulations show that the quality of the shape of the 3-D printed objects is determined by a fine interplay between elastic, capillary, and friction forces. Using simulation results, we identify the source of the object shape deformations and develop a set of rules for calibration of the parameters to meet the accuracy requirements. Comparison between different continuous 3-D printing setups shows that proposed modifications of the printing process could improve quality and accuracy of the printed parts.

■

within a region of a finite thickness. To explore the effect of the friction and capillary forces, we performed simulations of the modified CLIP process, as illustrated in Figure 1, that we named CLIP-I and CLIP-Ir printing techniques. In these simulations, we reverse design of the printing process in such a way that printed parts remain submerged in the precursor resin. Keeping parts submerged in the precursor resins reduces the effect of capillary forces acting on a printed part as it is pulled from the pool of monomers in the CLIP setup. Thus, our simulations show that with decreasing the printed feature size a competition between elastic, capillary, and friction forces is the main factor responsible for the process quality.

INTRODUCTION 3-D printing technology has revolutionized the parts-on-demand approach and found applications in tissue engineering,1−8 structural materials design,9−12 artisan and do-it-yourself 3-D printing, fabrication of microfluidic devices,13,14 and soft robots and actuators.15−17 These applications take advantage of the variety of additive manufacturing techniques such as fused deposition, stereolithography, laser sintering, and melting.18 A broad application of 3-D printing in manufacturing would require a new level of control over the quality and speed of the printing process. This could place 3-D printing technology beyond the research laboratories and prototype development niche to become economically viable alternative to conventional manufacturing. Recent development of the stereolithography 3-D printing technique such as the continuous liquid interface production (CLIP)19,20 addresses some of the outlined limitations by achieving new benchmarks in the speed of the printing process without sacrificing its accuracy. The CLIP method takes advantage of the oxygen inhibition21−23 of the free radical polymerization in photopolymerizing UV-curable resins, which maintains a layer of liquid at the bottom of the resin container during the polymerization via oxygen-permeable membrane, and this liquid layer, or the so-called “dead zone”, is crucial for a continuous supply of the precursor resins to a curing region. While there is a substantial interest in this printing technology, the fundamental understanding of the physical mechanisms controlling quality of the printing process is lagging behind. In order to elucidate factors influencing accuracy of the continuous 3-D printing techniques based on the photopolymerization, we have developed a coarse-grained molecular dynamics simulation approach to 3-D printing by modeling photopolymerization as a polymerization process occurring © XXXX American Chemical Society

■

RESULTS AND DISCUSSION In our coarse-grained molecular dynamics simulations24 of the 3D printing process we use a bead and bead−spring representation of the coarse grained monomers (uncured or precursor resin), cross-linked resin, printing platform, and surrounding walls (see Figure 1). The interactions between beads with diameter σ are modeled by the truncated-shifted Lennard-Jones (LJ) potential. The values of the Lennard-Jones interaction parameter εLJ between different species are given in Table S1 (see Supporting Information). All bonds connecting monomers are modeled by the bond potential represented by a sum of the FENE and LJ potential.25 During polymerization process (resin curing) each active monomer can form up to three bonds with its neighbors. In our simulations, to mimic photopolymerization of the UVcurable resins, monomer polymerization reaction can only occur Received: August 8, 2017 Revised: September 11, 2017

A

DOI: 10.1021/acs.macromol.7b01719 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 1. Snapshots of the simulation boxes showing system setups for molecular dynamics simulations of different techniques. (a) In the CLIP technique the printed part is attached to the platform and emerges from the pool of monomers during printing process. Two alternative CLIP-like printing techniques: CLIP-I (b) and CLIP-Ir (c). In these setups printed parts remain submerged in the pool of monomers. In the CLIP-Ir technique the printed part is attached to the bottom of the monomer pool and remains stationary during printing process. Monomers are shown in gray color, cured resin is shown in light gray color, and beads belonging to a platform are shown in yellow color.

Figure 2. Time evolution of the part growth during the printing process using three different setups.

within a window of thickness 2.0σ as shown in Figure 1. A monomer could form a bond with an active monomer belonging to a cured resin (network) if it is located within a reaction radius Rreact. A new bond could be added to an already connected bead as long as the number of bonds per bead is less than the maximum possible number and a bead remains within a polymerization window. By varying probability of chain propagation and cross-linking reactions, we can control modulus of the cured (cross-linked) resin. In our simulations we have implemented three different setups for 3-D printing: CLIP, CLIP-I, and CLIP-Ir as illustrated in Figure 1. In CLIP simulations (see Figure 1a) we reproduced the printing setup proposed by Carbon3D with polymerization of the resin occurring in a pool of monomers with a part being attached to a platform which is pulling it up from the monomer pool.19 The polymerization window is located between the bottom of the pool and the surface of the monomeric liquid. This is illustrated in Figure 2 that shows evolution of the part growth during printing process. In CLIP-I simulation setup (see Figures 1b and 2) the printing (resin curing) occurs within a monomer pool with printed part being removed from a pool after 3-D printing process is complete. CLIP-Ir simulation setup is similar to CLIP-I method (see Figures 1c and 2). The printed object is attached to the bottom of the simulation box and is growing up as polymerization window moves up with a top platform as shown in Figure 2. The dimensions of the polymerization zone and “dead zone” are the same in all our simulations. In our simulations the printing velocity v is defined as a velocity with

which polymerization window or platform to which the part is attached is displaced during printing process. If not specified, the printing velocity v is equal to v = 0.01σ/τLJ. Figure 3 shows final results of the simulations of the 3-D printing of parts by three different techniques. The first column shows shape of the original part. It follows from this figure that the best quality of part reproduction is achieved by the CLIP-Ir technique. This method minimizes the effect of the friction and capillary forces acting on the moving part as it is displaced to free a polymerization window for a new “layer” of monomers. The CLIP technique appears to be the second best and has comparable quality with CLIP-Ir for solid or bulkier objects such as cylinder, tube, mesh, and vessel. However, the quality of the printed spring and chain objects for the CLIP method demonstrates a large degree of distortion in comparison with similar objects printed by the CLIP-Ir method. The worst performance in terms of part’s quality is observed for the CLIP-I method. This is due the large friction force exerted by surrounding monomer liquid on the moving printed part. In the case of the CLIP technique the effect of the friction force is reduced since it is only acting on the section of the printed part which is still located inside the monomer pool. However, in the case of the CLIP method there are additional capillary forces acting on the part at the surface of the monomer pool. These forces should be considered as an additional source for reduction of the part’s quality. The CLIP-Ir technique is free from such defect generating forces since the printed part does not move and remains inside a monomer liquid pool during the entire printing B

DOI: 10.1021/acs.macromol.7b01719 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 3. Final shapes of the 3-D printed replicas obtained from molecular dynamics simulations. Target 3-D shape templates (first column), parts obtained by using CLIP (second column), CLIP-I (third column), and CLIP-Ir (fourth column) printing techniques printed with velocity v = 0.01σ/τLJ.

process. Thus, implementation of CLIP-Ir technique could offer an advantage for improvement of the quality of the printed parts. It is worth pointing out that utilization of this technique could only be beneficial for small size objects that could fit inside a monomer pool. For the large size objects the CLIP technique offers an advantage for printing parts of relatively high quality and at the same time avoiding using the large volumes of monomers. To highlight effect of the friction forces between printed part and monomeric liquid, we study effect of the printing velocity v on the final part’s quality. This is illustrated in Figure 4, which shows variation of the outer radius of the cylindrically shaped mesh along z-axis with the printing speed in the CLIP-I and CLIP-Ir printing techniques. There is a significant decrease in the radius in the middle with increasing the printing speed for parts produced by CLIP-I technique. These radius variations across the part height are due friction forces acting on the part surface causing shape deformations. This effect diminishes in the case of the CLIP-Ir technique. However, we still observe a small distortion of the radius with increasing the printing speed. This is due to monomer influx patterns generated close to the polymerization window. Geometric object dimensions are convenient to describe quality of printed parts of simple geometrical shapes. For parts of complex shapes to characterize accuracy of the part printing process we introduce the fidelity score

Figure 4. Dependence of the shape and radius of the object on the printing speed. Template has an outer radius RT = 35σ.

Q (z ) =

N (z ) −1 NT(z)

(1)

It is defined as the ratio of the number of pixels N(z) in a printed part at height z within a slice of thickness 2.0σ and a part cross C

DOI: 10.1021/acs.macromol.7b01719 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 5. Fidelity score for 3-D printed parts. (a) Illustration of procedure for calculation of the fidelity score by digitizing the beads in slice of printed object into pixelated image. (b) Variation of the fidelity scores along the z-axis for different templates.

Figure 6. Deformation of the cylindrically shaped objects during printing as a function of the cross-linking density.

section area to that number in the template (see Figure 5a). The distribution of parameter Q(z) along the z-axis is shown in Figure 5b. The perfect replication of the part corresponds to fidelity factor Q(z) = 0. Any deviation from zero points out on the part distortion. The larger the absolute value of the fidelity factor is the larger deviation of the printed part is from its template (see Figure 3). Note that the negative value of Q(z) (see Figure 5b) indicates shrinkage of the cross-section area of the printed part in comparison with that of the template.

In order to quantify the effect of the friction and capillary forces and establish how the effect of these forces could be optimized by increasing the cross-linking density in the cured resin, we have performed a set of simulations of printing a cylindrically shaped object by CLIP, CLIP-I, and CLIP-Ir techniques. Figure 6 shows dependence of the cylinder shape distortion on the probability of the cross-linking reaction vs linear chain growth during photopolymerization process. It follows from this figure that at high cross-linking density which D

DOI: 10.1021/acs.macromol.7b01719 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 7. Schematic representation of forces acting on the slice of the deformed cylinder during the printing process.

Taking a limit Δz → 0, we obtain the following differential equation connecting rate of change of the stress, radius R(z), and the velocity v(z)

corresponds to a larger modulus of the cured resin there is not much difference in the final shapes of the printed objects. However, for cross-linking densities smaller than 0.1σ−3, roughness start to appear on the surface of the cylinders. When the cross-linking density decreases further, there is a significant shape variation close to the platform for all printing techniques. This shape distortion is due to a fine interplay between elastic, capillary, and friction forces close to the printing platform. Note that this effect is more pronounced for the CLIP-I printing technique. Thus, a particular flow pattern at the beginning of the printing process could have a strong effect on the quality of the reproduction process. As printing shifts further away from the platform, the quality improves and cylindrical shape is recovered. This shape “healing” is driven by a balance between friction and elastic forces acting on the cylinder. Note that variations in the cross-linking density shown in Figure 6 could result in a 1000-fold decrease in the system modulus (see Supporting Information for details). To quantify contribution of different forces acting on the printed part as it increases in size, we consider forces acting on the slice of the resin with thickness Δz and radii R(z) and R(z +Δz) (see Figure 7). We set the zero of the z-axis at the surface of the moving platform. There are stresses σzz(z+Δz) and σzz(z) acting on two surfaces of a slice and capillary forces of net magnitude 2πγR(z+Δz) and 2πγR(z) acting along the circumference of the two boundaries of the slice, where γ is the surface tension of the monomer/resin interface or air/resin interface when the sample is pulled out of the monomer pool such as in the CLIP technique. In addition to internal stresses and capillary forces there is also a friction force acting on the side surface of the slice. This force is proportional to the relative velocity v(z) of the solvent with respect to that of the moving part with which the printing platform is pulled, F⃗(z) = −v(⃗ z)ξ2πR(z)Δz, where ξ is the friction coefficient. Here we will assume that the derivative |dR(z)/dz| ≪ 1 such that in calculating projections of the capillary and friction forces on the z-axis we can substitute cos θ(z) ≈ 1 in the leading term approximation. In the steady state, the sum of projections on the z-axis of all forces acting on the slice of the cured resin is equal to zero:

d(R(z)2 σzz(z)) d R (z ) + 2γ = −2ξv(z)R(z) dz dz

(3)

For a cured resin the stress σzz(z) of the uniaxially deformed slice in the limit of small deformations is equal to σzz(z) ≈ 3Gε(z)

(4)

where ε(z) is the strain of deformation of the cylinder along the z-axis and G is the shear modulus of the cured resin. For incompressible networks, the deformation along the z-direction is related to deformation in the radial direction as λR (z) = 1/ 1 + ε(z) ≈ 1 − ε(z)/2 . Therefore, taking into account the expression for network stress (eq 4) and assuming radius of the undeformed cylinder to be R0, we can rewrite eq 3 as follows: ⎛ 2ξv(z) γ ⎞ d ε( z ) ≈− ⎜1 − ⎟ 3GR 0 ⎠ dz 3R 0G ⎝

(5)

It follows from this equation that there are two different regimes for which solution of strain ε(z) as a function of z is positive or negative. The crossover between these two regimes is determined by the elastocapillary number γ/GR0.26−28 For small elastocapillary numbers, γ/GR0 ≪ 1, solution of eq 5 gives ε(z) < 0, which corresponds to increase of the radius of the cylinder R(z) ≈ R0 (1 + z/δG) with characteristic length scale δG ≈ ξv0/3R0G (where we substituted instead of v(z) a typical value of the velocity v0 at length scale δG). However, for large elastocapillary numbers, γ/GR0 ≫ 1, the solution of eq 5 gives ε(z) > 0. This promotes elongation of the cylinder along z-axis and shrinkage in the radial direction R(z) ≈ R0(1 − z/δγ) with a characteristic length scale δγ ≈ ξv0/γ being inversely proportional to the surface tension. The presented above analysis of the deformation of the cylindrically shaped objects qualitatively explains cylinder shape transformations seen in Figure 6. At low cross-linking densities close to the foundation of the cylinder the structure of the crosslinked resin resembles that of a weakly cross-linked brush layer with shear modulus approaching zero. In this case we can neglect the elastic term in eq 3 and consider only balance of the capillary and friction forces. Solution of the eq 5 gives ε(z) > 0 resulting in shrinkage of the cylinder in the radial direction. This regime

π (R(z + Δz)2 σzz(z + Δz) − R(z)2 σzz(z)) + 2πγ(R(z + Δz) − R(z)) + 2πR(z)Δzξν(z) ≈ 0 (2) E

DOI: 10.1021/acs.macromol.7b01719 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

(4) Murphy, S. V.; Atala, A. 3D bioprinting of tissues and organs. Nat. Biotechnol. 2014, 32 (8), 773−785. (5) Derby, B. Printing and Prototyping of Tissues and Scaffolds. Science 2012, 338 (6109), 921−926. (6) Atala, A.; Kasper, F. K.; Mikos, A. G. Engineering Complex Tissues. Sci. Transl. Med. 2012, 4 (160), 160rv12−160rv12. (7) Gross, B. C.; Erkal, J. L.; Lockwood, S. Y.; Chen, C.; Spence, D. M. Evaluation of 3D Printing and Its Potential Impact on Biotechnology and the Chemical Sciences. Anal. Chem. 2014, 86 (7), 3240−3253. (8) Xing, J.-F.; Zheng, M.-L.; Duan, X.-M. Two-photon polymerization microfabrication of hydrogels: an advanced 3D printing technology for tissue engineering and drug delivery. Chem. Soc. Rev. 2015, 44 (15), 5031−5039. (9) Coulais, C.; Teomy, E.; de Reus, K.; Shokef, Y.; van Hecke, M. Combinatorial design of textured mechanical metamaterials. Nature 2016, 535 (7613), 529−532. (10) Sydney Gladman, A.; Matsumoto, E. A.; Nuzzo, R. G.; Mahadevan, L.; Lewis, J. A. Biomimetic 4D printing. Nat. Mater. 2016, 15 (4), 413−418. (11) Kokkinis, D.; Schaffner, M.; Studart, A. R. Multimaterial magnetically assisted 3D printing of composite materials. Nat. Commun. 2015, 6, 8643. (12) Qin, Z.; Compton, B. G.; Lewis, J. A.; Buehler, M. J. Structural optimization of 3D-printed synthetic spider webs for high strength. Nat. Commun. 2015, 6, 7038. (13) Au, A. K.; Huynh, W.; Horowitz, L. F.; Folch, A. 3D-Printed Microfluidics. Angew. Chem., Int. Ed. 2016, 55 (12), 3862−3881. (14) Waheed, S.; Cabot, J. M.; Macdonald, N. P.; Lewis, T.; Guijt, R. M.; Paull, B.; Breadmore, M. C. 3D printed microfluidic devices: enablers and barriers. Lab Chip 2016, 16 (11), 1993−2013. (15) Wehner, M.; Truby, R. L.; Fitzgerald, D. J.; Mosadegh, B.; Whitesides, G. M.; Lewis, J. A.; Wood, R. J. An integrated design and fabrication strategy for entirely soft, autonomous robots. Nature 2016, 536 (7617), 451−455. (16) Rus, D.; Tolley, M. T. Design, fabrication and control of soft robots. Nature 2015, 521 (7553), 467−475. (17) Bartlett, N. W.; Tolley, M. T.; Overvelde, J. T. B.; Weaver, J. C.; Mosadegh, B.; Bertoldi, K.; Whitesides, G. M.; Wood, R. J. A 3Dprinted, functionally graded soft robot powered by combustion. Science 2015, 349 (6244), 161−165. (18) Gibson, I.; Rosen, D.; Stucker, B. Additive Manufacturing Technologies: 3D Printing, Rapid Prototyping, and Direct Digital Manufacturing; Springer: 2014. (19) Tumbleston, J. R.; Shirvanyants, D.; Ermoshkin, N.; Janusziewicz, R.; Johnson, A. R.; Kelly, D.; Chen, K.; Pinschmidt, R.; Rolland, J. P.; Ermoshkin, A.; Samulski, E. T.; DeSimone, J. M. Continuous liquid interface production of 3D objects. Science 2015, 347 (6228), 1349− 1352. (20) Janusziewicz, R.; Tumbleston, J. R.; Quintanilla, A. L.; Mecham, S. J.; DeSimone, J. M. Layerless fabrication with continuous liquid interface production. Proc. Natl. Acad. Sci. U. S. A. 2016, 113 (42), 11703−11708. (21) Ligon, S. C.; Husár, B.; Wutzel, H.; Holman, R.; Liska, R. Strategies to Reduce Oxygen Inhibition in Photoinduced Polymerization. Chem. Rev. 2014, 114 (1), 557−589. (22) Dendukuri, D.; Panda, P.; Haghgooie, R.; Kim, J. M.; Hatton, T. A.; Doyle, P. S. Modeling of Oxygen-Inhibited Free Radical Photopolymerization in a PDMS Microfluidic Device. Macromolecules 2008, 41 (22), 8547−8556. (23) Yagci, Y.; Jockusch, S.; Turro, N. J. Photoinitiated Polymerization: Advances, Challenges, and Opportunities. Macromolecules 2010, 43 (15), 6245−6260. (24) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: 2001; Vol. 1. (25) Kremer, K.; Grest, G. S. Dynamics of entangled linear polymer melts: A molecular-dynamics simulation. J. Chem. Phys. 1990, 92 (8), 5057−5086. (26) Cao, Z.; Dobrynin, A. V. Polymeric Droplets on Soft Surfaces: From Neumann’s Triangle to Young’s Law. Macromolecules 2015, 48 (2), 443−451.

continues until the modulus of the cross-linked cylinder reaches its bulk value changing a sign of the left-hand side of eq 5, γ/GR0 ≪ 1. After this takes place the cylinder radius begins to increase until it reaches the boundary of the polymerization window. The presented above analysis of the shape deformation during printing process indicates that we can introduce two dimensionless parameters γ/GR0 and ξv0/GR0 to control printing quality. These parameters should also be used for mapping results of the coarse-grained simulations of the 3-D printing process into real systems.

■

CONCLUSIONS We have developed a coarse-grained molecular simulation approach to model 3-D printing techniques based on photopolymerization of monomers. By comparing three different printing setups, we have established factors controlling fidelity of the 3-D printing process. Our simulations have shown a strong effect of the modulus of the cured resin (cross-linking density) on the quality of the printed cylindrically shaped objects. The observed shape deformation pattern could be explained by a scaling model that takes into account balance of elastic, capillary, and friction forces acting on a cylinder during 3-D printing process. This analysis highlights limitations of the described 3-D printing technique for printing of soft objects. The coarsegrained molecular dynamics simulation approach developed herein creates an opportunity for optimization of the 3-D printing setup and selection of the materials for efficient shape reproduction of 3-D printed parts.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.7b01719. Simulation details: polymerization rules, polymerization kinetics, and materials characterization (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (A.V.D.). ORCID

Zilu Wang: 0000-0002-5957-8064 Heyi Liang: 0000-0002-8308-3547 Andrey V. Dobrynin: 0000-0002-6484-7409 Funding

This work was supported by National Science Foundation (DMR-1624569). Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Hollister, S. J. Porous scaffold design for tissue engineering. Nat. Mater. 2005, 4 (7), 518−524. (2) Jakus, A. E.; Rutz, A. L.; Jordan, S. W.; Kannan, A.; Mitchell, S. M.; Yun, C.; Koube, K. D.; Yoo, S. C.; Whiteley, H. E.; Richter, C.-P.; Galiano, R. D.; Hsu, W. K.; Stock, S. R.; Hsu, E. L.; Shah, R. N. Hyperelastic “bone”: A highly versatile, growth factor−free, osteoregenerative, scalable, and surgically friendly biomaterial. Sci. Transl. Med. 2016, 8 (358), 358ra127−358ra127. (3) Kang, H.-W.; Lee, S. J.; Ko, I. K.; Kengla, C.; Yoo, J. J.; Atala, A. A 3D bioprinting system to produce human-scale tissue constructs with structural integrity. Nat. Biotechnol. 2016, 34 (3), 312−319. F

DOI: 10.1021/acs.macromol.7b01719 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules (27) Cao, Z.; Stevens, M. J.; Dobrynin, A. V. Elastocapillarity: Adhesion and Wetting in Soft Polymeric Systems. Macromolecules 2014, 47 (18), 6515−6521. (28) Carrillo, J.-M. Y.; Dobrynin, A. V. Molecular Dynamics Simulations of Nanoimprinting Lithography. Langmuir 2009, 25 (22), 13244−13249.

G

DOI: 10.1021/acs.macromol.7b01719 Macromolecules XXXX, XXX, XXX−XXX