《1. Introduction》

1. Introduction

Four-dimensional (4D) printing [1–7] is a manufacturing process that utilizes stimuli-responsive (smart) materials [8–10], mathematics, and additive manufacturing (AM). These three elements yield an encoded static smart structure that becomes a dynamic smart structure under the right stimulus through an interaction mechanism. 4D printing technology combines materials science and mathematics [11]. In 4D printing, in addition to smart materials and AM, mathematical models are further needed for accurately predicting shape-morphing behaviors over time [6,12]. In simpler cases, these models may not be necessary.

The structures obtained by 4D printing technology are 4D space–time dependent [6,13–15]. The time-dependent behavior of 4D-printed structures relies on the special spatial arrangement between active and passive materials in the three-dimensional (3D) space (encoded by printing). Therefore, the fourth dimension in 4D printing indicates predictable and desired evolution over time that depends on 3D space, and this 3D space is the specific arrangement of (physical) voxels determined by mathematics and realized by printing. Prior to 4D printing, researchers did not typically try to generate a specific printing path (mainly by an inverse mathematical problem) to obtain a predictable and desired evolution over time. 4D-printed products and structures have unique advantages, some of which are discussed here:

《1.1. Why 4D printing?》

1.1. Why 4D printing?

Using 4D printing, programming is performed at the material level [1,3]. Sensing and actuation are directly embedded into materials [3]. Random (free) energy is used to construct non-random structures [1]. Self-assembly is performed at the macro scale to eliminate post-fabrication assembly [3]. Buildings can be selfassembled on Mars [16]. The number of parts in a product is reduced [3]. Volume in shipping is reduced using 4D structures that are activated on-site [2]. The number of failure-prone components used in conventional robots and electromechanical systems is minimized [3]. 4D printing can be used for embedding bioinspired morphing features into materials [15]. More advantages can be mentioned for 4D printing in various fields.

《1.2. Why 4D by printing?》

1.2. Why 4D by printing?

First, the same reasons that motivate us to use AM instead of other manufacturing processes for conventional materials motivate us to use AM for smart materials. Some of these reasons are as follows: Molds and dies are eliminated, material waste is reduced, and complex geometries are obtained. Second, using printing, we can elaborate on (i.e., manipulate) the "internal structures” of multimaterials. In other words, 4D "printing” enables local anisotropy to be encoded [5]. This is uniquely important in 4D printing for obtaining various shape-morphing behaviors. Encoding the structures of materials is not achievable (or hardly achieved) by other manufacturing processes. As one example, bioengineers are primarily motivated to produce tissues and organs through bioprinting because a proper function of an organ depends on its internal structure, in which various cells are contained in a complex architecture and multi-scale vasculature [17]. Systematic yet local elaboration on the internal structure is not easily achieved by other fabrication processes. For instance, bioengineers do not use bioinjection molding or other methods to mimic biological structures that have highly complicated internal architectures.

Following this Technology, Entertainment, Design (TED) talk in 2013 that introduced 4D printing [1], Tibbits [2] published a paper on 4D printing entitled "4D printing: multi-material shape change.” Truby and Lewis [18], in their references section, indicated that the aforementioned paper by Tibbits [2] "describes the first embodiment of 4D printing.” Note the word "multimaterial” in the title of the work done by Tibbits [2]. This word has been used to introduce 4D printing. Although some studies have demonstrated single-material 4D-printed structures, the future of 4D printing resides in multi-material structures [2,19– 29]. Most current researchers consider 4D printing in a multimaterial fashion using smart (active) materials that are selectively arranged with conventional (passive) materials to obtain desirable shape-shifting behavior [22]. Utilizing a single (active) material for the entire 4D structure has been less considered [22]. The physics of 4D-printed structures usually requires the use of multimaterials [20]. A well-balanced 4D-printed component is intrinsically a multi-material structure [23]. In the most fundamental case, a multi-material 4D-printed structure has one active and one passive layer [21]. Performance improvement by allocating proper materials to related locations based on local necessities [30], multi-functionality by embeddable functions such as electronics [31], combining rigid and flexible sections in an integrated structure [32], and providing lightweight structures [33] are only some of the advantages of multi-material structures. Moreover, to enable the shape memory effect (SME) (which is not an intrinsic property of shape memory polymers (SMPs) [10]), a mechanical force is required in addition to heat (thermomechanical cycles). This force is usually provided externally by human intervention in single-material SMPs. However, 4D printing can help us arrange active and passive materials in a multi-material structure and utilize their internal mismatch-driven forces to enable the SME autonomously without requiring an external load (i.e., human intervention) for training [34,35]. Furthermore, combining SMPs with a passive layer (multi-material structure) can be one of the simplest means of producing reversible (two-way) SMPs [36–41]. For these reasons, multi-material structures have received growing attention in 4D printing. This means that when researchers perform general research on 4D printing, they typically elaborate on the active–passive combination (multi-materials) [21–24,42–46]. In 4D printing, it is best to arrange active and passive materials to obtain self-sensing and self-actuating operations, as no conventional sensors or actuators are used to enable predictable and complex shape-shifting behaviors. Self-sensing can be achieved by using active (smart) materials, and self-actuating can be achieved by arranging active and passive materials to generate some mismatch-driven forces and guide these self-actuation forces in the desired direction (as derived from mathematics). Nevertheless, the results of this study can be useful for both single-material (active) and multi-material (active–passive) structures, as discussed in a remark later in this article.

Various aspects of 4D printing have been explored in the literature. Several studies have considered beam and plate theories (e.g., see the supplementary information provided in Ref. [5]). However, no study has considered a general modeling of "time-dependent” behaviors (the fourth dimension) of 4D structures. Time-dependent behaviors are critical aspects of 4D (stimuli-responsive) materials, whether such materials are fabricated by printing (and thus called 4D-printed structures) or produced by other manufacturing processes. More importantly, numerous studies on 4D materials have used the Timoshenko bimetal model [47], which is linear with time, to analyze the time-dependent behaviors of their experiments. In this study, we observe that the Timoshenko bimetal model cannot capture the true time-dependent behaviors of 4D materials (except in some special cases or selected linear regions), although it does provide useful insights into time-independent behaviors, as many studies have correctly used the Timoshenko bimetal model for analyzing time-independent behaviors such as the effect of thickness on bending. In fact, the purpose of the Timoshenko bimetal model is not to model the time-dependent behaviors (the fourth dimension) of 4D materials. Therefore, an urgent need exists for qualitative and quantitative analyses of the fourth dimension of 4D materials. Some of the different formulas used or developed for modeling the shape-morphing behaviors of 4D-printed structures are organized and analyzed in Table 1 [5,20,21,37,48–60].

《Table 1》

Table 1 Different formulas used or developed in 4D printing for modeling shape-morphing behaviors.

The symbol definitions and notations in the table are relatively independent.

Now, we provide a discussion about Table 1:

(1) The formulas used or developed in 4D printing for the modeling of the shape-morphing behaviors can be divided into five categories: ① pure geometry, ② customized and case-specific models, ③ Euler–Bernoulli model, ④ Timoshenko bimetal model, and ⑤ extension of each of these models to plates (e.g., spherical and saddle-like structures), for which two curvatures are used (usually the mean and Gaussian curvatures, which are functions of the two principal curvatures). This is not a strict classification, but is instead for general understanding.

(2) Some studies have modeled the shape-morphing behaviors of 4D materials based purely on geometry. These models are basic and usually must directly measure one of the key parameters related to the shape, such as the curvature radius or bending angle. Therefore, these models cannot be used as standalone models for predicting shape-morphing behaviors.

(3) Some studies have developed customized models for their specific cases (e.g., [21,51]). These models are not typically used for general purposes and may not consider all of the important parameters.

(4) Several studies have directly used the Euler–Bernoulli model or customized it for their specific cases.

(5) Most studies have directly used the Timoshenko model (more specifically, the Timoshenko bimetal model [47]) or customized it for their specific cases. In addition, some studies have used Timoshenko (and Euler–Bernoulli) model without mentioning the name or work of Timoshenko (and Euler–Bernoulli).

(6) Finally, each of these categories can be extended to plates. One example that extended the Timoshenko bimetal model to 2D has been previously mentioned [5] (and more references could be cited for each of the aforementioned categories, particularly for the Timoshenko and Euler–Bernoulli models).

(7) The models that are based on pure geometry are a subset of all other categories. For example, Refs. [21,51] also used geometrical relationships in their derivations, in which they assumed that the energy generated in the contraction of the active layer would lead to the deformation of the passive layer. Those geometrical relationships are also intrinsically considered in the Euler– Bernoulli and Timoshenko models.

(8) Timoshenko’s theory is an improved and advanced version of the Euler–Bernoulli theory. The Timoshenko bimetal model is an even more advanced and relevant version.

(9) Even though most of the studies in 4D materials used the Timoshenko bimetal model, those based on pure geometry or the Euler–Bernoulli model can be regarded as a subset of the Timoshenko bimetal model. In addition, the customized (case-specific) models are not general and may also ignore some important parameters.

(10) Three main points can be summarized as follows: ① No general formula exists for properly modeling the time-dependent behaviors and that is suitable for various stimuli and materials; ② the Timoshenko bimetal model [47] seems to be the most trusted formula in 4D materials (mainly for modeling timeindependent behaviors); and ③ in our study, we consider the Timoshenko bimetal bimetal model [47] mainly for the initial step (i.e., the equilibrium (balances of forces and moments) and compatibility conditions). The Timoshenko bimetal model cannot capture the time-dependent behavior of various 4D materials. We introduce and embed "time” into the Timoshenko bimetal model systematically and in detail.

The main aspect of 4D-printed structures is the fourth dimension. However, no general formula currently exists for modeling and predicting this additional dimension. Here, by developing fundamental concepts and starting from the equilibrium and compatibility conditions, we derive a necessary bi-exponential formula for modeling and predicting the fourth dimension of any multi-material 4D-printed structure. We further validate our bi-exponential formula using various experimental data from previous independent studies and show that it is a general formula useful for various 4D (e.g., photochemical-, photothermal-, solvent-, moisture-, heat-, and ultrasound-responsive) multimaterial structures. This generality is observed (when analyzing various experimental data) because we build the bases of our biexponential formula comprehensively.

《2. Definitions, derivations, and discussions》

2. Definitions, derivations, and discussions

The time-dependent behavior (i.e., fourth dimension) of any 4D-printed structure must be predicted. A detailed but systematic study on 4D printing and related areas enables us to identify three general laws that govern the shape-morphing behaviors of nearly all 4D structures, even though many materials and stimuli exist (Fig. 1). The purposes of these laws are first, to understand, and second, to model and predict the fourth dimension.

《Fig. 1》

Fig. 1. Toward 4D printing laws. RH: relative humidity.

《2.1. First law of 4D printing》

2.1. First law of 4D printing

Nearly all shape-morphing behaviors (e.g., photochemical-, photothermal-, solvent-, pH-, moisture-, electrochemical-, electrothermal-, ultrasound-, enzyme-, hydro-, and thermoresponsivity) of multi-material 4D-printed structures originate from one fundamental phenomenon, which is "relative expansion” between active and passive materials.

This "relative expansion” is the origin of nearly all complicated 4D printing shape-morphing behaviors such as twisting, coiling, and curling, which are enabled by encoding various types of anisotropy between active and passive materials and fabricating different heterogeneous structures.

《2.2. Second law of 4D printing》

2.2. Second law of 4D printing

The shape-morphing behaviors of nearly all multi-material 4D-printed structures have four types of physics: mass diffusion, thermal expansion, molecular transformation, and organic growth, which are discussed, quantified, and unified in the subsections that follow.

All of them lead to the relative expansion between active and passive materials and consequent shape-morphing behaviors under stimuli. The stimulus is usually provided externally, but it can be internal.

2.2.1. Quantifying the second law

Here, we describe and quantify the four underlying physical concepts that lead to relative expansions between active and passive materials in multi-material 4D structures, resulting in various shape-morphing behaviors, with or without an SME.

2.2.1.1. Mass diffusion. In this category, a matter transport leads to relative expansion. Mass change due to absorption or adsorption of a guest matter (here, this is mainly a stimulus such as water or an ion) in a host matter can be modeled as follows [61–63]:

where t is time and m is mass. In addition, q and depend on swelling behaviors (host matter relaxation and guest matter diffusion) and are curve-fitting constant and time constant, respectively. Instead of curve fitting, models can be developed for q and . The exponential model given in Eq. (1) captures the correct behavior of mass diffusion for short- and long-time processes. One other model also exists called a power function, which is not accurate for long-time processes of mass diffusion [61]. However, we know that the main part of 4D printing is the intelligent behavior over time, which can be short or long.

The exponential model given in Eq. (1) has been mainly used for and (volumetric strain). However, we have

where V is volume, L is length, is strain, and x, y, and z are the Cartesian coordinates. In addition, the second- and third-order differential quantities are neglected. For an isotropic material,

In addition, mass and volume have a linear relationship. Therefore,

where C1 and are constants that depend on the previous parameters.

We quantify all four categories of the second law in terms of strain. One reason for this can be given by the following example: If we want to quantify the categories in terms of mass, then we do not have any mass change in the next category (thermal expansion). Similarly, we do not have temperature change in this category, whereas we do have it in the next category.

In this category, several stimuli and mechanisms can be used, including hydro-, solvent-, moisture-, pH-, enzyme-, photochemical-, and electrochemical-responsive mechanisms. All of these stimuli ultimately cause a matter transport, leading to relative expansion in multi-material 4D structures and resulting in various shape-morphing behaviors. It should be noted that these stimuli are some examples that can cause a matter transport. However, some of them are not unique to this category. The type of material and the interaction between the material and stimulus are important. For example, in the solid state, a photochemicalresponsive mechanism can lead to a shape-morphing behavior without any mass change. In this example, the photochemicalresponsive mechanism causes a molecular transformation (the third category).

2.2.1.2. Thermal expansion. In this category, a temperature change will increase (or decrease) the average distances between atoms and molecules (with constant mass), leading to relative expansion in multi-materials.

Strain due to temperature change is [64]

where α is the coefficient of thermal expansion (CTE) and T is temperature. Because we must predict the behavior of 4D-printed structures over time, we convert this temperature-based equation to a time-based form. By applying a thermal stimulus with temperature T2 to a structure with temperature T1 and assuming a uniform temperature in the structure, the temperature of the structure changes over time as follows [65]:

where T1( t = 0 ) is the initial temperature of the structure, is energy conversion in the structure, Q1 is the heat transferred between the structure and environment (other than the stimulus), and R1–2 is the (effective) thermal resistance between the thermal stimulus with temperature T2 and the structure with temperature T1. In addition, is a time constant that depends on density (ρ), heat capacity (cp), volume (V), and thermal (and thermal contact) resistance (R1–2), and can be modeled in a specific application [65]. It mainly takes the general form of [65]. Elaborating on Eq. (6), we obtain

By combining Eqs. (5) and (7), we obtain

where C2 and are constants that depend on the previous parameters. In this category, several stimuli and mechanisms can be used, including photothermal-, electrothermal-, and ultrasound-responsive mechanisms. All of these stimuli ultimately raise the temperature of the structure and consequently increase the average distances between atoms and molecules. For example, in electrothermal-responsive structures, the movement of a current through a resistance provides heat (so-called Joule or Ohmic heating). This heat increases the temperature, ultimately leading to expansion. Similarly, contraction is obtained by cooling. It should be noted that some materials will shrink (contract) upon heating. In these cases, α would be negative in the previous equations, and the form of the final equation (i.e., Eq. (8)) will remain intact.

2.2.1.3. Molecular transformation. In this category, a molecular transformation with constant mass and temperature leads to relative expansion. The molecular transformation in the structure can be triggered by an electric field, magnetic field, light, mechanical force, etc. For example, by applying an electric or magnetic field to appropriate materials, dipoles in the material will be aligned along the field, leading to a shape change; or by applying an ultraviolet (UV) light to a photostrictive material, as a result of trans-to-cis conversions, a shape change will occur. Similarly, by applying a mechanical tension to a polymer, the chains will be aligned along the force direction, resulting in a shape change. Regardless of the type of stimuli, these phenomena have some common features. First, these molecular transformations usually start from a randomly or differently oriented condition. Second, they require time to be fully realized. Third, they have a saturation limit (e.g., when all dipoles align, all trans-to-cis conversions occur, or all chains align). Based on the original definition of the derivative,

where  is a function of (i.e., ), and thus

The strain  and thus  are generally small numbers (this is a logical assumption, otherwise the stress–strain relationship becomes nonlinear). We approximate the function  around =0 using the Taylor series as follows:

The first term of the expansion is zero because , . In addition, we can ignore and so on ( is generally a small number). Moreover, this type of response (molecular transformation) has a saturation limit, as previously mentioned. Therefore, when is around zero (this usually occurs, as is small), we can approximate it as , where is the saturation limit strain. By considering these three points,

where ;   is a time constant that depends on the rate of the molecular transformation in the material, the strength of the stimulus, and other factors. Now, by combining Eqs. (10) and (12), we obtain

Eq. (13) is derived using a relatively similar approach as in Ref. [66]. However, in their derivation [66], the format of the function  was unknown and the authors indicated that the first term of the expansion, , was zero. According to the authors [66], this was because when time tends to large numbers, the strain approaches the saturation limit. Their justification is not valid because when tends to zero, an unknown function does not necessarily tend to zero. It can in fact tend even to a very large number. Therefore, we devised a different approach. In addition, Ref. [66] was only about the magnetic field. The solution of the differential Eq. (13) would be

This behavior is observed, for example, when: ① An electric field moves the charged regions in the related materials, leading to a strain [67]; ② a magnetic field moves the magnetized regions in the relevant materials, resulting in a strain [66]; ③ a light source causes photochemical conversions (that are different from a photothermal mechanism) in the appropriate materials, leading to a strain [68,69]; or ④ an external pressure moves polymer chains, thus producing a strain [70]. In all of these mechanisms, a molecular transformation with constant mass and temperature leads to a shape-morphing behavior. Therefore,

where C3 and are constants that should be modeled in a specific application based on the type and characteristics of the stimulus as well as material properties. For example, in a photochemicalresponsive mechanism, in crystals, C3 and depend on the irradiation flux intensity, quantum yield of the transformation, and other factors, and can be modeled in a specific case [69]. We introduced the constant-temperature feature in this category to indicate that the relative expansion of this category is not due to thermal expansion/contraction, even though some thermal fluctuation may occur in the structure due to bond cleavage or formation. In this category, various stimuli and mechanisms can be used as previously discussed.

2.2.1.4. Organic growth. In this category, a living layer (organism) exists and its growth over time can lead to relative expansion between active and passive materials. The organic growth can be defined as the increase in the weight or length of an organic system [71]. This usually occurs in areas of bioscience and bioengineering dealing with, for example, cells, soft tissues, organs, and scaffolds, which can be 4D printed; these are generally referred to as 4D bioprinting.

Kinetics of the organic growth is [71]

where L(t) is the length of the organic system, is the final length, and L0 is the initial length. In addition, is usually a curve-fitting time constant that depends on the metabolism of the living organism, environment, etc. Nevertheless, models can be developed for it. This formula shows the growth of individual organisms in a population, which is different from population growth. The population growth is the growth in the number of individuals in a population and is modeled by other formulas [72,73]. Next, based on the definition of strain, , we have

where C4 and are constants that depend on the previous parameters. In this category, one of the main stimuli that can trigger a living organism is the electrical signal (electrochemical mechanism). In addition, various stimuli such as pH, light, heat, and enzymes can be used to tune the growth rate.

2.2.2. Unified model of the second law

By quantifying the second law, we found that

where Cj and (j = 1, 2, 3, 4) are all constants. However, they depend on different factors as previously described.

《2.3. Third law of 4D printing》

2.3. Third law of 4D printing

Time-dependent shape-morphing behavior of nearly all multi-material 4D-printed structures is governed by two "types” of time constants. For the most fundamental case of a multi-material 4D-printed structure having one active and one passive layer (Fig. 2), the time-dependent behavior (in terms of "curvature,” which is a building block concept for shape-morphing behavior) is

where is the curvature induced under stimulus; r is the radius of curvature; h, , and are thicknesses identified in Fig. 2; E1 and E2 are Young’s moduli; and I1 and I2 are the second moments of area. The passive and active layers are denoted by numbers 1 and 2, respectively. HI is a constant that depends on Young’s moduli of the active and passive layers and the amount of the mismatch-driven stress generated at the interface. is a time constant that depends on the viscosity induced at the interface and Young’s moduli of the active and passive layers. The viscosity induced at the interface must be measured or modeled for a specific active–passive composite. The exact format of HI and can be developed in a specific application depending on the active–passive composite and by using parallel and series rules of springs (elasticity elements) and dashpots (viscosity elements). HII and are respectively equivalent to Cj and (j = 1, 2, 3, 4) of the unified model distilled in the previous section. Finally,  and  that are introduced as model coefficients for simplicity.

《Fig. 2》

Fig. 2. Toward the third law by analyzing the most fundamental multi-material 4D structure. (a) Before applying a stimulus; (b) after applying a stimulus; (c) cross-sectional view of an element of the bilayer after applying a stimulus. The passive and active layers are denoted by numbers 1 and 2, respectively. P1 : force in the passive layer; P2 : force in the active layer; M1 : moment in the passive layer; M2 : moment in the active layer; Ai : cross-sectional area of the active or passive layer; zi : distance from the neutral axis of the active or passive layer; σi : stress in the active or passive layer.

2.3.1. Proof of the third law

To derive Eq. (20), we start from the equilibrium and compatibility conditions that are the starting point for any related problem in mechanicsofmaterials[64].WealsoconsidertheTimoshenkobimetal model [47] and its basic assumptions. We consider the Timoshenko bimetal model mainly for the initial step of the equilibrium(balancesofforcesandmoments)andcompatibilityconditions.

Equilibrium:

First, we must have balances of forces and moments in Fig. 2. Therefore [47],

where P1 and P2 are forces; M1 and M2 are moments (as shown in Fig. 2).

Compatibility:

Second, at the interface of the two layers, the lengths of the two layers are the same after the stimulus has been applied. Because their initial lengths are also the same, their strains must be equal and thus [47,64],

The strain in each of the two layers has three main contributors as follows [47,74]:

Strain from curvature ( εcurvature ) [47]:

Strain from stress ( εstress ) [47]:

Next, in the following, we develop εstress and εexpansion for 4D multi-materials and incorporate them in the equilibrium and compatibility equations.

In this study and for 4D multi-material structures, we note that εstress is the strain due to the mismatch-driven stress at the interface of the active and passive materials. The mismatch-driven stress naturally leads to opposing resistive forces in the two layers (in general). We include this resistive effect by expanding the wellknown moment equation for each layer Mi . By considering Fig. 2, we obtain

where Mi  is moment in the active or passive layer, σi is stress in the active or passive layer, Ai is cross-sectional area of the active or passive layer, and zi is distance from the neutral axis of the active or passive layer as indicated in Fig. 2. Because the integral is the same as the summation (i.e., it is a continuous summation), we could separate the integral into two terms given as Eq. (28). The second term on the right-hand side of Eq. (28) shows the integral over an infinitesimal cross-sectional area dAi (Fig. 2) that is close to the interface. Therefore, Eq. (28) can be written as:

where Mm is moment due to the mismatch-driven stress. The first terms on the right-hand side of Eq. (29) are similar to those proposed by Timoshenko [47], and the second terms ( ) are introduced in this study. At this stage, are moments. Let us retain them as black-box terms.

Strain from expansion (εexpansion):

When nearly all types of shape-morphing mechanisms in multimaterial 4D structures are analyzed, the relative expansions induced under stimuli can be categorized into four main groups elaborated in the second law. We previously (in the second law) demonstrated that nearly all types of strains due to expansions induced by stimuli have the same format as follows:

We use as the time constant of εexpansion to distinguish it from , which will be introduced for the strain resulted from the mismatch-driven stress (εstress).

By combining Eqs. (21), (22), and (27), we obtain

where N =  is introduced for simplicity. We now assume that the active (stimuli-responsive) material is responsive under stimulus, whereas the passive material is not responsive under stimulus, as their names imply; nevertheless, the relative expansion is important for the shape-morphing behavior. Therefore, by applying Eq. (30) to Eq. (31), we obtain

Now, each term in Eq. (32) is strain. Therefore, and are strains. In addition, these two terms reflect the mismatch (viscosity) effect of the interface into each layer. The viscoelastic strain over time for a single viscoelastic material can be modeled as , where C0 is a constant that depends on Young’s modulus of the viscoelastic structure and the applied stress on the structure, and is a time constant that depends on the viscosity (η) and Young’s modulus of the viscoelastic structure [70,75]. Now, for the concept of viscoelastic behavior at the interface of the active and passive materials that we introduced in this study, we have

where G1 and G2 are constants that depend on Young’s moduli of the active and passive layers and the amount of the mismatchdriven stress generated at the interface. This stress is affected by the stimulus power (light intensity, pH value, etc.). B1 and B2 are time constants that depend on the viscosity induced at the interface (which is related to the active–passive composite) and Young’s moduli of the active and passive layers. It should be noted that Young’s modulus and viscosity are affected by the fabrication process [35] and its conditions, such as printing speed and printing resolution. In general, some of the properties of materials are affected by the printing conditions.

In addition, because the two layers are attached at the interface during the shape-morphing behavior, B1 and B2 (time constants of strains in each layer due to the mismatch-driven stress at the interface) are equal (B1 = B2 = ). Thus,

By using new uniform notations HI and HII , we have

Finally, by re-arranging, we obtain

It is the same as Eq. (20).

2.3.2. Stimulus-on versus stimulus-off

Eq. (20) is used when the stimulus is "on” and a curvature occurs in the structure. However, when the stimulus is "off,” the structure can return to its original shape by starting from the final curvature of the previous part (i.e., the stimulus-on region). Therefore, the governing equation for the second region can be obtained as:

We should clarify two points regarding Eq. (37). First, the time (t ) in the formula related to the stimulus-off region is measured from the point where the stimulus-off region begins. Otherwise, a transformation of time using t – t0 should be considered. Second, if the curvature at the beginning of the stimulus-off region (at the end of the stimulus-on region) is less than the maximum curvature, then the coefficients of the model for the stimulus-off region would be different from KI and KII. However, the format of the model will remain intact (i.e., it will be ). We should remember these two points, especially when we discuss the general graph and shape-morphing speed.

The stimulus-on and -off equations cover both reversible and irreversible materials. Fig. 3 [40,76,77] illustrates the concepts of reversible and irreversible materials. Some smart materials (whether fabricated by 4D printing or other manufacturing processes) intrinsically show reversible shape-morphing behaviors. These materials include some hydrogels [5,20,78], pH-responsive materials [79], liquid crystalline elastomers [76], and others. Some smart materials do not exhibit this intrinsic behavior. For example, SMPs are usually (but not always) irreversible (one-way). However, as mentioned in the introduction, they can be reversible if they are properly arranged with a passive material (as a multimaterial structure) [36–41].

《Fig. 3》

Fig. 3. Irreversible (one-way) and reversible (two-way) shape-morphing behaviors [40,76,77].

It should be noted that for materials that are intrinsically reversible, a self-locking mechanism can be devised by special arrangements of active and passive materials so that when the stimulus is off, the structure does not return to its original shape (if needed).

2.3.3. General graph

Based on Eqs. (20) and (37), the general graph is rendered as Fig. 4.

《Fig. 4》

Fig. 4. General graph exhibiting the time-dependent behavior of nearly all multimaterial 4D-printed (e.g., photochemical-, photothermal-, solvent-, pH-, moisture-, electrochemical-, electrothermal-, ultrasound-, enzyme-, and hydro-responsive) structures. (a) One cycle; (b) multiple cycles. Some applications require only one cycle; others require multiple cycles. In some applications, only one of the two regions of the graph is present, and in other applications, both regions are present. In some cases, the shape-morphing behavior can occur with the SME; in other cases, it can occur without the SME. It should be noted that we may need to swap the labels of "stimulus on” and "stimulus off” on this graph in some cases.

2.3.4. Validation

Here, we validate our derived bi-exponential formula. For completeness, we also consider the Timoshenko bimetal [47] and mono-exponential models. As Fig. 5 [37,80–86] shows, unlike the Timoshenko bimetal and mono-exponential models, the developed bi-exponential model perfectly captures the correct time-dependent behavior of various experimental data from previous studies, where each research group is usually expert in one specific type of stimulus or smart material (e.g., moisture or electrochemical-responsive materials). Therefore, in general, the time-dependent behavior of 4D multi-materials is nonlinear (with time) and has the specific format as given in Eq. (20). Some of the previous studies referenced in Fig. 5 [37,80–86] presented experimental data for time-dependent curvature, whereas others provided experimental data for time-dependent angle of rotation (deflection angle). However, the curvature and deflection angle have a linear relationship. Consequently, if one of them exhibits bi-exponential behavior, the other will exhibit bi-exponential behavior. All six items shown in Fig. 5 [37,80– 86] have one active and one passive material. The active component exhibits responsivity to the desired stimulus, whether with or without SME.

《Fig. 5》

Fig. 5. Validation of the proposed model by experimental data from independent studies in Refs. [37,80–84] for both on and off regions and various stimuli such as moisture [80], solvent [81], photochemical [82], photothermal [83], ultrasound [84], and heat [37]. We performed curve fitting using the Levenberg–Marquardt [85,86] method. We also applied other least-squares algorithms and obtained similar results.

We will discuss a case study as follows. Here, we show a stepby-step example regarding the bi-exponential formula derived in the third law. We generate the time-dependent behavior of the heat-responsive structure shown in Fig. 5 [37,80–86] using our bi-exponential equation in a step-by-step procedure and compare the generated behavior with the experimental data (actual timedependent behavior) proposed in that study (i.e., Ref. [37]). The study in Ref. [37] used polylactic acid (PLA) strips of 60 mm × 0.8 mm × 0.12 mm (length × width × thickness) in a bilayer structure with the paper membrane (60 mm × 0.8 mm × 0.1 mm) and a heating plate with a temperature of 90 C. The ambient temperature in those experiments was 20 C. They used a digital camera by which they measured the angle between the normals of the composite strip ends. The physics of this bilayer structure resides in Category 2 of the second law (thermal expansion/contraction).

To find HII, we have HII = , where the CTE of the active layer (printed PLA) has an average value of α ≈ 4 × 10–4 K–1 [87] (a more accurate approach is to find the effective thermal expansion coefficient in a specific case). The CTE of the paper membrane can be assumed to be zero [37], which is consistent with the assumption that the passive layer is not usually responsive under stimuli. In addition, T1( t = 20 °C, T2 = 90 °C, and the quantities and Q1 can be assumed to be zero in this case. Therefore, HII ≈ 0:028 (unitless). To find , we have , where = 1270 kg·m3 [37], cp = 1590 J·kg–1·K–1 [88], and V = 60 mm × 0:8 mm × 0:12 mm = 5.76 × 10-9 m3 . Therefore, = 0.01 J·K–1 . To calculate R1–2 , we model the heat transfer of this case study as Fig. 6, from which we obtain

where Rk is the conduction thermal resistance and Rku is the surface-convection thermal resistance. For a rectangular slab, Rk = , where Lk is the length of the slab in the direction of the heat flow, Ak is the area of the slab perpendicular to the direction of the heat flow, and k is the thermal conductivity of the slab [65]. In addition, kPLA = 0.13 W·m–1 ·K–1 [89] and kpaper = 0.05 W·m–1 ·K–1 [90]. Therefore, Rk,paper =  = 42 K·W–1 , Rk,PLA,horizontal =  = 19 K·W–1 , and Rk,PLA,vertical =  = 855 K·W–1.

For the surface-convection heat transfer, assuming thermobuoyant flow, we have [65]

where is the far-field fluid kinematic viscosity, is the fluid thermal diffusivity, g is the gravitational acceleration constant, is the fluid volumetric thermal expansion coefficient, is the far field fluid temperature, Ts is the solid surface temperature, Lku is the surface-convection characteristic length identified in Fig. 6(a), Ra is the Rayleigh number, Pr is the Prandtl number, Nu is the Nusselt number, kf is the fluid thermal conductivity, and Aku is the surface-convection area identified in Fig. 6(a). Air properties are obtained at an average fluid temperature of Tf =  = 330 K from Ref. [65], and thus = 18.37 × 10–6 m2 ·s–1 , = 26.59 × 10–6 m2 ·s–1 , Pr = 0.69, kf = 0.0287 W·m–1·K–1 , and = 0.003 K–1 . In addition, g = 9.81 m·s–2 , Lku = 0.8 × 10–3 m, and Aku = 48 × 10–6 m2 . It should be noted that a thermobuoyant flow between the horizontal heating plate and the air on top of it can be considered as well. However, without modeling the heat transfer between the heating plate and the air on top of it, we assume that after some time, the air on top of the heating plate has almost the same temperature as that of the heating plate.

By using these values, we obtain Ra = 2; which is less than 109 , and thus the flow is laminar [65]. Therefore, we use the Nusselt number formula related to laminar flow as given in Eq. (39). Furthermore, n1 = 0.5. Consequently, Nu = 1.6 and Rku = 363 K·W–1 . Therefore, R1–2 = 160 K·W–1 and finally, ≈ 1.6 s. We also model a case in which the composite strip is heated horizontally (rather than from the lateral side) so that the paper layer is between the heating plate and PLA layer. For this case, based on the physics of the new situation, we ignore the surface-convection heat transfer, but we consider the contact thermal resistance. A similar value of is obtained. Nevertheless, the first modeling is appropriate because the composite strip is heated from its lateral side in the experiment in Ref. [37] related to the bending angle measurement.

 《Fig. 6》

Fig. 6. Simplified thermal model of the case study. (a) Physical model; (b) thermal circuit model. g: gravitational acceleration constant; Ak,vertical: vertical conduction area; Ak,horizontal: horizontal conduction area; Aku: surface-convection area; Lku: surface-convection characteristic length; : far-field fluid velocity; TH: high temperature; TL: low temperature; Rk: conduction thermal resistance; Rku: surface-convection thermal resistance; Q: transferred heat.

It should be noted that the simplified thermal model shown in Fig. 6 is not unique and can be somewhat different depending on the type and level of simplifying assumptions.

To find , the following discussion is first required. PLA is a typical viscoelastic polymer. For this type of material, the time constant can be simply approximated by [70], where E '' is the loss modulus, E ' (equivalent to E) is the storage (Young’s) modulus, and f is the frequency in the dynamic mechanical analyzer (DMA) test. In a typical single-material SMP, to enable the SME, the required force (of the thermomechanical cycle) is applied externally. However, in a multi-material structure having active and passive elements, the force is generated internally between the active and passive materials. In this bilayer, we can assume that the passive layer (paper membrane) acts as a mechanism that applies the required force to the active layer (PLA). Therefore, the previous equation () used for a single-material SMP can be used for the PLA in the PLA/paper bilayer in this example. Thus,  of the general formula can be approximated by in this example. In this equation, we need the ratio of storage and loss moduli. Instead of using data from the literature that are related to molded and annealed PLA, we measure them for printed PLA. Because after the first cycle of shape shifting, a 4D-printed structure would be a treated (rather than fresh) printed structure, we go further and measure the storage and loss moduli of a treated printed PLA that has already completed the first shape-morphing cycle. The test conditions and results are illustrated in Fig. 7. To calculate the ratio between the storage and loss moduli, we consider the initial conditions; therefore, ≈  92. The value of f in the DMA test is 1 Hz, and thus ≈ 14.6 s. It should be noted that in a more accurate approach, can be obtained by  = ηinterface/Eeffective , where ηinterface is the viscosity induced at the interface of the active and passive materials, and Eeffective is the effective Young’s modulus of the active and passive layers.

《Fig. 7》

Fig. 7. DMA test. (a) DMA machine (RSA3, TA Instruments, USA); (b) test conditions; (c) storage (elastic) modulus; (d) loss modulus.

To find HI, we note that for a viscoelastic material with an elastic modulus of E under stress of σ0, we can approximate HI as σ0/E [70,75]. Now, for the viscoelastic behavior at the interface of the active and passive materials, this quantity should be adjusted to HI = σ0/Eeffective, where Eeffective depends on the E of active and passive materials, and σ0 is the mismatch-driven stress generated at the interface. Here, for PLA printed on paper, HI, which has the nature of strain as stated in the derivation of the third law, can be obtained by analyzing the strain–time curve of PLA printed on the paper membrane as measured by Zhang et al. [37]. They also measured the strain–time curve of PLA printed on a different material rather than paper [91]. In the strain–time curve presented in Ref. [37], when time tended to large numbers, the strain was approximately 0.01. This value (related to PLA) should be adjusted for the interface as follows. The stress generated at the interface can be obtained as σ0 = 0.01 × EPLA = 0.01 × 2.3 × 109 Pa = 0.023 × 109 Pa. In addition, considering the series and parallel rules for E of the active and passive layers, we obtain Eeffective = EPLA + Epaper, where Epaper = 5 × 109 Pa [37]. Therefore, HI =  ≈3 × 10-3.

The values of all the parameters in our bi-exponential equation are now set (Table 2).

《Table 2》

Table 2 Values of the parameters of our bi-exponential formula for a bilayer of PLA/paper composite.

The subscript number 1 indicates the passive layer (paper membrane) and number 2 indicates the active layer (PLA).

To convert the values of curvature to angle (bending angle), we have the following conversion factor: Angle =  = Angle = 3.4 × Curvature. The generated behavior is seen in Fig. 8 [37] and is in good agreement with the experimental data. It is also observed that both exponential terms in the bi-exponential formula are needed to capture the correct behavior, and each of them is responsible to tune one end of the whole behavior.

《Fig. 8》

Fig. 8. Comparison of theory and practice (experimental data from Ref. [37]) in the case study.

We presented three laws. Laws are different from theories. Law is the starting point and reveals something that exists. By further analysis, laws can lead to theories. Theories are usually more complicated expressions. A step-by-step example analyzed here can provide a basis for future theories in 4D printing. These laws capture and show the big picture. Currently, different approaches may exist to find and HII in a specific case, but the final bi-exponential equation in terms of and HII would be the same (Fig. 9). Future theories can examine more details of and HII  to find some general models for each of them.

 《Fig. 9》

Fig. 9. Road map of parameters. The four physical categories as well as each specific case in each category will require different approaches to obtain HI, , HII, and . Therefore, different primary parameters will appear on the road toward the final general bi-exponential equation.

《2.4. Remarks》

2.4. Remarks

We provide the following remarks for further insights, analyses, and clarifications.

Remark 1. True time-dependent behavior of 4D-printed multi-material structures. Our results show that two time constants generally govern the time-dependent shape-morphing behaviors of multi-material 4D structures. However, in some cases, the two time constants ( and ) may be approximately equal, and the time-dependent behavior can be modeled by a monoexponential equation. In some cases, the resistive (viscosity) effect at the interface of the active and passive materials (reflected in ) is negligible, and the first exponential term vanishes. In some other cases, the second exponential term may vanish. In addition, some times the two time constants are large, and the proposed bi-exponential formula tends to a linear model. In other words, if and approach large values, then , where b is a constant. This point can be realized by analyzing the related graph or by using Taylor series. Thus, in linear cases, both the proposed bi-exponential and Timoshenko bimetal models work.

Remark 2. Final shape versus instantaneous shape. We found that the Timoshenko bimetal model [47] as well as the monoexponential function cannot capture the correct time-dependent (instantaneous) behavior. However, when the final shape (maximum curvature) is achieved (t approaches large values), both the Timoshenko bimetal model [47] and our formula in Eq. (20) would be

Two points should be discussed in this regard. First, the aforementioned outcome implies that both the Timoshenko bimetal model and our formula provide similar analyses for time-independent behaviors, such as the effect of thickness on curvature. Second, some previous experimental studies reported a decrease in maximum (final) curvature with an increase in layer thickness, whereas some others reported an increase in maximum (final) curvature with an increase in layer thickness. We examine this point further by an analytical study such as follows. Eq. (40) depends on two quantities, namely, Young’s moduli ( Ei )and layer thicknesses , as h and Ii are functions of layer thicknesses (h = ; Ii = ; based on the basic assumptions of the Timoshenko bimetal model [47], the width of the strip is assumed to be small and specifically is taken as unity, as seen in Fig. 2). Therefore, we have

By analyzing Eq. (41), we find that the curvature first increases and then decreases with an increase in . Because Eq. (41) is symmetric in terms of and , the curvature exhibits a similar trend with respect to . The region in which the increase or decrease of the curvature occurs depends on the relative values of and as well as E1 and E2. In this article, we generate some possible scenarios as shown in Fig. 10. Because of the symmetry of Eq. (41), the same results as Figs. 10(a) and (b) are valid if we swap and in these two plots. A similar result to the scenario shown in Fig. 10(c) was proposed by Timoshenko [47].

《Fig. 10》

Fig. 10. Depending on the relative values of , , E1, and E2 (Eq. (41)), the relationship between curvature and layer thicknesses would be different (it can be decreasing, increasing, or a mixed behavior). (a) The relationship between curvature and thicknesses, when E1 ≤E2 and > or  < ; (b) the relationship between curvature and thicknesses, when E1 >>E2 and >  or < ; (c) the relationship between curvature and thicknesses, when E1 < or = or > E2 and = .

Remark 3. Shape-morphing speed. The shape-morphing speed is crucial in nearly any application performing dynamic intelligent behavior over time, and becomes critical in some applications such as autonomous deployment in space missions, drug delivery systems, and detection devices. Through derivatives of Eqs. (20) and (37), the magnitude of the shape-morphing speed for both on and off regions would be

As seen in Eq. (42), the shape-morphing speed is timedependent with the specific format shown. However, the Timoshenko bimetal model [47] gives a constant shape-morphing speed over time.

Based on Eq. (42), for a large amount of time (i.e., when the maximum curvature (final shape) is going to be achieved), the shape-morphing speed tends to zero for both the on and off regions. This point can also be captured from Fig. 4, where both the on and off regions become flat (constant) for t approaching large values, and the derivative of a constant function is zero.

Remark 4. Stimulus power. Here, we analyze the effect of stimulus power on the time-dependent behavior. By the expression ‘‘stimulus power,” we mean light intensity, temperature magnitude, pH value, moisture content, enzyme concentration, current magnitude, solvent concentration, and so on. By analyzing the various parameters of Eq. (20) and considering the concepts associated with the second and third laws, we find that the stimulus power will affect three parameters, HI, HII, and . Therefore, the time-dependent behaviors related to two stimulus powers would be

To illustrate the effect of stimulus power on time-dependent behaviors, we consider five stimulus powers. The related formulas would be

The general graph would be similar to Fig. 11(a). The general trend observed in Fig. 11(a) is consistent with the experimental data found in Refs. [82,83,92,93].

Some applications might exist in which a faster response without any change in the final (maximum) shape, unlike that shown in Fig. 11(a), is desirable. To this end, a smaller time constant (τ) with the same coefficient (K) is required. For five scenarios, the related formulas would be

《Fig. 11》

Fig. 11. (a) General effect of stimulus power (e.g., light intensity, pH value, and temperature magnitude) on the time-dependent behavior. This plot is based on Eq. (44); (b) tuning the response speed without changing the final shape (maximum curvature). This plot is based on Eq. (45).

The general graph would be similar to Fig. 11(b). To tune the response speed without interfering with the final shape, separate studies are needed. Carbon nanotubes [94] may be a possible solution, as they can be incorporated into stimuli-responsive materials to tune their response speed [95].

Remark 5. Extension of the developed concepts to complicated multi-material 4D structures. Based on the concepts developed in this study, the time-dependent behaviors of multi-material 4D structures can be analyzed, predicted, and tuned at nearly any level of complexity. Let us consider one simple example exhibited in Fig. 12. This case has four types of materials, two of which are active materials. For the top two layers, we have exponential terms type I and type II. Similarly, we can consider the middle two and bottom two layers. When these terms are superposed, this multimaterial structure has three exponential terms of type I and three exponential terms of type II; nevertheless, some of these exponential terms can be equal or negligible in a specific case, as discussed in Remark 1. In addition, this case will have the same general graph shown in Fig. 4, as its governing equation is a summation of exponential terms. However, the slope of its graph at a similar point would be different from that of the bi-exponential model.

《Fig. 12》

Fig. 12. 4D structure with more than two types of materials.

To realize complicated 4D structures, in addition to the multimaterial (rather than two-material) structures previously discussed, two additional important points should be considered. First, our model provides the time-dependent behavior of the curvature, which is a fundamental building block of shape shifting in multi-material structures. Other higher-level shape-morphing quantities, such as curling, twisting, coiling, and their combinations, originate from this quantity. Second, we discussed the time-dependent behavior of the curvature in one direction (narrow strip). For a plate having the same materials as those of the original narrow strip, the curvature would be the same in any direction (this point can be concluded by analyzing Ref. [47]), and can be modeled by the same two exponential terms of the original (parent) narrow strip developed in this article. For a plate that has different materials in different directions, the curvature would be different in each direction. However, the same two types of time constants and exponential terms proposed in this study can be used to model the time-dependent curvature in a specific direction accordingly. Future studies may incorporate the proposed two time constants (and exponential terms) in extensions to plates and so on.

It should also be noted that these concepts have been developed for multi-material structures that do not necessarily need to be multi-layer. In some cases (e.g., graded multi-materials with gradients), the boundary between active and passive layers may not be as clear as those shown in Fig. 2. However, in these cases, eventually, the active and passive materials will have contact in some regions, and the same concepts developed in this study will be useful. In addition, multiple stimuli may also be used in some cases. In those applications, the fundamental individual concepts developed in this study can be useful as well. As one example, a recent study [96] found the Timoshenko bimetal model useful for modeling some parameters of electro–thermo–hygro reversible 4D materials.

Remark 6. Other manufacturing processes. We distilled three laws that govern the shape-morphing behaviors of nearly all the 4D multi-materials, whether fabricated by printing and so-called 4D-printed structures or produced by other manufacturing processes.

Stimuli-responsive multi-materials can be produced by various manufacturing processes. However, AM has some benefits. First, the same reasons that motivate us to use AM for conventional (passive) materials will motive us to utilize AM for stimuli-responsive (active) materials. In other words, 4D printing conserves the general advantages of AM (such as material waste reduction, elimination of molds, dies, and machining [97], and providing complex geometries) that are not present in other manufacturing processes. Second, printing helps us elaborate on the internal structures of multi-materials precisely to enable various shapemorphing behaviors. In other words, 4D printing enables encoding local anisotropy [5] in multi-materials. In fact, before the introduction of 4D printing technology, researchers did not typically try to find a specific printing path by mathematics to yield a predictable and desired shape-morphing behavior after printing. 4D printing is a new manufacturing paradigm that combines smart materials, mathematics, and AM, as illustrated in Fig. 13 [6].

《Fig. 13》

Fig. 13. 4D printing process [6]. Readers can also see the differences between 3D and 4D printing processes shown in Ref. [6].

Remark 7. SME. As explained previously, SME is not an intrinsic property [10]. The four shape-morphing mechanisms discussed, quantified, and unified in the second law can occur with or without an SME. These four shape-morphing mechanisms are the underlying physical concepts for the relative expansion and subsequent shape-morphing behavior in 4D structures.

For example, the "heat-responsive structure” shown in Fig. 5 illustrates a typical SMP (i.e., a polymer with an SME) in a bilayer with a passive material, and its shape-morphing behavior is enabled by the relative thermal expansion of the active and passive materials [36,37,91,98]. The desired shape during the recovery step is induced by thermal expansion, and the SME determines (guides) the route of the shape change [98]. In addition, the following analyses on the SME can be useful for both single- and multi-material structures.

Bonner et al. [99] indicated that some of the previous models for SMPs are very complex and require many experiments to identify their parameters. This hinders the practical aspects of those models. However, researchers in Refs. [91,99–101] reported that the strain–time relationship of an SMP could be modeled by the exponential formula . Next, let us consider the following systematic discussion.

Usually two main approaches are used to model the behaviors of SMPs [102–109]. The first approach is thermoviscoelastic modeling; the second is based on the micromechanics of the material and is usually called phase-transition modeling [102–109].

According to the previous studies in the field of SMPs, in the thermo-viscoelastic modeling approach, the macro-scale strain of the SMPs has two parts: mechanical strain and thermal strain [102–114]. The mechanical strain, regardless of whether it is modeled by one branch or multiple branches (for multiple relaxations), is directly related to stress. The thermal strain, which is related to thermal expansion is enabled by the temperature change. During the recovery process of an SMP, the stress is zero (σ = 0) [10,115,116] (obviously, this is not a case of constrained recovery, in which the strain is kept constant. Here, a shape-shifting behavior is desired and thus the strain changes. Therefore, the other case, which is free recovery is of interest). In addition, the stress at the beginning of the recovery step is zero [10,115]. Now, consider that in the characterization studies of polymers, a creep test is used in which a constant stress σ0 is applied at time t0 and removed at time t1. Immediately after t1, a recovery of strain occurs. This shape-morphing behavior can be modeled by the Kelvin–Voigt equation as  (for the stimulus-on region, for instance). This behavior and its physics are consistent with the discussions and formulations in the molecular transformation category of the second law. However, in the thermomechanical cycle of the SMPs, the so-called recovery process is different from the recovery in a creep test of a polymer. After the programming step in the thermomechanical cycle of an SMP, we obtain a structure that has a specific shape and can maintain its shape for as long as necessary. Its external stress has already been removed and the structure is static for a finite time (in contrast to the recovery of the creep test). In addition, this structure has a residual strain (embedded during the programming step), as also mentioned by Zhang et al. [37,91]. However, it has its own properties, such as a CTE. When we apply heat to this structure during the recovery step, it will expand with its own CTE through the thermal strain. Interestingly, Zhang et al. [37] considered a similar concept implicitly. They defined an equivalent CTE, namely, αeff , and assumed that the strain of the SMP εp, after the programming step and during the recovery step was equal to its thermal strain through the equivalent CTE αeff (i.e., εp(t') = , where was the heating rate). Therefore, the strain of the SMP during the (free) recovery step in the thermomechanical cycle can be modeled by its thermal strain and there is no (need for) mechanical strain (regardless of whether it is modeled by one or multiple branches). The thermal strain has already been quantified over time in the second category of the second law and has the exponential formula εthermal(t) =C2 [1 – exp() . We should also add that for the multi-material case, the mismatch-driven mechanical stress (force) will be reflected in its related exponential term.

Readers are encouraged to read Ref. [116] for a better understanding of the thermal (and mechanical) strain and stress. For example, Lagace [116] mentioned that the term "thermal stress” is a misnomer, as it is actually "stress due to thermal effects.” Stresses are always "mechanical” [116].

Also useful is a discussion of the other modeling approach—the phase-transition modeling that is based on micromechanics modeling of the material (as mentioned by Diani et al. [102]). One of the main studies on this approach was conducted by Liu et al. [117] to model the evolution of microstructures in SMPs [76]. In addition, note that Hu et al. [76] considered the work of Liu et al. [117] as meso-scale modeling of SMPs. For small strains, Liu et al. [117] introduced an internal state variable called storage strain (stored strain) (εs), which is a "history” variable, to model phase changes. Some researchers analyzed the work of Liu et al. [117] and obtained an analytical equation for the storage strain. In their equation, the storage strain is directly proportional to stress (εs σ) [118,119]. Thus, according to their analytical model, if stress is zero in a process, εs would be expected also to be zero in that process. In addition, some other researchers such as Qi et al. [120] proposed a means of phase-transition modeling with no need to define a storage strain (stored strain).

Based on the aforementioned analyses, we did not consider SME as a specific category of fundamental physical concepts in the second law.

Remark 8. Applicability of the results of this study to singlematerial 4D structures. In the introduction, we discussed singlematerial versus multi-material structures in 4D printing and thus aimed at multi-material (active–passive) structures. However, our results can also be applied to single (active) materials. The first and second laws for single-material 4D structures would be similar to those for multi-materials because the shape-morphing behavior of nearly all smart materials is due to the expansion or contraction in the material. This expansion falls into one of the aforementioned four categories. This point is clear for most groups of smart materials and stimuli. Only for SMPs is further discussion needed. The key mechanism of the SME is the temperature-dependent mobility in SMPs [109]. In the previous remark concerning the SME, we clarified that the shape-morphing behavior of SMPs upon stimulus during recovery is (or can be) identified by the thermal strain εT, which is related to thermal expansion (one of the four categories of the second law). It should be remembered that the unified equations of the second law are all based on "strain.” It might be interesting to note that when polymer chains are programmed mechanically, a molecular transformation (chain orientation in a specific direction) will occur as explained and quantified in the molecular transformation category of the second law. However, the most important part of the shape memory behavior is the recovery step [107], as the shape shifting is obtained through this step [12]. The second law is mainly concerned with single-material 4D structures because the passive element in a multi-material 4D structure is not responsive (in general). Otherwise, it would be stimuli-responsive (active) rather than passive. Nevertheless, we used the word "relative.” Finally, the third law is a transition from single (active) materials to multi (active–passive) materials. In addition, some of the results from and discussion about the third law can also be useful for single-material 4D structures. As one example, the general graph of Fig. 4 can also be used for singlematerial SMPs. A single-material SMP can have a curvature depending on the direction in which it has been programmed. It also does not possess the exponential term related to the mismatch effect between active and passive elements. For a monoexponential equation, the general trend of on and off regions would be similar to that of a bi-exponential equation. However, the slope of its graph at a similar point would be different from that of the bi-exponential model.

In addition to the previous note, from one perspective we can say that for single-material 4D structures, we mainly need the first and second laws, and for multi-material 4D structures, we need all the three laws.

Remark 9. Scope and exceptions. Throughout this work, we did not make any specific assumption regarding the types of materials, stimuli, and length scales (macro, meso, etc.) for which these three laws are valid.

However, two points should be made about the scope of this study. First, these three laws are about shape-morphing behaviors that are currently the focus of studies in the 4D printing field. In addition, as we mentioned in the introduction, following the TED conference [1], the first embodiment of 4D printing has been entitled "4D printing: multi-material shape change” [2]. Note the word "shape” in the title of the work done by Tibbits [2]. This word has been used to introduce 4D printing. The evolutions of other properties such as color or thermal resistance have not been discussed in our study. Nevertheless, the shape-morphing behavior can alter other properties or functionalities in the desired manner. As an example, the researchers in Refs. [121,122] demonstrated smart patterned surfaces that can alter their geometries in a manner that leads to changes in their effective emissivity. This results in control of the satellite temperature without using conventional controllers and energy supplies. Second, in simple linear expansion/contraction shape-morphing behaviors, no curvature exists, where these behaviors may be considered as zero-curvature shape shifting for which the radius of curvature is infinity.

In science and engineering, laws are flexible and can have exceptions. We have reviewed a few hundred published works, including those on stimuli-responsive structures fabricated by printing or produced by other manufacturing processes. However, we could not find any exception (i.e., counter-example) to our laws. Nevertheless, we used the expression "nearly all” in the three laws for possible exceptions in the future. As a side note, we should mention that our results are general and target the fourth dimension.

Some existing works at first glance may seem to present counter-examples to these three laws. However, based on an in-depth analysis and with their fundamental physics considered, their compliance with these three laws will be understood. For example, the built-in (direct) 4D printing proposed by some researchers [123,124] has the same "underlying physics” of conventional SMPs. As we noted earlier, to enable the SME of SMPs, a mechanical force is required in addition to heat (the thermomechanical (thermo + mechanics) cycles of SMPs indicate this point as well). In the built-in (direct) SME [123,124], this mechanical force (used for programming) is provided during the printing (i.e., the printing and programming steps are integrated). A similar concept was discussed in the introduction and has also been mentioned in other studies [34,37]. The underlying physics of the shapemorphing behavior of these examples is the same as the "heatresponsive structure” shown in Fig. 5 and discussed in Remark 7.

Sometimes, a formula is derived (and validated). However, it is valid and applicable only for a specific range of cases. The biexponential formula derived and validated in this article is a universal governing equation that can model and predict the fourth dimension of nearly any 4D structure. In our study, this generality was possible because we underpinned its bases comprehensively. We used the word "law” in this article because our results are general and are also required for understanding, modeling, and predicting the fourth dimension of 4D-printed structures.

Remark 10. Future works. By itself, 3D printing is considered a multi-disciplinary field. Thus, more research areas will be involved in the 4D printing field. This diversity can increase the strength of 4D printing.

Future studies can consider several topics. One main topic is the compatibility of materials in multi-material structures. The materials should form a strong bond at their interface. Their bond should also remain strong under stimuli. Other topics include measurement and modeling of some of the parameters in the biexponential formula, such as time constants. Some parameters should be measured for active–passive materials, whereas models (whether case-specific or general) can be developed for other ones. Experimental studies can be conducted to find the exact values or the ranges of parameters for categories of materials. Other topics include software and hardware development and their integration. Future 4D printing software should have some levels of prediction. It can also provide the means of tuning the behavior over time. Future 4D printing hardware development requires some control strategies that can handle multi-material printing. Assume that ten materials are going to be encoded in a single-piece structure in various locations (voxels), with ten nozzles working simultaneously. Developing new smart materials with new applications is also a main research area. The printability of smart materials is the next topic. Many smart materials exist. However, they must become printable. Tuning the response speed is also a major topic. Still another topic is related to product development by 4D printing. 4D printing is not just a concept; it is also a manufacturing paradigm. Therefore, new products or applications that can obtain unique features by 4D printing should be regularly considered.

Remark 11. Closing. Fig. 14 presents a summary of the proposed three laws.

《Fig. 14》

Fig. 14. Summary of the proposed 4D printing laws. RH, PT, ET, UL, PC, PR, EC, and EN denote relative humidity, photothermal, electrothermal, ultrasound, photochemical, pressure, electrochemical, and enzyme, respectively. Additional stimuli and mechanisms can be indicated in each of the four physical categories.

《3. Conclusions》

3. Conclusions

We examined the fourth dimension of 4D structures and revealed three general laws that govern the time-dependent shape-morphing behaviors of nearly all 4D structures. The main purpose of the first and second laws was to understand the fourth dimension, and that of the second and third laws was to model (and predict) the fourth dimension. Through a detailed systematic qualitative and quantitative study, we derived and validated a universal bi-exponential formula that governs the fourth dimension. The results of this study can serve as general design principles for the future. They can also be incorporated into future software and hardware development in 4D printing.

It should be noted that a pure experimental study may not be able to generate a general conclusion concerning the relationship between two quantities, as it may not cover all possible regions of the relationship in various cases. Drawing a systematic conclusion is the strength of an analytical study (validated by experimental data) as conducted here.

Recently, many exciting demonstrations of 4D printing have been revealed that can hardly be achieved by other processes. However, more collaborations between scientists and engineers from various fields are needed to move from laboratory to fabrication and unveil the full potentials of 4D printing.

Historically, some topics are developed, enter research communities, and become popular; then, after a short period, they lose their broad interest. However, it is expected that 4D printing will remain attractive and useful for a long period because stimuliresponsive materials have already demonstrated their promising applications in various fields. 4D printing further helps us encode and elaborate on the structures of stimuli-responsive multimaterials by leveraging the strengths of AM (more specifically, printing) and mathematics, which are the other main elements of 4D printing.

《Acknowledgements》

Acknowledgements

Farhang Momeni wants to thank Shien-Ming Wu Manufacturing Research Center in the Department of Mechanical Engineering at the University of Michigan, Ann Arbor for support, and Andrea Poli for assistance in DMA tests in the Mechanical Testing Core directed by Ellen Arruda in the Department of Mechanical Engineering at the University of Michigan, Ann Arbor. For drawing the galactic shape of Fig. 14, we were inspired by a display designed by Rod Hill showing advancements in Reconfigurable Manufacturing Systems and installed on the wall of the Engineering Research Center (ERC) for Reconfigurable Machining Systems (RMS) at the University of Michigan, Ann Arbor. The authors hope that their results will be beneficial in ethical applications. Finally, the authors’ contributions are as follows: Farhang Momeni conceived the study, obtained and analyzed the results, and wrote the manuscript. Jun Ni discussed the results, reviewed the manuscript, and provided feedback.

《Compliance with ethics guidelines》

Compliance with ethics guidelines

Farhang Momeni and Jun Ni declare that they have no conflict of interest or financial conflicts to disclose.