Self-folding Structural Design Using Multiscale Analysis on the Light-absorption Folding Behaviour of Polystyrene Sheet

Concentrated light-absorption on specific areas of polystyrene (PS) sheet induces self-folding behaviour. Such localized light-absorption control is easily realized by black-coloured line pattern printing. As the temperature in the line patterns of PS sheet increases differently due to the transparencies in each line pattern, localized thermal contraction generates folding deformation of the PS sheet. The light-activated folding technique is caused by the shape memory effect of PS sheet. The shape memory creation procedure (SMCP) is described by using molecular dynamic (MD) simulation, and the constitutive model of PS sheet is identified. This study employs the shell/cohesive line element for the folding deformation of PS sheet, and utilizes the constitutive model obtained from the MD simulation. Based on the continuum-model analysis of the PS sheet folding deformation activated by light, various self-folding structures are designed and manufactured.


Multiscale continuum model analysis of self-folding structures
The self-folding structure using the PS sheet is activated by light including infrared ray, and its folding behaviour can be predicted from the continuum model analysis mentioned in the next chapter. The folding behaviours of some self-folding structure models are estimated by the continuum model analysis. The constitutive deformation models of polystyrene obtained from MD simulations are applied to the continuum model analysis. Figure 1 shows the simple hexahedral shape using PS sheet that is first performed. The self-folding structure with  26 . The folding angle of the hexahedron is 90°, so the line pattern width is determined to be 2.1 mm. The rectangular cell is set to the cell size 2 × 2 cm, and the PS sheet has the thickness 0.23 mm. Because each line pattern simultaneously starts folding deformation, the transparencies of the line patterns are the same as 0%. The surrounding temperature is fixed as 40 °C, and the folding deformation starts at the glass transition temperature of polystyrene sheet of 102 °C. Figure 2 shows the actual model and deformation process of the second model, which is a self-folding structure with twisted shape. The line patterns are arranged in parallel, and the PS sheet with rectangle shape is cut in the diagonal direction. As for the hexahedral shape model, the line pattern width is set to 2.1 mm, due to its folding angle of 90°, and the model is designed to the size 2 × 8 cm.
The last model is of hexahedral shape using sequential folding behaviour. As aforementioned, the self-folding structure with complex shape requires precise folding sequences between the respective line patterns. Accordingly, although the self-folding structure with hexahedral shape such as Fig. 1 shows successful simulation result, the actual model can sometimes fail to reach the planned shape. For stable self-folding success of the structure with complex shape, a sequential folding technique such as Fig. 3 is considered 30 . The final shape in Fig. 3(c) is the same as the model in Fig. 1(a), but its stability is higher than the simple hexahedral shape, due to the collar faces attached to the edge line of the PS sheet, and the sequential folding design. Because the faces of the PS sheet have to be folded at their folding sequences, the respective transparencies are differently applied to each of the line patterns, as in Fig. 3(b). In the line patterns, the time to approach glass transition temperature is determined by the transparency, so the faces can be sequentially folded.

Conclusion
The self-folding behaviour of a PS sheet activated by light irradiation is designed by fabrication of the line pattern and corresponding localized light absorption. Especially, in order to maximize the applicability of non-contact actuation, sequential folding structure is manufactured by modifying transparency of the pattern. Because the self-folding functionality is induced by shape recovery of microscopic polymer chain, MD-based multiscale approach is developed for precise prediction and design. The shape-memory cycle for polystyrene sheet is performed, and then the behaviour of the pre-strained polystyrene sheet is estimated at the nanoscale. The influence of pre-strain on thermal in-plane shrinkage is identified, and it is implemented to the continuum model. Also, elastic and hardening moduli along the pre-stretched direction are estimated as a function of internal and external parameters. These are used to provide the simplified linear elasto-plastic deformation model. In continuum scale folding model, the cohesive line element is employed to describe the discontinued rotation angle at the folded line. The final deformed configuration obtained via multiscale framework shows a good agreement with that of actual experiment. We expect our results can be used to design various complex shape and precisely control sequential deformation of the self-folding structure. Figure 4 shows a schematic diagram of the multiscale framework which connects the molecular shape recovery motion and continuum scale folding deformation of the PS. First, the shape memory cycle is described by the MD simulation. The dependence of thermal shrinkage on initially applied strain is quantified by monitoring change in the shape of unit cell during thermal transition. The microscopic information is transferred to continuum scale analysis which implements the finite shell and cohesive line elements to predict the folding angle. Then, uniaxial compressive loading is applied to each deformed MD unit cell to get the stress-strain curve. We propose simple equations to estimate the elastic and hardening moduli as a function of the temperature and initially defined shape. These are utilized to construct the simplified bilinear elasto-plastic constitutive deformation model. This section introduces the computational methodology of both molecular dynamics and continuum simulation to provide the prediction model of mechanical response of the self-folding PS.

Methods
Multiscale constitutive model construction using molecular dynamic simulation. This section describes the atomistic modelling and MD simulation methodology to describe the shape memory creation cycle to impose the shape memory effect to the PS polymer chains, and implement the microscopic design parameters into the continuum-scale model. Our previous study 42 revealed that only pre-stretched PS shows self-folding behaviour at the glass transition temperature. This is because the main mechanism of the thermal contraction is  the shape recovery from anisotropic to isotropic state. The previous results also showed that the thermoelastic properties and final recovered configuration of the self-folding PS are determined by the shape-programming path. In this study, the linear elasto-plastic deformation model of the oriented PS sheet is expressed as a function of two key parameters, which are the pre-imposed strain, and the temperature. Furthermore, the relation between the pre-strain applied to the PS, and its local shrinkage, is quantified, using the MD calculations and experimental values. This relation is implemented to the continuum scale prediction of the macroscopic folding angle and final deformed configuration.
Atomistic modelling of the initial bulk PS unit cell and relaxation simulation are performed by using Materials Studio 5.5 (Accelys, Inc.), which uses the Teodorou-Suter method to construct the amorphous structure 44 . The intramolecular and intermolecular interactions are described by the polymer consistent force field (PCFF) 45 . This can identify the thermomechanical properties and glass transition behaviour of the polymeric materials, because the potential coefficients are derived from the experimental properties of the organic compounds. The atom-based summation with a cutoff distance of 9.5 Å is used to calculate the non-bonded interaction between atoms. In treating the Coulombic interaction, the distance-dependent dielectric method with a dielectric value of 2.6 46 is used.
The MD unit cell consists of atactic PS polymer chains, and the periodic boundary condition is applied to every surface of the unit cell for modelling the bulk system. In order to ensure the physical crosslinking between chains, the molecular weight of the PS should be larger than its entanglement molecular weight (M e = 19,000 g/ mol). The entanglement is required to realize the shape recovery deformation, because it plays the role of a net point for fixing the original shape. Therefore, two polymer chains with a molecular weight of 20,831 g/mol each (400 repeating styrene units) are used to construct the unit cell. Energetically favourable structure is made by applying a sequential equilibration procedure: energy minimization using the conjugated-gradient method, NVT ensemble for 500 ps at 300 K, and the NPT ensemble for 2.5 ns at 300 K and 0.1 MPa. The time step increment is 1.0 fs. In order to ensure that the polymer chains are distributed without any anisotropy, the annealing simulation, including five cycles of the NVT ensemble for 250 ps, is performed. The temperature reduces from 900 to 725 K at each cycle, and the scalar orientation order parameter of the PS backbone chains (this will be explained below, see also Eq. 1) converges to zero after the annealing process. Figure 5 shows the configuration of the amorphous bulk unit cell.
Prior to applying the shape memory creation treatment, our model is utilized to calculate the thermoelastic properties to confirm sufficient equilibration and validate the reliability of simulation methods. The glass transition temperature and linear coefficient of thermal expansion (CTE) of relaxed unit cell are determined by monitoring the variations of specific volume and mean square displacement (MSD) of the system within the temperature range of 300-500 K. The relaxed structure at 500 K and 0.1 MPa is cooled to 300 K with a cooling rate of 10 K/ns under isobaric condition. Then, the specific volume-temperature plot is obtained by averaging the specific volume for the last 100 ps of the NPT equilibration steps (1 ns) at each temperature. The candidate of T g within the simulated temperature range is approximated as the temperature at which the MSD of the system abruptly increases. It is because the glassy-rubbery transition induces huge increase in the chain mobility. The approximated point is used to divide the glassy and rubbery states of the unit cell, and T g is defined as the intersection of the linear regression lines at each region. The linear CTE can be obtained by using the slope of specific volume-temperature plot. The procedure to characterize the glass transition and calculate the thermoelastic properties is described in our previous study in detail 42 . These properties showed good agreement with experimental values.
The thermomechanical cycle is applied to the initial isotropic polymer chains, in order to produce the shape memory effect. All MD simulations describing the shape memory cycle are carried out using the LAMMPS code developed by Sandia National University 47 . The relaxed unit cell is heated up to T h = 550 K, above the thermal transition temperature (T g = 383 K and T m = 513 K), by applying an equilibration step, which consists of energy minimization, NVT, and NPT ensemble simulations for 1.0 ns, respectively. After the heating up process, the thermomechanical test begins by biaxial stretching the PS chains along the in-plane directions (x and y) with a nominal strain rate of 10 /s 8 . The deformed states with six different programmed strains (ε m = (15, 30, 45, 60, 80, and 100) %) are prepared to estimate the relation between the initial molecular anisotropy and corresponding mechanical deformation behaviour, and the extent of thermal shrinkage. Then, each MD cell is cooled down to T l = 300 K using the NVT ensemble for 1 ns, and the external loading is removed by using the NPT ensemble for 1 ns, respectively. During the stress relaxation, the nominal strain slightly decreases to the fixed strain (ε u ). Final in-plane thermal shrinkage of the oriented PS is realized by re-heating it with a rate of 5 K/ns under isobaric condition. The unit cell is heated up to the temperature at which the molecular chains recover to the isotropic state ((550-650) K). During contraction to the recovered shape (ε n ), the cell lengths along all three dimensions are controlled independently to reflect the anisotropic movement of the chains. The time step for the entire thermomechanical cycle is 0.5 fs. Figure 6 shows the entire scheme of the thermomechanical cycle conducted in the presented MD simulation study.
Entropic shape recovery is caused by microscopic conformational transformation of the polymer chains. This microscopic motion and corresponding anisotropic/isotropic transition during shape memory cycle can be effectively captured using orientational order parameter (P i ) 42 , which can be calculated via: Here, θ i is the angle between the i-axis of the Cartesian coordinate, and the vector of the polymer backbone segments. The polymer chains are perfectly isotropic when P i equals to zero, and they become in aligned state as P i approaches 1.0. As the PS unit cell is stretched up to ε m = 100%, the order parameter along the loading direction increases gradually up to 0.1. All the unit cells with different ε m are heated and Fig. 7 shows corresponding in-plane thermal contraction strain of the oriented PS with respect to temperature. The lengths of unit cells are monitored, until the polymer structure is fully recovered to the isotropic state (P i = 0). As the pre-applied strain increases, the shape recovery ratio ( decreases. This is because the unrecoverable strain and disintegration of the polymer network occur with larger deformation. However, Fig. 7 shows that the total shrinkage increases as the unit cell is more deformed. Table 1 shows microscopic order parameter and shape memory properties of the oriented PS with different pre-strain. Previously, we have established a relation between the initial programmed strain and total in-plane shrinkage by using MD data points, and this can predict the experimental results with larger ε m 48 . This relation is used to predict final folding angle of the localized heated region of the oriented PS sheet.  The operating temperature for shape recovery deformation estimated in the MD simulation is slightly higher than the experimental T g ((382-388) K). The temperature shifting is due to the limitation of short time-scale, which cannot observe all of the conformational recoveries at the experimentally suggested temperature. However, our simulation schemes can clearly observe the shape changes induced by thermal transition of the polymer chain molecules. Furthermore, in contrast to the instantaneous heating technique 49,50 , the long-time gradual heating used in this paper has the advantage of characterizing the dependence of the microstructural parameters on the macroscopic self-folding behaviour of the PS.
The macroscopic self-folding angle and final deformed configuration of the PS sheet are considerably varied by the local temperature distribution and intrinsic mechanical anisotropy of the irradiated region. Therefore, the compressive mechanical response of the pre-stretched PS is parametrically studied in terms of the operating temperature and pre-strain using the MD simulation. Each pre-deformed MD unit cell (ε m = (0, 30, 60, 80) %) after stress relaxation at 300 K is compressed along the in-plane directions with a constant strain rate of 10 8 /s under the NPT ensemble. The applied strain ranges up to 15%, and it is known that both elastic and plastic deformations of the PS can be clearly seen within this strain range 51 . The virial formula is used to calculate the stress at each deformed state, and the operating temperature ranges from (300 to 440) K. The results from three MD runs with different initial atomic velocity distribution are averaged to present the stress-strain relationship of the PS. Figure 8 shows the change of stress-strain curve with temperature and programmed strain. Figure 8(a) shows that the mechanical response of the PS to the uniaxial compressive loading can be assumed as linear elasto-plastic behaviour. The linear elastic regime lasts within a strain range of ε = (0-3) %, and yielding occurs at ε = 4%. To simplify the model, the dependence of the yield point on the temperature and the morphology of the oriented PS are not considered. As the temperature rises, the slope of curve at the elastic regime (elastic modulus (E)) is almost linearly decreased. In the case of the non-oriented system, E decreases with a rate of 10.3 MPa/K at low temperature below T g . This decreasing rate is reasonable, compared with other literature 51 , which is of the same order as the experimental value (5 MPa/K). With heating up above the glass transition temperature, the modulus drops to about 1.38 GPa at the viscous rubbery regime. The decreasing tendency of the modulus in this measurement can also be clearly seen in the plastic deformation. The slope of the s-s curve at the plastic regime (K) is about 0.6 GPa at 300 K, which is 77% smaller than the elastic modulus. As the temperature rises within the (300-380) K range, K decreases with a rate of 3.9 MPa/K, and drops to 0.25 GPa at the rubbery state.
The influence of mechanical anisotropy on the deformation behaviour, dominantly determined by the extent of the pre-applied strain, is also quantified. Figure 9(a) and (b) show variations of tangential stiffness at the elastic and plastic regions at 300 K, respectively. The mechanical properties (E and K) are increased with increasing the pre-strain. The biaxial in-plane stretching applied during the manufacturing process of the oriented PS imposes molecular orientation on the macromolecular chains. Therefore, the increase of mechanical anisotropy with greater extent of pre-stretching is evident, even within the fully plastic regime, as well as the linear elastic regime. E and K are fitted by an exponential rise to the optimal model, because they converge to specific value. The solid lines in Fig. 9(a) and (b) are the fitted curves, and the square of the correlation coefficients is about 0.81.
In order to construct the constitutive model with general expression, two slopes of the s-s curve at different regimes are expressed as a function of the temperature and pre-strain. In order to make a simple model, the influences of these two input variables are independently considered, and hence it is assumed that the decreasing rates of modulus with respect to the temperature are all the same, regardless of the pre-strain. The final prediction equation for E and K at temperature below T g can be expressed as: For precise prediction of the self-folding deformation, it is important to predict how much in-plane contraction will occur at the stage of glassy-rubbery transition of the polymer. As discussed above, the extent of shape recovery has close relation to the degree of anisotropy of the microstructure. Our previous work 42 provides a model to expect total in-plane shrinkage of the oriented PS as a function of the programmed strain with practical range (ε m = (150-300) %). For example, the polystyrene sheet used in our study shows shrinkage of about 50% compared with the original shape, and its pre-strain can be estimated to be 200% using this model. Because the folding behaviour of the oriented PS occurs without any external mechanical loading, the stress-strain curve is shifted by the amount of the predicted contraction at = T T g . The deformation at temperature above T g is regarded as plastic deformation, because the self-folding of PS is not recoverable, unless another thermomechanical cycle is applied to it. Combining the MD simulation-based prediction of mechanical parameters and in-plane shrinkage together, the constitutive deformation model with various pre-strain and temperature is constructed just as shown in Fig. 10: Shell finite element. This section employs the shell element used in Brank et al.'s research 40,41 . The shell director is defined using an Euler rotation matrix, and it is assumed as the vector with magnitude 1, which is defined as a unit sphere. The shell director is expressed by the combination of two Euler angles, which are the rotation angles α and β. The rotation angles α and β are defined from the x-axis and y-axis rotation, respectively.  In the undeformed configuration, it is assumed that the x-y axis rotation angles are initially zero. Then, the relationship between the undeformed shell director → T and the deformed shell director → t follows as: The rotation matrix Q is determined by the rotation angles α and β, and is separated into two parts, α Q ( ) 1 and β Q ( ) 2 . Each matrix has rotation angles on its respective rotation axis (x,y). Also, if the undeformed shell director → T is assumed to 0 0 1, which is orthogonal to the surface, the deformed shell director → t can be expressed as: First, the undeformed/deformed position vectors are expressed as: and mean the position vector, mid-surface position vector, and shell director defined in the undeformed configuration, respectively; → → → r x t , and denote the position, mid-surface position vector, and shell director in the deformed configuration, respectively; and θ is the curvilinear coordinates. Greek indices are 1 and 2, while Latin indices run from 1 to 3. The displacement vector defined from the mid-surface position vectors has the following relationship: The displacement vector → w is induced from the position vectors as: where, the vector  → is expressed as → = → g i i   in the curvilinear coordinate system, and  i is the contravariant component.
The covariant basis vectors are defined at the material point as: where,  → i means the covariant base vector. In addition, the metric tensors in the undeformed/deformed configuration are obtained. Using the metric tensor, the Green-Lagrangian strain tensor separated into the membrane, bending and transverse shear strains, is obtained as: ,  ,  ,  ,  ,  ,   ,  ,  ,  ,  ,  ,  ,  ,  ,  , , , Total energy can be expressed by the Green-Lagrange strain and the second Piola-Kirchhoff stress tensor S. Then, by the approximated shell theory, the internal energy If the high-order term of θ 3 is neglected in the → g and → a formulation, the assumption = + g a / 1 ≈1 is considered. In Eq. 11), as the isotropic and linear elastic material condition is applied, the resultant force, moment, and transverse shear force can be respectively expressed as: where, k is the transverse shear correction factor that is determined to be 5/6. By applying the virtual work principle, the variational term of the internal energy can be obtained as: A The strain measure variational terms are obtained from the directional derivative formulation of virtual displacement and virtual rotations. Then, they are calculated as: , , The folding deformation of the PS sheet shows large displacement and rotation angle, so it requires non-linear deformation analysis. The iterative Newton-Raphson method is employed for the non-linear problem in this paper, and the tangential stiffness matrix and the internal force vector should be constructed from the variational form of the internal energy. The differentiation for the kinematic variables x and t is performed, and it follows that: where, the matrix Λ is defined as: sin( )sin( ) cos( )cos( ) cos( ) 0 sin( )cos( ) cos( )sin( ) (16) The tangential stiffness matrix is obtained from the differentiation of the variational internal energy form. The tangential stiffness matrix is usually separated into the geometric part and material part, as: The individual derivative terms of the membrane, bending, and transverse shear strain measures in the geometric part are expressed as: , , When the finite element discretization is applied to the geometric tangential stiffness matrix, the terms are formulated as:  where, B ijkl is the 4th order constitutive tensor expressed in Eq. 12). B ijkl is determined by the constitutive relationship between the strain and stress, so it is changed via the stress that the material is given. However, when the polystyrene sheet is folded by light-absorption, the major deformation occurs at the black-coloured line pattern. The non-printed area of the polystyrene sheet, where the shell element is applied, does not generate large strain. Therefore, the linear elastic model is assumed in the non-printed area of the polystyrene sheet.
Cohesive line element. The folding deformation has been expressed by using the bending deformation considering the cohesive line element, damage model, crack and hinge line, due to its localized deformation and discontinuous rotation 38,39,[52][53][54][55] . In the folding deformation analysis, this section employs the cohesive line element, among the various analysis methods for folding deformation. The cohesive line element provides additional nodes at the folded line, and then the discontinued rotation angles generated at the folded line can be handled without difficulty. This section introduces the cohesive line element used in Giampieri et al. 's study 38,39 .
In the cohesive line element, the additional nodes are defined at the previous existing nodes. Figure 11 shows that accordingly, the nodes have the same coordinates.
The coordinate system → e i in the cohesive line element is the co-rotation frame defined in the reference interface surface. In addition, the kinematic variables at the line are interpolated by the linear shape function ζ N( ); the position vector at the line is then expressed as: Also, the covariant base vector in the undeformed configuration can be defined as: where, → G n is obtained from the other two orthogonal covariant vectors. The cohesive line element allows the displacement jump for the discontinuous rotation angles. The displacement jump at the interface line is expressed as: 3 Also, the deformed position vector and shell director at the line are defined as: 3 From Eq. 24, the covariant basis vector for the deformed interface surface is expressed as: The total internal energy is separated into the shell element part and the cohesive line element part, and it is shown as: where, → P and → w are the internal stress and strain, respectively, at the cohesive line. The variational form of internal energy can be easily formulated. The variational term of the displacement jump is expressed as: 3 where, the virtual shell director jump is: In Eq. 28, ϑ → denotes the rotation axis vector, which is expressed as: where, Ω denotes the skew-symmetric tensor.
are induced from Eqs 28 and 29: 3 3 Then, the variational displacement in the interface line is derived Substituting Eq. 31 into Eq. 26, the variation form of the internal energy δU Line is separated into the in-plane part and out-of-plane part, as:  where the drilling moment M t is assumed to be zero. As the shell element, the linearization process of internal energy is performed in the cohesive line element for the non-linear deformation problem. The tangential stiffness is directly obtained from the derivative of the virtual internal energy, and it is separated into geometric part and material part, as: g Line A L L and each term of the material tangential stiffness matrix is formulated as: