A novel directional freezing based three-dimensional (3D) printing technique is applied to fabricate graphene aerogel (GA). Thermal property of the graphene ink is one of the key impacts on the material morphology and process efficiency/reliability. We develop a heat transfer model to predict temperature evolution of the printed materials and then estimate layer waiting time based on it. The proposed technique can not only improve the process efficiency and reliability but also serve as a flexible tool to predict and control the microstructure of the printed graphene aerogels. Both the simulation and experiment results demonstrate the efficiency and effectiveness of the proposed approach.

## Introduction

Graphene, an emerging two-dimensional (2D) nanomaterial that offers extraordinary mechanical, electrical, thermal, and acoustic properties [15], possesses the great potential to assemble multiscale, multifunctional three-dimensional (3D) monolith structures with widespread applications in flexible electronics, catalysis, energy storage, biomedical tissues and scaffolds, and thermal/acoustic insulation [616]. However, in order to fully unlock its exotic physicochemical properties and explore its multifunctional performance, it is highly desirable to design and fabricate engineered 3D graphene structure, such as graphene aerogel, in a controllable manner. Three-dimensional printing or additive manufacturing (AM) is an emerging technology that can fabricate physical objects directly from a computer-aided design (CAD) model without part-specific tooling and fixture, such unique property opens up great opportunities to spatially manipulate graphene nanomaterial and form them into desired engineering structures [16]. However, due to the complex graphene morphology and nontrivial material forming mechanism, it is until very recently two research groups have investigated 3D-printing graphene aerogel through micro-extrusion technique. The key concept of these techniques is to develop printable graphene ink by increasing the viscosity and form it as shearing thin non-Newtonian material [7,17,18]. Though the 3D structures are successfully printed, these approaches suffer from several drawbacks including the necessity of fillers, possible poor bonding and limited structure complexity. In order to fundamentally address these challenges, we introduced a new 3D-printing technique to fabricate 3D graphene structures with pure graphene material, better material bonding, higher structure integrity, and complex architecture printing capability. This technique is based on inkjet printing and directional freezing, where the ejected aqueous nanomaterials are frozen, and then, the growing ice crystals squeeze the graphene sheets into 3D network. Since the macrostructure can be controlled by the 3D-printing process and the microstructure can be controlled by the freezing condition and material composition, such new technique has great potential to fabricate multiscale and multifunctional 3D graphene architectures. Similar technology has been utilized to develop rapid freezing prototyping (RFP) 3D-printing process, which builds 3D ice structure by depositing and rapidly freezing water droplet layer-by-layer. Extensive work has been done based on this process [1922], however, the application of this technology is limited to visualization and molding purpose. In this research, we integrate directional freezing based 3D‐printing with freeze drying technique to form a seamlessly integrated novel process for printing multiscale multifunctional 3D graphene architectures.

As one type of crystalline material, the aqueous graphene suspension experiences phase changing during the solidification process. Therefore, the excessive latent heat must be released before the solidified material further decreases its temperature toward the environmental temperature. The buildup latent heat not only prevents proper freezing of the new material which consequently results in printing failure, it also affects the structure integrity and the material morphology due to the nonuniform temperature gradient. In order to effectively release the latent heat and super heat for ice printing, Leu's group investigated the heat transfer mechanism during the printing process with customized simulation model. The optimal waiting time was studied by using thermal finite element analysis (FEA) technique [20,23,24]. However, all the works focused on rather simple models (vertical column and vertical wall) with fixed waiting time for all the layers. The actual models for 3D printing are usually quite complicated with different cross-sectional shapes for different layers. In order to further improve the 3D-printing efficiency, it is natural to extend the constant waiting time to adaptive waiting time to accommodate the shape changes between the layers. Furthermore, in order to fundamentally understand the microstructure formation process which is primarily governed by temperature gradient, it is very desirable to develop a framework based on analytic and/or numeric thermal model to predict the solidification behavior (temperature and time) for each single droplet. Li and Leu developed FEA-based tool to simulate the heat transfer process, however, this work highly relies on commercial FEA software which lacks the flexibility to explore rather complex 3D-printing process. Furthermore, the inefficient data transfer mechanism as well as the generic algorithm itself severely increase the computational cost, both time and space. In our research work, we optimized the simulation process and developed an efficient in-house tool by using matlab programming language. Some researchers modeled other 3D-printing processes based on FEA simulation technique. Costa utilized matlab to solve 3D filament based free form extrusion process, which is also a temperature-reduction model with element mutual effect [25,26]. Shen and Le formulated models to study the thermal behavior under the effect of powder sintering in electron beam AM [27,28]. However, these research works are quite different from ours in that all these simulation models are formulated for heating-up process, where the raw material has no fixed solidification point. In our freezing-based model, the raw material experiences phase-changing process, which makes the model more complicated than other techniques. In order to address these research challenges, we propose a thermal model based FEA approach to investigate the relationship between the temperature history and fabrication process and ultimately suggest optimal settings to control and assure the process reliability, part quality, and structure integrity.

The remainder of the paper is organized as follows: We will first introduce the system setup of the directional freezing based 3D-printing process and overview the mechanism of this new technique. Second, we will discuss problem description and model formulation, and a detailed algorithm framework will also be introduced. Third, based on the analytic model and the predication framework, several test cases will be studied to demonstrate the temperature revolution and optimized waiting time, and heat transfer behavior will be analyzed based on the simulation results. Finally, we conclude and discuss the future work.

## System Setup and Motivation

### System Setup.

The experimental setup used in this project is shown in Fig. 1(a). A personal computer controls the whole system, including the three-axis motion driving system, material jetting system, and temperature and pressure control system. Besides, the computer also stores the software that can import 3D digital model and convert them into machine commands for the motion trajectory, material extrusion, and other process parameters. The material jetting system consists of a pressure regulator, multiple material reservoirs, and piezo-based jetting devices. A temperature control system is applied to the print head such that the temperature of the material is constant and just above the melting point in order to keep reliable printing property. A two-nozzle platform is setup to print graphene and supporting material (ice), as shown in Fig. 1(b). The whole system was placed in a freezer with inside temperature (−20 °C) well below the material's freezing point (∼0 °C).

For inkjet printing, the graphene nanomaterial is hydrophobic and thus segregates in water even at very low concentration unless surfactants added for their surfaces are functionalized. The segregation of graphene leads to the nozzle clogging and thus severely affects the printability. Due to the presence of hydrophilic functional groups, graphene oxide (GO) is hydrophilic and can be easily dispersed in water at relatively high concentrations. Although GO is not electrically conductive, it can be thermally, chemically, and photothermally reduced to graphene. Based on the excellent properties stated above, GO is suitable to serve as raw material. Low concentration (10 mg/ml) aqueous GO ink was selectively ejected onto an aluminum platform with a certain amount of flow rate and moved following a predefined path. The ink is then frozen through the heat conduction with the platform and heat convection with the ambient. During the frozen process, the phase changing and separation can force the nonfreezable solid nanomaterials (GO) to accumulate between the growing ice crystals. The nanomaterials trapped by the ice crystals form bridges between the crystals. If the loading of nanomaterial is enough, the entrapped nanomaterials will form a continuous 3D network molded by crystals. The 3D-printed “green” part will then be freeze dried and the water will be removed to achieve the final graphene aerogel (GA) with porous structures. Figures 1(c)1(e) show multiple 3D-printed GA structures. The interested readers are referred to our previous work for detailed explanation of the 3D-printing GA process [29].

### Motivation of the Proposed Work.

The aqueous GO ink in the aforementioned 3D-printing process is a low-viscosity Newtonian fluid, allowing it to be printed in a drop-on-demand (DOD) mode, where the material is jetted drop-by-drop only if needed. These discrete droplets serve as the basic elements in jetting-based printing process. The discrete ejection mode differs our technique from other continuous 3D GA printing approaches. The DOD-based technique can achieve much higher accuracy and flexibility. In this process, the layer bonding and structure forming are based on the freezing-induced solidification; therefore, the stability of the ink solidification and the efficiency of the printing process are two crucial but contradictory criteria.

The solidification of the newly deposited droplet is powered by the conduction with the previous solidified layers and the convection with the low-temperature environment. Because the phase-changing process takes non-negligible time due to the release of latent heat, the newly deposited droplet will not be properly solidified if the previous layer is not yet frozen. Therefore, certain amount of waiting time is desired to solidify the most recent layer of material, which otherwise will damage and even fail the printing process. Figure 2 shows a failure case where inadequate waiting time is applied. Theoretically, the longer waiting time is preferable to allow for the complete heat dissipation and assure the reliability of the solidification process, however, longer waiting time, on the other hand, prolongs the fabrication time. In addition, the aqueous GO ink is prone to dry or freeze at the tip of the printing nozzle and subsequently clogs the print head if the waiting time is too long. Thus, the optimal waiting time is expected to balance the process reliability, efficiency, and part quality. The thermal management, therefore, plays an important role in controlling the macroscopic structure of the 3D-printed part.

The final product of the proposed 3D-printing process is graphene aerogel, which is obtained by freeze drying the 3D-printed ice green part. According to the mechanism of the freeze drying process, the integrity and morphology of the graphene aerogel are highly dependent on the temperature gradient. At low temperature, the nucleation rate of the ice is much higher than the crystal growth rate, thus crystals do not have enough time and space to grow and a large number of small ice crystals form. At higher freezing temperatures, crystal growth dominates and larger pores form after freeze drying [30,31]. Therefore, besides the important role in controlling the macroscopic structure, the thermal management is also critical to control the microscopic structure (pore size, density, and morphology) of the 3D-printed GA part.

## Problem Description and Approach Overview

### Problem Description.

As discussed in the previous section, “Motivation of the Proposed Work,” certain waiting time is expected between two adjacent layers in order to assure the printing quality and success rate. However, balance has to be made between the process efficiency and reliability. Therefore, finding out the optimal waiting time for a specific print is the key research question. The major criterion of the waiting time is to assure that the previously deposited layer is completely solidified and/or reaches certain low-temperature level. In the jetting-based 3D-printing process, the new droplet is ejected from the nozzle, lands on the ice surface, deforms, spreads, freezes, and bonds with previous layer through hydrogen bonding. The driving force of solidification is from the low-temperature ambient through convection and the cold ice surface through conduction. The thermal behavior and thermal management are the keys to identify the optimal waiting time in this dynamic process. However, as the materials are deposited sequentially and the layer geometry varies for typical parts, it is very challenging if not impossible to solve this problem through analytic approach with closed-form formulation. Instead, numerical approach such as FEA based modeling and simulation can effectively characterize and analyze the complex dynamic process, we will therefore investigate the thermal behavior of the directional freezing based 3D-printing process through FEA simulation based on the thermal modeling described in the following sections, “Process Modeling” and “Parameter Settings and Algorithm.”

### Process Modeling.

The process modeling for thermal analysis is based on the heat transfer theory, and the general governing equation [32] is
$∇2T+q˙k=1α∂T∂t$
(1)
where $T$ is the temperature, $q˙$ is the internal heat generation rate, $k$ is the material conductivity, and $α=k/ρc$ is the thermal diffusivity of the material and related to material conductivity $k$, material density $ρ$, and specific heat $c$. In this case, no internal heat is generated inside the elements, thus $q˙$  = 0 and Eq. (1) is simplified as
$∇2T=1α∂T∂t$
(2)

In Cartesian coordinate system, the Laplacian term $∇2T$ can also be expressed as $(∂2T/∂x2)+(∂2T/∂y2)+(∂2T/∂z2)$, which represents heat flux from three directions into the current unit volume.

Shape simplification: Solidified part is built up by ink droplets, which are deposited drop-by-drop. When a droplet is deposited, it will spread and contact with the adjacent droplets. Since the spreading speed is much higher than the freezing speed [20], the droplet can be modeled as thin rectangular slab shape without loss of accuracy. Considering the small volume (<100 pL) of each droplet, we believe this shape simplification is valid and it can significantly improve the computational efficiency.

Thermally smallelement: Even with the shape simplification, calculating the temperature for every element (droplet) is still computationally challenging. In addition, solving Eqs. (1) and (2) is not trivial. Based on our observation, homogenizing the internal temperature for each element can significantly reduce the computational cost and increase the simulation efficiency. If the heat conduction inside the element body is significantly faster than its surface heat exchange with outside environment, we can consider the element to be thermally small. To verify this assumption, Biot number is introduced in this work. Biot number is a dimensionless number which indicates a ratio of heat transfer resistance inside a body and on the surface of the body. If a Biot number is much smaller than 1, heat transfer resistance inside the body is much smaller than heat transfer resistance on the surface, which also means heat transfers much faster inside the body than on the surface. Generally, if Biot number is smaller than 0.1, the body then can be considered “thermally simple” or “thermally small”. The Biot number is defined as the following equation [33]:
$Bi=hLckb$
(3)

where h is the surface heat transfer coefficient. For frozen GO ink, h is around 6–30 . In this case, the contact heat transfer is considered as one type of the film heat transfer, we approximately set it as a relatively large value, 50 $W/(m2 K)$. $Lc$ is the characteristic length, a typical value is 0.05 mm based on the real size of a droplet. $kb$ is the thermal conductivity of the body. For liquid GO ink, thermal conductivity is around 0.6 $W/(m K)$, thus it follows that $Bi$  = 0.00416 $≪$ 0.1. For frozen GO ink, thermal conductivity is around 2.0 $W/(m K)$; therefore, the corresponding $Bi$  = 0.00125 $≪$ 0.1. Based on this calculation, one droplet, or one element, can be reasonably considered to be thermally small.

With the thermally small property, the internal conduction can be ignored without losing computational accuracy. Therefore, $(∂2T/∂x2)+(∂2T/∂y2)+(∂2T/∂z2)$ can be replaced by $∑MλiAiqi″$, as they both represent heat flow into the element in the equivalent way. Equation (2) can be further simplified as
$ρcδVdTdt+∑MλiAiqi″=0$
(4)
where $T$ is the temperature, $t$ is the time, $ρ$ is the density of the ink, $c$ is the specific heat of ink, $δV$ is the volume of a liquid droplet, and $λi$, $Ai$, and $qi″$ are the boundary condition (BC) coefficient, working area, and heat transfer of the $i$th surface, respectively. Equation (4) can be expanded into specific heat transfer equation as follows:
$ρcδVdTdt+∑M1λiAihsur(T−Ta)+∑M2λjAjhn(T−Tn)+∑M3λkAkhs(T−Ts)=0$
(5)
where $Ta$ is the ambient temperature, $hsur$ is the surface heat transfer coefficient, including several types of film heat transfer, such as convection, radiation, and evaporation, $hn$ is the contact heat transfer coefficient with the $n$ th adjacent element, $hs$ is the reciprocal of thermal contact resistance between elements and substrate, $Tn$ is the temperature of the $n$th element, and $Ts$ is the temperature of substrate, which is identical with the ambient temperature in this simulation. Equation (5) can be further derived into
$ρcδVdTdt+C1T−C2=0$
(6)
where the constant numbers are denoted as $C1$ and $C2$ for the sake of simplicity. Solving Eq. (6), we obtain
(7)
where $tc$ is the exact time at which the current droplet is deposited. When calculating the temperature of an element, the temperatures of adjacent elements keep unchanged within one step. Based on this approximation, Eq. (7) can be further written as
$Ti+1=(Ti−C2C1)exp[−C1·Δt]+C2C1$
(8)

where $Ti+1$ and $Ti$ are the temperature in current and previous step, respectively, and $Δt$ is the step size. Equation (8) is the final explicit analytic solution and can be directly utilized for temperature update during the iteration-based dynamic thermal evolution process.

### Parameter Settings and Algorithm.

During the printing process, a moving printing head continuously deposits GO droplet on a cold plate or previous cold layer. The newly deposited droplet can be modeled as a moving heat source in the simulation, and thus, the BCs have to be iteratively updated when a new droplet is deposited. Because of considerably small volume, each droplet being modeled as one element can still have a reasonable accuracy. The aforementioned shape simplification is shown as follows: the two-dimensional (2D) model of one element can be simplified as rectangular shape as shown in Fig. 3. With such simplification, it is more computationally efficient to define and update the BCs of the model. For all the elements, only three types of heat transfer modes need to be considered as shown Fig. 3: (1) heat transfer with substrate (Fig. 3(a)); (2) heat transfer with adjacent elements (Fig. 3(b)); and (3) heat transfer with ambient (Fig. 3(c)). As substrate has a considerably large conductivity compared with ink, the boundary condition for the bottom of first layer is set as constant temperature same as the ambient temperature.

The initial conditions, material properties, and process parameters are listed in Tables 13, respectively. Based on the heat transfer model and parameter settings, an efficient thermal FEA-based algorithm is developed to study the heat transfer mechanism of the directional freezing 3D-printing process. The algorithm is illustrated in the following framework (Fig. 4). The algorithm takes the tool path data based on the input geometry and converts them into information for the FEA model. After setting the material properties, process parameters, and the initial conditions of every element, substrate, and the ambient, the algorithm updates the temperature of each element (droplet) iteratively until the entire printing process is completed.

At the beginning of each iteration, one droplet (element) will be deposited and the related BCs will be updated unless the layer-change condition is met. After that, the system will judge whether the current element is in its phase-changing process, if not, the temperature of every deposited element will be simply computed and updated. If the current element is under phase-changing process, the latent heat of the element will be updated and the next iteration starts. If the current layer is completed, i.e., the last droplet of the layer has been deposited, the waiting time will be reset and updated and the criterion of “whether it is ready to print next layer” will be evaluated. We set a target temperature for all the deposited elements especially the ones in the most recently deposited layer, when the temperature of all the deposited elements reaches the threshold, the system is ready for the next layer. If it is ready to print next layer, the waiting time register will be reset and a new layer will be ready for process. If not, the temperature of all the elements will be computed for another step time and the waiting time will be updated. The whole process will terminate when the preset step number is met, normally the whole simulation process consists of the whole material building process and a postcooling process. Note that the liquid ink and solid frozen ink have different material properties, and their material properties are temperature dependent, when the temperature of an element is updated, its material properties need to be updated if necessary.

## Results and Discussion

The aforementioned algorithm is implemented and tested in matlab environment, with parallel computing strategy applied to improve computation efficiency. Both 2D and 3D models with different geometries are tested. The effect of the waiting time is studied for various cases including no waiting time, constant waiting time, and adaptive waiting time. The detailed results are discussed in the following sections, “Two-Dimensional Model Study” and “Three-Dimensional Model Study.”

### Two-Dimensional Model Study.

Case 1: No waiting time study. A simple 4 × 4 square shape was first used to study the heat transfer during the printing process and verify the predictive model. The tool path is shown in Fig. 5. Four vertically aligned elements (1, 8, 9, and 16) are selected to study the heat transfer and temperature evolutions. In this study, no waiting time is assigned between two consecutive layers.

In Fig. 6, four curves are plotted to indicate the temperature evolution of element 1, 8, 9, and 16 within 1 s. During this time period, only elements 1 and 8 completed phase-changing process. Latent heat releasing process has considerable impact on temperature evolution of adjacent elements. For instance, temperature of element 1 keeps at around −16 °C before element 8 released all its latent heat during phase-changing stage. After element 8 released all the latent heat, temperature of element 1 continues to drop to a lower value, which indicates the influence from element 8 on element 1. Besides, temperature of element 1 drops more quickly than elements above it, showing the influence from substrate to element 1 is more significant than elements above. Other interesting results can also be found in Fig. 6. When droplet 16 is deposited on top of droplet 9, because of the high temperature of the new droplet, temperature of droplet 9 starts to increase above 0 °C during its phase-changing stage, indicating that droplet 9 starts to melt. This case is common in real fabricating process. Figure 2 shows a real melting case, several top layers melt during printing, resulting in a rough surface and inevitable imprecision.

Case 2: Constant threshold temperature waiting time. To avoid melting during fabrication process, one straightforward way is to wait until the deposited materials are completely solidified and/or reach to a certain threshold temperature (lower than freezing point) before a new layer is deposited. Simulation is performed to test the effect of waiting time before depositing new layers using the same test case. Figure 7 shows the results of same test case with a threshold temperature of −19 °C in Fig. 7. Clearly, no melting issues are observed as existed materials are more prepared for a new layer.

Table 4 shows the waiting time for every layer. As can be seen from the table, to reach the same threshold temperature, it takes less and less waiting time as the layer number increases. The major reason is the effect of the heat conduction with the platform is attenuating as the part grows, and the thermal conductivity of the frozen ink is much lower than the metal platform.

As waiting time is inevitable for every layer except for the first one, it becomes a critical factor for the fabrication efficiency of freezing-based 3D-printing process. In this work, waiting time is studied with the help of heat transfer simulation model. In order to better understand the property of waiting time, the previous 4 × 4 case is extended to 40 × 4 case. The simulation results of the waiting time are shown in Fig. 8.

As can be seen from Fig. 8, the waiting time stays at an approximately constant value after the number of layers reaches a certain level. In this case, the waiting time stays as 4.5220 s after layer #35 is deposited. As the temperature of substrate is considered as a constant at relatively low value, the influence from substrate to temperature evolution of the first deposited layers is significant. But as the layer number increases, this influence is attenuating, and after certain layers (#35 in this case) this influence disappears and the waiting time for the following layers will keep unchanged. We believe the waiting time is related to layer thickness, ambient temperature [24], and the geometry and dimensions of the part.

Case 3: Adaptive waiting time. Constant waiting time cannot assure both the process reliability and efficiency. An adaptive waiting time is desirable to maximize the fabrication efficiency without losing the reliability. A 20 × 20 square case, as displayed in Fig. 9(a), is simulated and the waiting time is collected for every layer. In this case, the threshold temperature of depositing a new layer is still −19 °C. Figure 10(a) shows the relationship between the waiting time and layer number. This plot is evidently different from the first 20 layers plot line in Fig. 8, verifying our hypothesis that the part dimension affects the waiting time.

Figure 10(b) shows the waiting time trend for the 20 × 20 right triangle in Fig. 9(b). Different from the square shape cases, the waiting time decreases after layer #15, which is mainly caused by the decrease of layer dimension. In conclusion, without geometry change, the waiting time will gradually reach a constant value with the attenuation of substrate effect. However, if there is geometry change, the waiting time will vary accordingly after certain number of layers. For other more complicated cases, the waiting time is expected to vary depending on layer location, ambient temperature, layer thickness, and the part geometry and dimension.

From the above scenarios, we can conclude that for the first several layers, the impact of substrate is the most significant factor, and the waiting time linearly increases with the layer number. After certain number of layers, the effect of substrate will gradually attenuate and eventually disappear. In order to further verify the impact of substrate, two models with same top half part (Fig. 11) are simulated. Figure 12(a) clearly shows the influence of the substrate for the case in Fig. 11(a). The waiting time continuously increases till layer #16 and then decreases because of the reduction of layer size. The part in Fig. 11(b) has 12 more layers duplicating the first layer compared with the part in Fig. 11(a), and the simulated result is shown in Fig. 12(b). As can be seen from Fig. 12(b), the first several layers are significantly influenced by the heat drainage of substrate, and the waiting time increases linearly with the number of layers. From layer #20, the waiting time suddenly decreases because of sudden reduction of layer size, and it then stays at approximately constant level because of constant layer size. Finally, the waiting time gradually decreases due to the decrease of the layer size. The comparison between these two test cases clearly demonstrated that the waiting time is influenced by the part geometry and the cold substrate.

### Three-Dimensional Model Study

#### Model Formulation.

The 2D model in the previous section, “Two-Dimensional Model Study,” is effective to investigate and understand the thermal behavior of the freezing-based printing process. In practice, the real parts are printed in 3D space, and it is more likely that layers printed with different toolpath and printing strategies have different optimal waiting times. However, it is computationally more expensive to simulate the 3D model than the 2D model. For the sake of computational simplicity, a typical cubic model is formulated to study the thermal behavior of the 3D model. Similar to 2D model, the element shape is simplified as flat cube for 3D model. The droplet (element) will be deposited one by one with the predefined printing parameters and tool path to build the large cubic part, as shown in Fig. 13.

The thermal evolution equations (6)(8) are also applicable to 3D models. The major difference between 3D and 2D model is the heat transfer conditions and boundary conditions, which is updated according to toolpath strategy, speed, and waiting times. In this model, as every layer is same in size and geometry, toolpath strategies for each layer will be identical.

#### Simulation Results.

Figure 14 shows the history of the waiting time for the 3D cubic model. It can be clearly seen that the waiting time for 3D model shows the same trend of “increase-slow-stable” as its counterpart 2D model and the substrate effect also plays a crucial role in the solidification process of the first few layers. It is expected that the waiting time is likely to vary with different layer geometries and tool path strategies.

In Sec. 2, Fig. 2 demonstrated a failed test case caused by inadequate waiting time during fabrication process. In order to obtain an acceptable part, the proposed thermal model is applied and the optimal waiting time is identified for each layer. The waiting time is then applied to the actual printing process, and the fabricated GA part (Fig. 15) shows substantially improved quality and reliability.

## Conclusions

Three-dimensional graphene aerogel has widespread applications. We presented a novel 3D-printing technique based on directional freezing to fabricate complex graphene aerogel architectures with multifunctional and multiscale properties. The heat transfer is vital to the structural integrity and morphology of the graphene structure, as well as the reliability and efficiency of the process. This paper studied the temperature evolution and explored the adaptive waiting time during printing. An analytic and numeric model is formulated and implemented in matlab to dynamically predict the temperature of each droplet and waiting time between consecutive layers. The experimental results show that the temperature evolution and waiting time of the first several layers are significantly affected by substrate. After a certain number of layers, the effect of substrate on the waiting time gradually attenuated. For a uniform geometry, the waiting time stays constant after the first several layers. While for the complex geometry, the shape of the layers significantly affects the waiting time. Our proposed technique can efficiently and robustly predict the waiting time for any geometry, which consequently reduces the fabrication time. Both the simulation and experimental results demonstrated the effectiveness of the proposed technique. Future work includes integrating the predictive model into tool path planning framework to further optimize the efficiency and reliability.

## Acknowledgment

The authors acknowledge the seed funds support from the Sustainable Manufacturing and Advanced Robotics Technologies Community of Excellence (SMART CoE) and the start-up funds from the State University of New York at Buffalo.

## References

References
1.
Balandin
,
A. A.
,
Ghosh
,
S.
,
Bao
,
W.
,
Calizo
,
I.
,
Teweldebrhan
,
D.
,
Miao
,
F.
, and
Lau
,
C. N.
,
2008
, “
Superior Thermal Conductivity of Single-Layer Graphene
,”
Nano Lett.
,
8
(
3
), pp.
902
907
.
2.
Lee
,
C.
,
Wei
,
X.
,
Kysar
,
J. W.
, and
Hone
,
J.
,
2008
, “
Measurement of the Elastic Properties and Intrinsic Strength of Monolayer Graphene
,”
Science
,
321
(
5887
), pp.
385
388
.
3.
Neto
,
A. C.
,
Guinea
,
F.
,
Peres
,
N.
,
Novoselov
,
K. S.
, and
Geim
,
A. K.
,
2009
, “
The Electronic Properties of Graphene
,”
Rev. Mod. Phys.
,
81
(
1
), p.
109
.
4.
Novoselov
,
K. S.
,
Geim
,
A. K.
,
Morozov
,
S.
,
Jiang
,
D.
,
Zhang
,
Y.
,
Dubonos
,
S. A.
,
Grigorieva
,
I.
, and
Firsov
,
A.
,
2004
, “
Electric Field Effect in Atomically Thin Carbon Films
,”
Science
,
306
(
5696
), pp.
666
669
.
5.
Geim
,
A. K.
, and
Novoselov
,
K. S.
,
2007
, “
The Rise of Graphene
,”
Nat. Mater.
,
6
(
3
), pp.
183
191
.
6.
Cong
,
H.-P.
,
Wang
,
P.
, and
Yu
,
S.-H.
,
2013
, “
Stretchable and Self-Healing Graphene Oxide–Polymer Composite Hydrogels: A Dual-Network Design
,”
Chem. Mater.
,
25
(
16
), pp.
3357
3362
.
7.
Jakus
,
A. E.
,
Secor
,
E. B.
,
Rutz
,
A. L.
,
Jordan
,
S. W.
,
Hersam
,
M. C.
, and
Shah
,
R. N.
,
2015
, “
Three-Dimensional Printing of High-Content Graphene Scaffolds for Electronic and Biomedical Applications
,”
ACS Nano
,
9
(
4
), pp.
4636
4648
.
8.
Leigh
,
S. J.
,
,
R. J.
,
,
C. P.
,
Billson
,
D. R.
, and
Hutchins
,
D. A.
,
2012
, “
A Simple, Low-Cost Conductive Composite Material for 3D Printing of Electronic Sensors
,”
PLoS One
,
7
(
11
), p.
e49365
.
9.
Maiti
,
U. N.
,
Lim
,
J.
,
Lee
,
K. E.
,
Lee
,
W. J.
, and
Kim
,
S. O.
,
2014
, “
Three-Dimensional Shape Engineered, Interfacial Gelation of Reduced Graphene Oxide for High Rate, Large Capacity Supercapacitors
,”
,
26
(
4
), pp.
615
619
.
10.
Menzel
,
R.
,
Barg
,
S.
,
Miranda
,
M.
,
Anthony
,
D. B.
,
Bawaked
,
S. M.
,
Mokhtar
,
M.
,
Al-Thabaiti
,
S. A.
,
Basahel
,
S. N.
,
Saiz
,
E.
, and
Shaffer
,
M. S.
,
2015
, “
Joule Heating Characteristics of Emulsion-Templated Graphene Aerogels
,”
,
25
(
1
), pp.
28
35
.
11.
Wicklein
,
B.
,
Kocjan
,
A.
,
Salazar-Alvarez
,
G.
,
Carosio
,
F.
,
Camino
,
G.
,
Antonietti
,
M.
, and
Bergström
,
L.
,
2015
, “
Thermally Insulating and Fire-Retardant Lightweight Anisotropic Foams Based on Nanocellulose and Graphene Oxide
,”
Nat. Nanotechnol.
,
10
(
3
), pp.
277
283
.
12.
Xu
,
X.
,
Li
,
H.
,
Zhang
,
Q.
,
Hu
,
H.
,
Zhao
,
Z.
,
Li
,
J.
,
Li
,
J.
,
Qiao
,
Y.
, and
Gogotsi
,
Y.
,
2015
, “
Self-Sensing, Ultralight, and Conductive 3D Graphene/Iron Oxide Aerogel Elastomer Deformable in a Magnetic Field
,”
ACS Nano
,
9
(
4
), pp.
3969
3977
.
13.
Ye
,
S.
,
Feng
,
J.
, and
Wu
,
P.
,
2013
, “
Highly Elastic Graphene Oxide–Epoxy Composite Aerogels Via Simple Freeze-Drying and Subsequent Routine Curing
,”
J. Mater. Chem. A
,
1
(
10
), pp.
3495
3502
.
14.
Vickery
,
J. L.
,
Patil
,
A. J.
, and
Mann
,
S.
,
2009
, “
Fabrication of Graphene–Polymer Nanocomposites With Higher-Order Three-Dimensional Architectures
,”
,
21
(
21
), pp.
2180
2184
.
15.
Estevez
,
L.
,
Kelarakis
,
A.
,
Gong
,
Q.
,
Da'as
,
E. H.
, and
Giannelis
,
E. P.
,
2011
, “
Multifunctional Graphene/Platinum/Nafion Hybrids Via Ice Templating
,”
J. Am. Chem. Soc.
,
133
(
16
), pp.
6122
6125
.
16.
Bourell
,
D. L.
,
Leu
,
M. C.
, and
Rosen
,
D. W.
,
2009
, “
,” The University of Texas at Austin, Austin, TX, pp. 11–15.
17.
García-Tuñon
,
E.
,
Barg
,
S.
,
Franco
,
J.
,
Bell
,
R.
,
Eslava
,
S.
,
D'Elia
,
E.
,
Maher
,
R. C.
,
Guitian
,
F.
, and
Saiz
,
E.
,
2015
, “
Printing in Three Dimensions With Graphene
,”
,
27
(
10
), pp.
1688
1693
.
18.
Zhu
,
C.
,
Han
,
T. Y.-J.
,
Duoss
,
E. B.
,
Golobic
,
A. M.
,
Kuntz
,
J. D.
,
,
C. M.
, and
Worsley
,
M. A.
,
2015
, “
Highly Compressible 3D Periodic Graphene Aerogel Microlattices
,”
Nat. Commun.
,
6
, p. 6962.
19.
Barnett
,
E.
,
Angeles
,
J.
,
Pasini
,
D.
, and
Sijpkes
,
P.
,
2009
, “
Robot-Assisted Rapid Prototyping for Ice Structures
,”
IEEE International Conference on Robotics and Automation
,
ICRA'09
, May 12–17, pp.
146
151
.
20.
Bryant
,
F. D.
, and
Leu
,
M. C.
,
2009
, “
Predictive Modeling and Experimental Verification of Temperature and Concentration in Rapid Freeze Prototyping With Support Material
,”
ASME J. Manuf. Sci. Eng.
,
131
(
4
), p.
041020
.
21.
Ossino
,
A.
,
Barnett
,
E.
,
Angeles
,
J.
,
Pasini
,
D.
, and
Sijpkes
,
P.
,
2009
, “
Path Planning for Robot-Assisted Rapid Prototyping of Ice Structures
,”
Trans. Can. Soc. Mech. Eng.
,
33
(4), pp.
689
700
.
22.
Zhao
,
X.
,
Landers
,
R. G.
, and
Leu
,
M. C.
,
2010
, “
Adaptive Extrusion Force Control of Freeze-Form Extrusion Fabrication Processes
,”
ASME J. Manuf. Sci. Eng.
,
132
(
6
), p.
064504
.
23.
Sui
,
G.
, and
Leu
,
M. C.
,
2003
, “
Thermal Analysis of Ice Walls Built by Rapid Freeze Prototyping
,”
ASME J. Manuf. Sci. Eng.
,
125
(
4
), pp.
824
834
.
24.
Liu
,
Q.
, and
Leu
,
M. C.
,
2007
, “
Finite Element Analysis of Solidification in Rapid Freeze Prototyping
,”
ASME J. Manuf. Sci. Eng.
,
129
(
4
), pp.
810
820
.
25.
Costa
,
S.
,
Duarte
,
F.
, and
Covas
,
J. A.
,
2011
, “
Using MATLAB to Compute Heat Transfer in Free Form Extrusion
,”
MATLAB—A Ubiquitous Tool for the Practical Engineer
,
C. M.
Ionescu
, ed., InTech, Rijeka, Croatia, p.
453
.
26.
Costa
,
S.
,
2013
, “
Free Form Extrusion: Extrusion of 3D Components Using Complex Polymeric Systems
,”
Ph.D. thesis
, Universidade de Minho, Braga, Portugal.
27.
Bellini
,
A.
,
Shor
,
L.
, and
Guceri
,
S. I.
,
2005
, “
New Developments in Fused Deposition Modeling of Ceramics
,”
Rapid Prototyping J.
,
11
(
4
), pp.
214
220
.
28.
Shen
,
N.
, and
Chou
,
K.
,
2012
, “
Thermal Modeling of Electron Beam Additive Manufacturing Process: Powder Sintering Effects
,”
ASME
Paper No. MSEC2012-7253.
29.
Zhang
,
Q.
,
Zhang
,
F.
,
Medarametla
,
S. P.
,
Li
,
H.
,
Zhou
,
C.
, and
Lin
,
D.
,
2016
, “
3D Printing of Graphene Aerogels
,”
Small
,
12
(
13
), pp.
1702
1708
.
30.
Moner-Girona
,
M.
,
Roig
,
A.
,
Molins
,
E.
, and
Llibre
,
J.
,
2003
, “
Sol-Gel Route to Direct Formation of Silica Aerogel Microparticles Using Supercritical Solvents
,”
J. Sol-Gel Sci. Technol.
,
26
(
1
), pp.
645
649
.
31.
Baumann
,
T. F.
,
Gash
,
A. E.
,
Chinn
,
S. C.
,
Sawvel
,
A. M.
,
Maxwell
,
R. S.
, and
Satcher
,
J. H.
,
2005
, “
Synthesis of High-Surface-Area Alumina Aerogels Without the Use of Alkoxide Precursors
,”
Chem. Mater.
,
17
(
2
), pp.
395
401
.
32.
Kakac
,
S.
, and
Yener
,
Y.
,
1993
,
Heat Conduction
,
Taylor and Francis
,
Washington, DC
, Chap. 2.
33.
Bergman
,
T. L.
,
Incropera
,
F. P.
, and
Lavine
,
A. S.
,
2011
,
Fundamentals of Heat and Mass Transfer
,
Wiley
,
Danvers, MA
, pp.
260
261
.