Cold-programmed shape-morphing structures based on grayscale digital light processing 4D printing

Shape-morphing structures that can reconfigure their shape to adapt to diverse tasks are highly desirable for intelligent machines in many interdisciplinary fields. Shape memory polymers are one of the most widely used stimuli-responsive materials, especially in 3D/4D printing, for fabricating shape-morphing systems. They typically go through a hot-programming step to obtain the shape-morphing capability, which possesses limited freedom of reconfigurability. Cold-programming, which directly deforms the structure into a temporary shape without increasing the temperature, is simple and more versatile but has stringent requirements on material properties. Here, we introduce grayscale digital light processing (g-DLP) based 3D printing as a simple and effective platform for fabricating shape-morphing structures with cold-programming capabilities. With the multimaterial-like printing capability of g-DLP, we develop heterogeneous hinge modules that can be cold-programmed by simply stretching at room temperature. Different configurations can be encoded during 3D printing with the variable distribution and direction of the modular-designed hinges. The hinge module allows controllable independent morphing enabled by cold programming. By leveraging the multimaterial-like printing capability, multi-shape morphing structures are presented. The g-DLP printing with cold-programming morphing strategy demonstrates enormous potential in the design and fabrication of shape-morphing structures.


Supplementary Note 2: Finite element modeling 2.1 Constitutive model
The rubbery material B3 is modeled as an incompressible neo-Hookean solid as it is elastic in the studied temperature range (25-100 o C).The cold-draw and shape memory behaviors of the glassy material (B1 and B2) are captured using an incompressible, stress-and temperature-dependent, multi-branch viscoelastic model, which consists of one equilibrium branch (neo-Hookean model assumed) and multiple nonequilibrium branches (modeled by visco-hyperelastic Maxwell elements) in parallel.In this model, the Cauchy stress σ is where F is the deformation gradient, b = FF T is the left Cauchy-Green tensor, the subscript "t" or "s" denotes the deformation from time 0 to t or s.
where Eeq is the modulus of the equilibrium branch, Ei the instantaneous modulus of branch i, Eins the total instantaneous modulus; they are taken to be small-stress moduli and assumed to be temperature-independent.τi R denotes a reference relaxation time of the branch i.
The shift factor, α(t), written as a function of time t in Eq.( 1) for brevity, is taken to not only depend on temperature but also stress at time t.Therefore, α can be written as [2] T S α α α = where α T denotes the temperature-dependent shift factor, α S the stress-dependent shift factor.
The temperature-dependent shift factor α T is taken to follow the Williams-Landel-Ferry (WLF) approximation for T>Tref, i.e., ) and the Arrhenius form for T<Tref, i.e., ( ) In Eqs.( 4 The stress-dependent shift factor α S is taken to follow [2] ( ) where M is the equivalent shear stress, M is the Mandel stress, M'=dev(M), ES is the activation energy, and s is the athermal shear strength.To capture the softening effects, s where s0 is the initial value of s, ssat is the saturation value of s, h0 is a prefactor, γ v is the viscous strain, and v γ is the viscous flow rate.For a specific branch, ( ) with μi denoting the instantaneous shear modulus and R i i τ ατ = the relaxation time.Thus, τi R can be interpreted as the relaxation time under small stress and at the reference temperature Tref.Here, to capture essential physics while facilitating the numerical implementation, we use to define an effective viscous flow rate, instead of that of a specific branch.μins=Eins/3 due to the incompressibility.τ0 is treated as a fitting parameter.

Numerical implementation
We implement the material model into Abaqus through the built-in viscoelastic model in which the Prony-series and user-defined shift factor can be specified.The Prony-series parameters are specified as those under small stress and at the reference temperature Tref (at which α=1).The shift factor Eq.( 3) is implemented through user subroutines USDFLD and UTRS.USDFLD is invoked to access stress values of an integration point and update α S following Eqs.( 6)-( 9), which will be passed into UTRS to update α T and finally α.The updates are in an explicit (or forward Euler) manner.
We also implement the material model into a MATLAB script for the uniaxial tension to assist the parameter identification.We derive the following integration algorithm for the stress updates.The stress at an arbitrary time t, i.e., Eq.( 1), can be rewritten as Then, the stress at t+δt can be updated as where The trapezoidal integration rule has been used to derive Eq.(12).

Parameter identification
The SMP parameters including Eeq, Ei, τi R , C1, C2, and AFc/kb are identified using DMA results from multi-temperature and multi-frequency tests.[3] Temperature ramps from 10℃ to 130℃ in 5℃ increment.At each temperature, the material is kept isothermal and stress-free for 5 min and then subjected to a cyclic loading at multiple frequencies (0.1Hz~20Hz).The SMP parameters are identified in two steps.First, the storage modulus E' versus test frequency ω data at different temperatures can be shifted to a master curve at Tref, obtaining discrete data points of α(T).Performing nonlinear numerical fitting to these data points using Eqs.( 4)-( 5) yields the material constants C1, C2, and AFc/kb (see Supplementary Figure 8a for B1 and 9a for B2).Second, with the multi-branch model, the storage modulus E' and loss modulus E'' can be written as '' 1 Using Eqs.( 13) and (14), performing nonlinear numerical fitting to the shifted E' and E'' versus ω data yields the material constants Eeq, Ei, τi R (see Supplementary Figure 8b for B1 and 9b for B2).The obtained parameters are listed in Supplementary Table 1 and 2.
The cold-draw parameters s0, ssat, h0, τ0, and ES/kb are identified using uniaxial tensile data for different strain rates at room temperature.The parameters are roughly estimated following the procedure of [2] and further optimized by numerical fitting.The obtained parameters lead to simulated stress-strain curves in good agreements with the experiments at multiple strain rates for the two glassy materials (Supplementary Figure 10).The obtained parameters are listed in Supplementary Table 1 and 2.
All parameters identified above are further used to simulate a complete process of the cold-drawing (i.e., stretching, holding, releasing) and heat-recovering.The simulated stress and strain curves, along with the temperature history, achieve good agreements with the experiments for the two glassy materials (Supplementary Figure 11).denote the temporarily fixed strain of the glassy fiber due to stretch, which is the mismatch strain between two material phases.Then the longitudinal strain (perpendicular to cross-section) for an arbitrary point of position y can be written as

Supplementary
)-(5), C1, C2, and AFc/kb are material constants, T is the current temperature, Tref represents a structural transition temperature.T and Tref are in Kelvin scale.Note that α T = 1 at T =Tref.

where
εNP is the longitudinal strain of a neutral plane (NP), representing the in-plane stretching due to balance of forces; κy' represents the bending strain due to balance of moments, with κ denoting the curvature of NP (positive value means curling up) and y'=y-yNP the relative position to NP (position yNP).It is seen that NP has zero longitudinal strain due to bending.The balance laws for forces and moments are 0 Ef and Em are the Young's modulus of the fiber and matrix, respectively.Inserting (15) into (16) gives the curvature κ

Table 1 .
Material parameters for the B1 material.