Analytical Study and Experimental Verification of Shear-Thinning Ink Flow in Direct Ink Writing Process

The direct ink writing (DIW) process is one of the most broadly used 3D printing techniques. It is capable of assembling feedstock materials to complex 3D shapes on a digitally defined deposition path. DIW is a versatile and cost-effective additive manufacturing process, offering unparalleled opportunities in fabricating multiscale and multifunctional products. Although DIW has easy accessibility and the process is straightforward, interests in material variability, shape complexity, printing resolution, and structure fidelity are continuously leading the ever fast development of the DIW process [1]. In fact, the rapid advancement of DIW has been heavily adopted across many research areas and found diverse applications of the printed samples, including energy storage [2,3], catalyst [4], bone repair [5], artificial organ [6], water purification [7], and thermal insulation [8–10].

Controlling the extrusion performance of ink is a complex challenge. Two factors are of importance for successful DIW, the ink printability and the printing resolution. Inks with high printability shall possess non-Newtonian shear-thinning properties to ensure continuous filament extrusion and maintain the free-standing 3D shape after printing. The shear-thinning property depicts that the ink viscosity decreases with the increasing shear rate. Shear rate, as a measure of the deformation rate of the ink, increases as the ink is extruded along the axial direction of the syringe (Fig. 1(b)). Shear stress is induced during the extrusion process, and the relation between the shear stress and the shear rate is a downward concave curve or called a flow curve (Fig. 1(c)). The slope of the flow curve is the viscosity, which decreases with the increasing shear rate. In the extrusion process, shear stress elongates the structural units (e.g., untangled polymer chains, stretched fibers, deformed bubbles) within the colloidal inks and thus originates the shear-thinning behavior. For instance, particle suspension is one of the widely used DIW inks. At a quiescent state, the particles within the ink could be either aggregated due to van der Waals attraction or randomly separated by the repulsive charges on surfaces. In addition, fibers within the ink could orient randomly, following the minimum energy principle [11]. However, when subjected to shear stress, the particles are disintegrated from the particle–particle interaction, and the fibers are stretched along the flow direction. As a result, these changes facilitate the shear flow, thus lowering the viscosity with the increasing shear stress [12]. Guo and Zhou [1] comprehensively summarized the ink composition and the resulting rheology property, providing an effective ink formulation strategy. A rotary rheometer is often used to characterize the shear-thinning property. However, the viscosity data are not always applicable to inks in the real DIW process. On the other hand, the in situ and real-time measurement of the viscosity–shear rate relation in DIW is still a challenge. Such investigation cannot be carried out using approaches like particle imaging velocimetry. The use of tracing particles does not provide sufficient resolution for numerical studies in a three-dimensional space [13]. Furthermore, increasing the tracing particle concentration could change the nature of the ink dynamics. As such, understanding the DIW process and advancing the ink extrusion process is rigorous science and relies on gathering sufficient data from experimental and analytical techniques.

As the computational power rapidly increases and more empirical models are established, computational fluid dynamics (CFD), a finite element analysis approach, becomes a powerful tool to elucidate the physical phenomena during the shear-thinning extrusion process. In fact, researchers rely on simulation results to guide the physical experiments. For example, simulation studies can reveal the shear rate and shear stress of non-Newtonian fluid and provide essential information to prevent shear-induced damage during oil drilling operations [14] or affect the cell viability in extrusion-based bio-fabrication [15]. The numerical simulation is also employed to study the printability of the food DIW process. Yang et al. [16] identify that the inlet volume flowrate is a primary parameter to determine the extrusion velocity field. Liu et al. [17] and Guo et al. [18] find that the syringe-based process is favorable for continuous and stable extrusion. On the other hand, the screw-based DIW process could suffer an overly high shear rate and incur issues such as backflows and discontinued extrusion. Although these results are informative and offer qualitative instructions, the printing resolution is not of a major concern. As a matter of fact, the printing resolution and shape fidelity in both the lateral and vertical directions contribute to acquiring the high performance of the printed samples. However, to the best of our knowledge, comprehensive simulation works have yet to be carried out to exhibit the ink flow behavior in the DIW process. Therefore, a comprehensive simulation study of the DIW process to reveal the viscosity, velocity, and shear stress profiles is highly desired.

Although CFD methods have been adopted in studying the shear-thinning property of non-Newtonian fluids, however, the studies of non-Newtonian fluid theories are still being established, and it is far more challenging than that of the typical Newtonian fluid. Consequently, a research gap exists in using numerical simulation to investigate the DIW process. Herein, aiming to investigate and quantify the effects of syringe-nozzle geometry on the ink flow behavior in the DIW process, we employed three types of syringe-nozzle geometries for CFD simulation studies. The conical nozzle, cylindrical nozzle, and conical body with cylindrical tip nozzle are the most widely used configurations in DIW. The shear-dependent viscosity variation is verified by the CFD simulation and validated by the experimental test. The shear stress and velocity profiles of each syringe-nozzle model are reported and compared. The shear stress and velocity values, distributions, and variations are critical factors toward a stable and high-resolution DIW process.

Furthermore, we found high consistency in our simulation results with both our and other published experimental studies. To the best of our knowledge, our work is the first comprehensive physics-based simulation study dedicated to demonstrating the shear-thinning ink property and comparing the ink flow behavior in different syringe-nozzle geometries. Moreover, we examined the strengths and weaknesses of each syringe-nozzle model. On the basis of these comparative studies, we proposed an improved solution for high-resolution and high-shape fidelity DIW. In parallel with the simulations, experimental verifications of the shear-thinning rheological properties and the printing qualities using different syringe-nozzle geometries are carried out. This research contributes to the development of the direct ink writing process by providing an in-depth investigation of the shear-thinning flow behavior in DIW and a guideline for future improvement for high-quality printing. The rest of this article is organized as follows: the CFD methods and experimental procedures are explained in Sec. 2. Section 3 presents the simulation case studies using different syringe-nozzle models. The ink viscosity, shear stress, and velocity profiles are revealed, and a comparison of each model is provided to illustrate the strengths and limitations. The experimental validations including ink rheological characterizations and printing tests using different syringe-nozzle geometries are presented. An in-depth and comprehensive discussion is provided in Sec. 4.

The yield stress τ_yield is the shear stress at zero shear rate, k is the consistency coefficient that measures the average viscosity, and n is the flow index, which is a measure of deviation from a Newtonian fluid. When 0 < n < 1, the ink is a non-Newtonian fluid. The more the significant shear-thinning of ink, the closer n is to 0. In addition, when n = 1, the viscosity value becomes constant, denoting a typical Newtonian fluid. In practice, it is found that when n ≥ 0.8, ink behaves as a Newtonian fluid [19]. The parameters used in the simulation are collected from several references and are summarized in Table 1.

We adopted the standard $K−ε$ model in ansys for all the simulation studies.

Three types of syringe-nozzle geometries are incorporated for the case studies, namely, conical nozzle, cylindrical nozzle, and conical body with cylindrical tip nozzle, with nozzle orifice diameter of 0.6 mm for all case studies. These three types are the most common geometries used in the literature. Due to the axisymmetric structure of syringes and nozzles, many literature utilized the 2D model [21,34–36] or the half-plane model [25] for faster convergence and reduced calculation complexity. However, we specifically constructed 3D computer-aided design (CAD) for demonstrating the dynamic viscosity change, ink flow velocity, and wall shear stress. The 3D structured mesh is assembled using the sweep-axisymmetric method with hexahedra grids in the center. This algorithm guarantees that the simulation result does not deviate at the nozzle region where the cross section is small compared to the syringe body. We also performed convergence studies to verify the grid independence of the results. Based on the convergence study, the element (grid) size is set at 0.05 mm, and the skewness is of 0.39. The total elements of the structured mesh are 1,254,084 with 5,045,659 nodes for the cylindrical nozzle model.

The constant pressure (5 psi) condition is imposed at the syringe inlet for all the models. The normal and axial forces at the outlet are set at zero for all the case studies. The wall is subjected to no-slip boundary conditions.

The rheological characterization is carried out to verify the shear-thinning properties of the ink flow used in the simulation. We prepare three printable inks for the rheological tests, including the silica aerogel-based ink, the silicon dioxide-based ink, and the alumina-based ink. The silicon dioxide nanoparticle is purchased from SkySpring Inc., and the alumina powder is purchased from Inframat Inc. The ink is prepared by mixing the particles (70 wt%) in deionized water followed by homogeneously mixing for 5 h, at 600 rpm on a magnetic stirring hot plate (Thermo Scientific, Waltham, MA) [37]. Cellulose nanopowder is added to the ink in 1.1 wt% as a viscosity modifier and mixed for 6 h. All of the chemicals are as-purchased without further purification.

The rheological behavior of the printable ink is characterized at room temperature using a rotational rheometer (Anton Paar MCR 72, Ashland, VA). The measurement setup is equipped with a parallel plate–plate geometry. The diameter of the plates is 25 mm, with a 1 mm gap in between. The apparent viscosity measurement is performed at a shear rate that ranges from 10⁻² to 10² s⁻¹.

A customized pneumatic pressure-based DIW printer is used for the experimental tests. An air compressor (P1.5IU, Ingersoll Rand) is used to provide pneumatic pressure of 5 psi regulated by a Performus V controller (Nordson EFD, Westlake, OH). Inks are loaded into 5 cc luer-lock syringe (Nordson EFD, Westlake, OH). Helix-locked 800-µm nozzles, including conical nozzle, cylindrical nozzle (Nordson EFD, Westlake, OH), and conical body with cylindrical tip nozzle (Fisnar, Germantown, WI), are attached to the syringe. The extruder path is generated using the open-source slic3r software (slic3r.org). Samples are printed onto an aluminum substrate at room temperature. The as-printed samples are dried at ambient pressure and at room temperature for 12 h.

We divide the syringe-nozzle model into three regions for analyzing the simulation results. Region 1 denotes the flow inlet area. Region 2 is the syringe body with a simple cylindrical geometry. Region 3 describes the nozzle outlet area. As outlined in Section 1, during extrusion printing, pressure is exerted at the inlet of the syringe (in region 1), forcing the ink to flow along the longitudinal direction (z-axis in Fig. 2) of the syringe (in region 2). Finally, ink is deposited through the nozzle outlet (region 3). Figure 2 shows the simulated result of ink viscosity profiles. In general, high viscosity is observed in the inlet area, where the shear deformation is small. A more significant shear deformation occurs at the position below the inlet area, leading to a sharply decreased viscosity (as shown at the interface of regions 1 and 2). At region 2, the ink viscosity reaches an equilibrium state as indicated by the uniform color, denoting the flow shear rate $(γ˙)$ remains a constant. As the syringe-nozzle geometry necks down at the outlet region, the viscosity again decreases sharply (refer to the interface of regions 2 and 3). The lowest viscosity appears at the nozzle tip, where the shear rate is the largest, and the radial dimension is the smallest. The viscosity contours in Fig. 2 reveal the “shear-thinning” property of the non-Newtonian ink. When the ink flow experiences the shear deformation, it becomes thinner than that in the quiescent state. The term “apparent viscosity (Pa s)” is often used to describe the shear-dependent viscosity or used as a measure of resistance to the shear deformation of such non-Newtonian fluids.

Understanding the flow shear stress during extrusion is vital in many engineering fields. Low shear stress is favorable in biomedical applications for avoiding shear-induced degradation. In contrast, high shear stress could be beneficial to enhance the sample properties after extrusion printing. Figure 3 reveals the shear stress changes during extrusion. The shear stress is low in the inlet area, then quickly increases, and remains stable at the syringe body regions. Finally, at the narrower nozzle region, the shear stress again quickly increases and reaches the maximum value at the nozzle tip. From the color contour, we can see that the shear stress can be stabilized within the cylindrical tips (Figs. 3(b) and 3(f)), and this trend is consistent with the viscosity results. Figure 3(g) presents the curves of shear stress values along the longitudinal direction. It is worth noting that the cylindrical nozzle has more prominent stress at the nozzle tip than the other cases. This finding provides valuable guidance for the syringe-nozzle design. The high shear stress could facilitate the alignment of filler materials such as rod-like fibers [38,39] and 2D materials like graphene [26,40] during the extrusion 3D printing. By aligning the filler materials in the axial direction, it enhances properties such as mechanical strength, electrical conductivity, and thermal insulation [41]. Conversely, low shear stress could be critical to bioprinting to ensure the cell viability, which is a measure of live and healthy cells within the ink [36,42,43]. Overly high stress could lead to severely degraded cell viability. In this case, the conical nozzle or the conical nozzle with the cylindrical tip nozzle would be ideal tools.

The velocity contours at the cross-sectional plane (x-z plane) are shown in Figs. 4(a)–4(i), and the velocity profiles are summarized in Fig. 4(j). The velocity at the outlet region is instructive for calibrating the extrusion printing parameters (e.g., printing speed). We noticed that velocity could be stabilized within the cylindrical tips (Figs. 4(c), 4(i), and 5), and this trend is consistent with the viscosity and shear stress results. The stabilized ink velocity at the nozzle is critical to ensure uniform extrusion and the high fidelity of the 3D-printed specimen. Unstable velocity variations are exhibited at the interface of regions 2 and 3 in the conical body with a cylindrical tip nozzle and a conical nozzle (shown in Figs. 4(b) and 4(e)). These variations lead to locally sudden-increased velocity and asymmetric velocity distribution. The reason could be linked to the acute necking of the model geometry. Considering this abnormal variation, we further examined the shear stress values within the necking region and concluded that a “dead zone” may exist as a result of the acute necking geometry. The ink flow is aggregated within the “dead zone,” resulting in a locally increased shear stress and velocity. On the other hand, such local variation is not observed in the cylindrical nozzle (refer to Figs. 3(f), 4(h), and 4(i)). In this case, considerably large shear stress (five to six times greater than in other cases) is presented in the cylindrical nozzle. The large shear stress facilitates the flow and thus eliminates the flow aggregation. A similar result was also reported by Shao et al. [25]. This finding is instructive in designing an optimized model for uniform extrusion (in Sec. 3.5). Our previous work [37] conducted a comparative experimental study using three types of nozzles for extrusion printing. We find that the conical body with a cylindrical tip nozzle has the most uniform extrusion and thus ensures the best fidelity of printed parts. The result is also consistent with the simulation studies as we presented in this work.

The outlet velocity value and contours are summarized in Fig. 5. From the contours, we observe that the center velocity has a plateau region and gradually decreases to the periphery. The decreasing velocity from the inner area to the outer is because the models are subjected to the no-slip boundary condition, denoting that the periphery velocity shall be zero. In this regard, a gradually decreased velocity from the core to the periphery would be preferable. Furthermore, a greater percentage of the plateau region is favorable for extruding the cylindrical-shaped filament, ensuring a high printing resolution. The simulated plateau velocity has the value of 5 mm s⁻¹ to 10 mm s⁻¹, which is consistent with our previous experimental studies [8] and other published works [27,34,44–46].

The aforementioned case studies provided a valuable comparison of each syringe-nozzle model. Based on the simulation case studies, it can be summarized that with a cylindrical tip at the nozzle outlet, the ink velocity and shear stress can be stabilized, providing a favorable condition for high-resolution printing. On the other hand, the sharp geometry change could lead to locally fluctuated flow and should be addressed. We summarized the detailed pros and cons in Table 2.

In light of the comparison, we proposed a preferable solution to integrate the merits from each model and optimize the geometry to eliminate the unsteady flow. Specifically, the optimization targets are as follows: (1) eliminate the “dead zone” area to create smooth ink flow at the syringe necking region and (2) maintain the steady velocity at the outlet and reduce the shear stress at the tip.

The proposed model extended the basic geometry of the conical body with a cylindrical tip nozzle. The original sharp edges at the syringe necking region were changed to curved and smoothly converged edges (Fig. 6(a)). The tip length was elongated from 1 mm to 2.5 mm. The simulated velocity results are shown in Fig. 6. The contour (Figs. 6(a)–6(c)) and the profile (Figs. 6(d) and 6(e)) revealed no unsteady flows throughout the entire model. The curved edge facilitates the ink flow by mitigating the stress concentration, and the velocity can gradually increase at the necking region and nozzle region. In addition, steady velocity holds at the nozzle outlet. Together with the shear stress results in (Figs. 6(f) and 6(g)), it shows the “dead zone” area was eliminated in the proposed model. The maximum shear stress was much lower than the cylindrical nozzle and marginally higher than the other two cases as summarized in Fig. 6(i). This improved solution could be further applied to multinozzle and multimaterial 3D printing, where the stable extrusion behavior is at the heart of these research areas [47].

We perform the rheology measurements for the ceramic inks with silica aerogel, alumina, and silica nanoparticle. The apparent viscosity (Pa s) is plotted against the shear rate (s⁻¹) in Fig. 7. The curves show that all the inks exhibit the shear-thinning behavior under shear stress, as evidenced by a rapid decrease in the apparent viscosity with respect to the increasing shear rate from 10⁻² to 10² s⁻¹.

Physically, the apparent viscosity during the extrusion is very low, and ink exhibits typical fluid-like behavior. The viscosity is thinned to be low enough for the continuous filament extrusion. Upon exiting the nozzle, as shear deformation vanishes, the apparent viscosity increases and instantly recovers the elastic (solid-like) behavior. Therefore, the ink can form and maintain the self-supported 3D shape without wetting or spreading during the direct ink writing process.

To validate the findings from simulation studies, we customized a pneumatic pressure-based direct ink writing printer (Fig. 8(a)) and printed the silica ink using three types of nozzles, namely, the conical body with the cylindrical tip nozzle (Fig. 8(b)), the conical nozzle (Fig. 8(c)), and the cylindrical nozzle (Fig. 8(d)). The geometry of each nozzle type corresponds to the CAD model in simulation studies. From the printing results, the most intuitive distinction between each sample is the printing quality and shape fidelity. (i) Conical body with the cylindrical tip nozzle: From the simulation results, we learn that this nozzle has stable extrusion velocity owing to the cylindrical tip geometry, and large plateau velocity distribution at the nozzle outlet, which are preferable conditions for high quality and uniform extrusion. As can be seen in Figs. 8(f) and 8(j), the printed samples have high quality in terms of layer uniformity, shape fidelity, and resolution. (ii) Conical nozzle: Owing to the gradually narrowed cone shape, this nozzle has the lowest shear stress in ink extrusion. However, the velocity distribution at the nozzle outlet (Fig. 5) is a bell-shaped curve. The large velocity difference between the nozzle center (0.01 m s⁻¹) and the peripheries (0.004 m s⁻¹) may lead to nonuniform extrusion. This can be verified from the experiment tests shown in Figs. 8(g) and 8(k), where the printed samples have poor surface quality. Geometrical details such as twisted contours and spherical curves cannot be observed. (iii) Cylindrical nozzle. Similar to the conical body with a cylindrical tip nozzle, the cylindrical shape can stabilize the shear stress and extrusion velocity and leads to uniform extrusion behavior. High surface quality can be found on the printed samples shown in Figs. 8(h) and 8(l).

In practice, the three types of syringe-nozzle models included in our simulation study are all widely used for the 3D DIW process. Pilot DIW works by Lewis et al. [48–51] employ the cylindrical nozzles for DIW. As this field rapidly develops, the low-cost conical nozzles made from plastics are becoming popular [27,52,53]. On the contrary, the conical body with cylindrical nozzles is costly, but it could play a vital role in high-resolution DIW [9,37]. Recently, as more research endeavors are invested, and methodologies delved deeper, the attention has been directed toward the DIW resolution and fidelity. The high-resolution and high-shape fidelity 3D printing assembles the feedstock materials into continuous and conformal 3D structures. Consequently, it ensures a more promising performance, such as electron transfer, energy storage, thermal conductivity, and optical transmittance of the printed specimen [1]. For example, Chortos et al. [54] present the multinozzle DIW for soft and flexible dielectric elastomer actuators, possessing significant high-energy densities and fast actuation rates. The printing fidelity is particularly important since any defects could cause thinning dielectric segments that reduce the breakdown field when bending or folding the devices. In this regard, the outcomes of this work could serve as an informative guide for experimental studies in the following two aspects.

1. Understand the ink flow behavior during DIW and select the appropriate nozzle geometry accordingly. For example, from our comparative simulations, we learned that the conical nozzle and conical body with a cylindrical tip nozzle have low shear stress. Shear stress is a major contributing factor to cell viability during 3D bioprinting. Therefore, it is making the conical nozzles favorable candidates to maximize cell viability. On the other hand, high shear stress could present on the cylindrical nozzle, generating a favorable condition for aligning the 2D filler materials (e.g., fiber) in ink. The orientation of fibers within a composite is a significant factor that defines the ultimate strength. When subjected to high shear stress, the fibers can overcome the random orbital rotation and align in the flow direction [55,56]. This mechanism also works with other 2D materials such as graphene, boron nitride, and MXenes. In addition, the nozzles with cylindrical tips have stable extrusion shear stress and velocity at the outlet. Thus, it can be used to precisely control the printing resolution and shape fidelity, as shown in the experimental verification. Moreover, the velocity at the nozzle outlet can be referenced for defining the printing speed, an essential value that shall be specified while slicing the CAD model, saving efforts for the time-consuming trial-and-error process.

2. Identify the potential “dead zone” area in the model and propose the improved geometry towards stable and high-resolution printing. As shown in Sec. 3.3, the dead zone appears at the syringe necking region. However, it could be eliminated by changing the shape corner to curved corners in the proposed model. In fact, avoiding the dead zone has more valuable meaning for physical experiments, primarily because most of the fluids used for DIW are multiphase fluids [1], such as particle suspensions [8,27,37,39,51], fiber-induced suspensions [26,40,57], oil–water emulsions [58,59], and foamed inks [8,37,60]. The different phase materials may not be uniformly dispersed within the ink, or the uniformity may degrade during extrusion. In such cases, any nonuniformity could lead to unsteady flow, resulting in clogging issues at the “dead zone” area and inconsistent ink dispensing as the DIW process progresses. The stable extrusion is even more pivotal when deploying DIW for large-scale 3D printing. A larger printing scale increases the likelihood of shape deviation from the predefined model caused by unstable or nonuniform extrusion. Therefore, the proposed improved solution is a promising candidate for stable and high-resolution DIW.

In summary, we conducted the CFD simulation of shear-thinning inks in the DIW process. The results reveal the ink flow behavior throughout the syringe-nozzle model. Comparative case studies demonstrate the ink shear stress and flow velocity in three different syringe-nozzle geometries and found high consistency of our simulation results to the physical experiments. The results offer readers a foundation for understanding the DIW process and guidance on selecting the nozzle types, economizing the efforts of preliminary trial-and-error experiments. Moreover, we identified the “dead zone” that could cause flow clogging. On the basis of the findings, we proposed the improved syringe-nozzle geometry design for smooth and uniform ink flow that could lead to high-resolution DIW and enhanced performance of the printed samples. To validate the simulation results, we use the corresponding DIW configuration to print 3D models and compare the printing quality. Consistent with the simulation analysis, we find the cylindrical nozzle and the conical body with a cylindrical tip nozzle have high printing quality and shape fidelity. To the best of our knowledge, our work is the first comprehensive physics-based simulation study of the direct ink writing process. We believe the simulation results discussed here can motivate further developments of the direct ink writing process.

The authors would like to gratefully acknowledge the supports from the National Science Foundation (NSF) through the award CMMI-1846863, FM-2134409, and the U.S. Department of Energy (DOE), Office of Energy Efficiency and Renewable Energy (EERE) through the award DEEE-0008675.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

k =

consistency coefficient (Pa s)

n =

flow index

K =

kinetic energy (kg m² s⁻²)

$γ˙$ =

shear rate (s⁻¹)

$ε$ =

dissipation rate (m² s⁻³)

µ =

apparent viscosity (Pa s)

$v$ =

velocity vector (m s⁻¹)

ρ =

density (kg m⁻³)

τ =

shear stress (Pa)

τ_yield =

yield shear stress (Pa)