Process fundamentals and quality investigation in extrusion 3D printing of shear thinning materials: extrusion process based on Nishihara model

Due to the complexity of extrusion 3D printing objects and lack of universal parameters, there are few classified statistics or quantitative description for the “extrudability” of material with shear thinning properties. In order to systematically study the rheological properties of extrudable materials and establish a unified criterion or quality evaluation standard from the basic process theory, a Nishihara rheological constitutive model based on visco-elastic and visco-plastic characteristics of materials is explored to predict and investigate the extrusion performance of four typical dispersed-continuous phase binary mixture. The creep model result shows that, expect for the gelatin network-water system with almost complete elastic behavior, the material with shear thinning property based on Nishihara model is in good agreement with the experimental results. In the extrusion process, the shear stress is related to the advancing speed, material viscosity, and nozzle size. The effects of advancing speed and nozzle size on shear stress show antagonistic characteristics in a certain range; that is, the velocity gradient is the dominant factor at lower extrusion speed, and the dynamic viscosity is the dominant factor at higher extrusion speed. In terms of extrusion properties, the material system with smaller yield strain/stress has the least obvious extrusion delay characteristics, and is easier to extrude under the condition of the same material strength.


Introduction
Extrusion printing is the most common printing technology at present, which is very similar to the traditional FDM 3D printing in extrusion and molding process. This kind of printing technology needs to put ink material into the container, and then use pneumatic pressure or mechanical load as the driving force to extrude it to the predetermined position on the manufacturing platform through the nozzle. Extrusion printing is suitable for a relatively wide range of ink materials with corresponding process under different parameters based on rheology analysis of ink properties [1]. It has high flexibility and can easily switch between various materials to realize printing of different 3D configurations, which has special significance for simulating artificial tissues or organs [2][3][4][5][6]. Although the precision of extrusion printing is relatively low (the resolution is difficult to reach 20 μm), due to the simple device and low cost, it is still the most suitable 3D printing method in printing tumor models, prostheses, medical implants, and other tissues with low demand for structural accuracy [7][8][9][10][11][12]. Therefore, by systematically studying the process basis of extrusion 3D printing and realizing the programmed control of printing accuracy, it will greatly broaden the application field and achieve more high-end and accurate application requirements.
Due to the high shear stress, it is generally believed that the cell survival rate of extrusion 3D printing (the cell survival rate is about 60%) is often lower than that of other printing methods such as ink-jet printing (85%), stereo lithography appearance (80%), and laser induced forward transfer (95%) [13]. In order to solve this problem, many ink materials nowadays have the property of shear thinning, which makes the apparent viscosity of ink decreases with the applied high shear stress, so as to reduce the damage of shear stress to cells. Therefore, in the design of extrudable bio ink which can be used at room temperature, the mixing mode of dispersed-continuous phase tends to be used to ensure both of the mechanical strength and biocompatibility of printing matrix.
The materials used for extrusion 3D printing mainly include organic materials and inorganic materials. Organic polymers and composites usually exist in the form of filamentous thermoplastic polymers, active monomers, resins, and powders. They not only have their own rheological properties to meet the needs of extrusion printing process but also can be used as the dispersed phase or continuous phase in the two-phase mixed system to form composite ink according to the solidity. In addition, new printing properties of the system are also developed by the two-phase interaction. For example, Serra et al. mixed polylactic acid with bioactive CaP glass to print porous tissue engineering scaffolds with biocompatibility [14]. Tekinalp et al. studied the mechanical bearing capacity of the printed parts composed of carbon fiber and ABS resin [15]. Postiglione et al. prepared polymer nanocomposites based on PLA and MWCNTs, and developed a low-cost liquid deposition molding technology, which can print 3D conductive microstructures with arbitrary shapes [16]. In contrast, inorganic nanoceramics are mainly printed by powder or ink. The powders are sintered by laser or bonded together by auxiliary adhesive. If the adhesive has the property of light curing, stereolithography can be used for direct 3D printing of ceramic materials [17]. On the other hand, in order to realize the curing process of ink-jet or extrusion, post-treatment such as high-temperature sintering is used for ceramic particle suspension after printing [18]. By compounding with flow continuous phase, even ceramic particles that are difficult to flow under normal circumstances can be extruded for 3D printing. For example, in order to increase the bending strength of materials, Li et al. added calcium sulfate and dextrin into porous alumina ceramics to enhance the connection between different phases [19]. Maurath also pointed out that by improving the rheological properties and uniformity of ceramic suspension, the printing process and sintering process can be optimized to print the honeycomb structure with high strength, which had no cracks and better dimensional stability [20].
In terms of the model research of extrusion stage in the 3D printing process, the creep model commonly used at present has undergone the following three stages of evolution. In the initial stage of development, Tavenas and Grutin et al. used empirical/semi empirical models to express the linear relationship between strain or strain rate and time by semi logarithmic or double logarithmic coordinate systems [21,22]. After that, Feda et al. started from the basic mechanical concepts and mathematical theories to simulate the creep characteristics of soil by using the combination of model elements (Hooke elastic body H, Newton viscous body H, Saint Venant plastic body H) [23]. On this basis, Ollier, Güneyisi and Ou further applied the composite rheological model formed by the combination of these models to the polymer and clay materials with similar special properties. In addition, the viscoelastic plastic model is developed on the basis of the classical plastic theory [24][25][26]. Pronk et al. believe that the plastic potential is consistent with the creep potential, and establish a comprehensive model to describe the rheological properties of materials [27]. Duty, Ibrulj and Hossain respectively carried out rheological analysis on extruded 3D printing materials with nonlinear mechanical properties from the perspectives of theory, modeling and experimental verification, and made application prospects for the shear thinning characteristics of the materials [28][29][30].
Due to the complexity of printing objects, multi material strategy has become the important way for the development of extrusion printing. Because of the different composition mechanism of various types of inks, there are few classified statistics and general description of properties. In addition, the materials of each printing model are highly targeted, which leads to the lack of universal parameters. In terms of the extrudability of materials, it is difficult to describe the shear stress distribution in the nozzle which is directly related to the cell survival rate. Due to the different rheological properties of materials, the extrusion state and difficulty mostly depend on empirical judgment, and there are no quantitative fluid mechanics calculation results currently. Therefore, it is urgent to break through the limitation of qualitative description and establish a unified criterion or quality evaluation standard from the basic process theory.

The technological process of extrusion 3D printing: extrusion process
The technological process of extrusion 3D printing shown in Fig. 1 includes the process of extruding materials from the container to form fibers and the process of stacking fibers on the substrate according to the predetermined trajectory. The general extrusion 3D printing uses heat source to pretreat the material, and cooling after printing to complete the phase transformation and curing process. In contrast, the extrusion 3D printing at a room temperature has special rheological properties and a certain degree of shape retention ability. In the process of extruding and flowing out from the nozzle, the material has experienced high to low shear, that is, high shear through the finer nozzle and low shear under gravity after extrusion. Therefore, extrusion 3D printing requires that the target material has the mechanical stimulation response characteristics of shear thinning and the effect of fast and reversible modulus on shear stress, so as to ensure that the material can be extruded out of the nozzle to become fibers smoothly and keep the configuration from collapsing in the process of stacking and deposition.
In the whole process, the key extrusion process depends on rheology to describe and evaluate the flow and deformation behavior of materials, especially in the transient region between solid and fluid [31]. Extrusion 3D printing material system is a viscoelastic body between ideal fluid and ideal solid in rheology. When a kind of viscoelastic material is subjected to shear, it shows both viscous and elastic properties; that is, it can flow and deform under the action of external force at the same time [32]. Elastic behavior means that the deformation characteristics can be restored after the load is released, which just like a spring that obeys Hooke's law. Generally speaking, the process of restoring deformation after the material is stressed requires certain recovery energy. During the deformation, this part of energy is stored in the deformed material, which is generated by expanding the inner superlattice structure. When the applied stress does not exceed the structural strength of the material, the stored energy will be used as the driving force to restore the structure to its original state after the load is released. On the other hand, the viscous behavior comes from the friction between the components in the flowing fluid, which just like a dashpot that obeys the Newton's law of viscosity. With the heat generated by friction, the deformation energy will be converted into heat consumption, also known as energy loss.
The viscoelasticity of materials can be characterized by the storage modulus G′ and loss modulus G″. The storage modulus is the response to the elasticity and mechanical strength of the material system, which represents the internal stress during the recovery process and the energy stored by the elastic deformation, while the loss modulus corresponds to the viscosity characteristics of the system, which represents the energy dissipated in the form of heat due to viscous deformation, reflecting the viscous nature of the material. The ratio of loss modulus to storage modulus (G″/G′ = tanθ) is called loss tangent or damping factor, which can characterize the mechanical loss of materials. The smaller the loss tangent is, the closer the material is to the ideal elastic material [33].
In a certain range of low shear and small deformation, there is a linear relationship between stress and strain (also between stress and strain rate), which is called linear viscoelastic region (LVE). The structure change of the material in this region is reversible elastic deformation; that is, after the load is removed, the structure of the material can be generally restored to the original state [34]. When the applied strain is beyond the LVE, the structure of the material will be irreversibly changed with relatively complex response of stress and strain. This critical point is called yield point (or yield stress), which refers to the shear stress at the boundary of LVE [35]. The yield stress is mainly due to the 3D network structure with strength formed by the interaction between molecules in the material. When the shear stress exceeds the yield point, the network structure will be destroyed. For different 3D printing materials, the structural strength can be characterized by synthesizing the G′ value and yield stress value of the LVE [36].

Material and methods
Materials based on extrusion 3D printing generally need to have thixotropic properties; that is, the strength and viscosity of materials are sensitive to shear force, showing the characteristics of "shear thinning." When the external loading is removed, the material can quickly restore the modulus and maintain the solid shape. In order to ensure both of the fluidity and formability at the same time, the two-phase mixed system is selected as the basic composition of the printing material (inorganic-inorganic, organic-inorganic and inorganic-organic systems are selected as the basic types of dispersed-continuous phase). The materials and characteristics used in this work are as follows ( Fig. 2a-d): Laponite: Magnesium lithium silicate is a kind of synthetic inorganic flake silicate additive, with the diameter of about 30 nm and thickness of about 1 nm. When the film is evenly dispersed in water, the flakes are peeled off with negatively charged surface and positively charged edge. If the content is higher than the critical concentration, a "card house" structure is formed due to strong electrostatic attraction between the layers of the clay, which is connected to the polymer chain by hydrogen bonding and coordination as shown in Fig. 2a [37]. The mechanical strength of laponite hydrogel is relatively high with regular lamellar network, which can effectively reduce the stress concentration. Moreover, the molecular chain with high shrinkage can also slow down the stretching during the elongation process, thus showing strong tensile properties. It is noteworthy that the non-covalent bonds are formed again through the diffusion of polymer chains in the damaged laponite hydrogel, which means that the gel has good self repairing property [38]. Laponite has good water swelling property and strong adsorption capacity. It is widely used in food, medicine and cosmetics as suspension, paste thixotropic agent, emulsion and ink stabilizer, thickener, and so on [39].
Alumina NPs: Nano-alumina is a kind of ultrafine particles with the size of 1 ~ 100 nm, which has strong volume effect, quantum size effect, surface effect, and macroscopic quantum tunneling effect. Due to its high mechanical strength and biological inertia, nano-alumina particles are often used as 3D printing materials for artificial bones or scaffolds [40]. In general, the agglomeration of alumina nanoparticles depends on the surface charge shown in Fig. 2b, so even if they are fully dispersed in aqueous solution, they can only get suspension and cannot be self-assembled. In order to carry out the molding process independent of high temperature sintering, the dispersed phase with the ability of fixing nanoparticles must be used to make it into a semi fluid state for extrusion. In addition to the high-temperature resistance and friction resistance characteristics of ordinary ceramic materials, alumina NPs also have porosity, high dispersion, and high activity, which can be used in the coating of high-quality inkjet printing, toughening filler of mobile phone back plate, chemical catalyst carrier, and other fields [41].
Poloxamer (pluronic): Poloxamer is a block copolymer composed of hydrophilic polyoxyethylene and hydrophobic polyoxypropylene ether. Its commercial name is pluronic, and it is generally used as a polymer nonionic surfactant. Because of its amphiphilic properties, when the content in aqueous solution is higher than the critical micelle concentration (CMC), these copolymers will automatically assemble into micelles. The center of the micelle is hydrophobic polyoxyethylene, and the outer layer is hydrophilic polyoxyethylene chain (as shown in Fig. 2c) [42]. Pluronic is a temperature sensitive material, which means that all kinds of products in liquid state at room temperature do not have molding ability. Therefore, it is necessary to anchor the dispersed phase particles in its hydrophilic block to provide mechanical strength and stable molding. The amphiphilic property of pluronic can be used to increase the water solubility of hydrophobic oily substances or increase the miscibility of different hydrophobic substances, and it is commonly used in industrial products, cosmetics and pharmaceuticals. In the biomedical field, it is also used to evaluate the delivery effect of many drugs, and can also be used in cell culture media by using its buffering effect [43].
Gelatin: Gelatin is a kind of white or yellowish translucent powder formed by partial degradation of collagen in animal skin, bone, sarcolemma, and other connective tissues, which belongs to macromolecular hydrophilic colloid. The gelatin hydrogel mainly depends on intramolecular hydrogen bonds and hydrated hydrogen bonds to maintain stability. Gelatin hydrogel is sensitive to environmental stimuli, such as temperature, pH, and ionic ambience, which results in phase transition behavior shown in Fig. 2d [44]. It can be soluble in hot water and has good biocompatibility, which makes it widely used in biomedical materials. Gelatin hydrogel has a solid morphology with good elasticity and poor plasticity. This is caused by partial gelatin molecules returning to collagen structure in solution. In order to enhance the mobility of the hydrogel, gelatin can be made into nanoparticles to weaken its intramolecular hydrogen bonding interactions. As the hydrolysate of natural protein, the biodegradation of gelatin can play a significant role in the field of medicine. After gelatin is prepared into nanosized solid particles, it can be used as an excellent drug and gene delivery carrier for controlled or targeted release to enhance the cellular uptake ability of some materials and avoid toxic/side effects, improving the therapeutic effect of drugs [45].
In this work, the DHR-1 type rotary rheometer produced by TA Instrument Company is used to analyze the rheological properties of different material systems. The instrument adopts magnetic suspension thrust bearing, porous carbon radial bearing, and optical coding displacement sensor with axial oscillation force range of 0.1-50 N, oscillation displacement range of 1-100 μm, axial displacement resolution of 20 nm, frequency range of 6 × 10 −5 -100 rad/s (10 −5 -16 Hz), step strain response time of 10 ms, step rate response time of 15 ms, normal force range of 0.005-50 N, and normal force resolution of 0.001 N.
In order to reveal the role of materials in the mixed printing ink, the rheological properties of several typical dispersed-continuous binary systems (laponite-water, alumina NP-pluronic, gelatin network-water, and gelatin NP-water) with different ratios of parameters are analyzed comprehensively, including frequency scanning to measure the modulus and composite viscosity, amplitude scanning to measure the yield stress/strain, flow testing to measure the shear thinning behavior and hysteresis behavior, and step testing to evaluate the creep recovery and stress relaxation behavior, which is used to describe the extrusion performance of various materials.

Basic rheological models of materials
The extrusion process of 3D printing material includes the process that the material reaches the yield shear stress/strain through the syringe thrust and the steady-state extrusion process that the material enters the plastic deformation stage. Therefore, the compressibility and non-Newtonian fluid characteristics of materials in the whole process need to be described by a corresponding constitutive model. In rheology, all material models can be composed of three basic elements: elastic element (H), viscous element (N), and plastic element (S). In the practical engineering application of extrusion 3D printing, the properties of ink composed of binary mixture are not single viscous, elastic, or plastic. In order to accurately describe the rheological properties of extrusion ink, it is necessary to combine the three elements. Basically, these three elements can be combined to form four initial time-dependent models, namely, viscosity, viscoelasticity, viscoplasticity, and viscoelastic plasticity, which are called basic rheological models. Taking one to four of the basic rheological models in series, 15 rheological models (including 4 basic rheological models and 11 composite rheological models) can be obtained. That is to say, there are 15 rheological behavior states that can be described by the theoretical rheological mechanics model [46]. Several typical composition modes are shown in Table 1 (where "-" means series connection and "║" means parallel connection).

Constitutive equation analysis of Nishihara model
As mentioned above, Nishihara model is composed of an elastic element, a Kelvin body, and ideal viscoplastic body in series, which can best reflect the viscoelastic-viscoplastic characteristics of extrudable materials. The composition is shown in Fig. 3.
The method of replacing linear elements with nonlinear elements (mostly viscous elements) can reflect the nonlinear characteristics of rheological behavior of binary mixtures (non-Newtonian fluids). However, each single component in the model element method is one-dimensional. When deriving the three-dimensional constitutive relationship, the method similar to the one-dimensional Hooke's law extended to three-dimensional in elastic theory is adopted at present, and the following assumptions are required: 1) The volume deformation of the material is completed at the moment of stress, and it is elastic deformation, which does not change with time.
2) The creep is only caused by the deviatoric stress tensor, but the material does not creep under the action of the spherical stress tensor (in the creep state, the object only changes in shape, but not in volume). 3) In the rheological process, Poisson's ratio does not change with time. According to the mechanical model of the elastomer, it can be expressed as follows: Under the condition of constant load , that is, = 0 , the deformation develops at a constant rate, and the creep model of the elastic element can be obtained as follows: A Kelvin body is composed of an elastic element and a viscous element in parallel, and its constitutive equation is: When the constant load is unloaded at t = t 1 , that is, = 0 , and = 1 at this time. However, with the increase of time, the strain gradually decreases and tends to zero when the t increases to a certain value. Therefore, the creep equation is as follows: An ideal viscoplastic element consists of a plastic element and a viscous element in parallel. Under constant load , the deformation increases uniformly with the increase of time, and the deformation does not recover after unloading at t = t 1 . Therefore, its constitutive equation is different before and after reaching the yield stress: According to the composition of each element, the constitutive equation of the above model can be expressed as follows: When the applied stress does not reach the yield stress: When the applied stress reaches or exceeds the yield stress: where is the tress of model, is the strain of model, s is the yield stress of viscoplastic element, k 1 is the elastic coefficient of elastic element, k 2 is the elastic coefficient of viscoelastic element, 1 is the viscous coefficient of viscoelastic element, 2 is the viscous coefficient of viscoplastic element, ̇ and ̈ are the time derivative of model stress, and ̇ and ̈ are the time derivative of model strain. According to the above equations, the creep equation is obtained as follows: where 0 is the transient stress and t is the loading and unloading time.
The creep loading and unloading curves are shown in Fig. 4.
It can be seen from the figure that when the stress level is low (not reaching the yield stress), the material deforms rapidly at the beginning ( t is very small), and gradually tends into stable creep (the strain is constant when t tends to infinity). When the stress level is equal to or exceeds a critical value ( s ), the material will gradually transform into unstable creep (the strain increases with the increase of t when it enters the plastic deformation region). At the moment when the stress is removed, the material will rebound instantaneously, and then, the deformation will gradually recover with the increase of time (related to the yield strength of the material, the deformation will not recover completely). The model has the properties of instantaneous deformation, decay creep, steady creep, elastic aftereffect, and relaxation, which can well and universally reflect the state of most materials with both viscoelasticity and viscoplasticity.

Acquisition and significance of parameter variables
Generally speaking, the state of typical Newtonian fluid in nozzle only needs to consider its density, viscosity and bulk modulus, but as a non-Newtonian fluid with shear thinning characteristics, it has compressibility and viscoelasto-plastic properties simultaneously. In the simulation based on ANSYS-Fluent module, the yield stress/strain of the material can be obtained according to the minimum shear strength required by steady-state extrusion, and the accuracy of the model can be verified by comparing with the amplitude scanning and step creep results in rheological test. In addition, we can also obtain the wall shear strength distribution of the material during stable extrusion to consider the influence of nozzle morphology, advance speed, and other parameters on the extrudability of the material. The investigation objectives, characterization results, and verification methods are shown in Table. 2.
It should be noted that the steady-state simulation analysis can be used to reduce the calculation difficulty according to the parameters of the current working condition in the stable extrusion state, and the results are equivalent to the transient calculation results in the unsteady process (that is, the extrudable properties are only considered without the delay or deformation). However, due to the compressibility of materials, there is a time delay effect before reaching the yield stress/strain. Therefore, it needs to determine the response speed of materials to shear rate by flow test.

Extrusion process -results and verification of elastic deformation stage
In the extrusion nozzle, the non-Newtonian compressible fluid first experiences the process of gradually applying pressure. In this process, the material only has elastic deformation, and the pressure of the propulsion section increases with the increase of the propulsion distance. There is a certain time difference from the beginning to the yield stress/strain of the material, which leads to the delay of the extrusion material.

Laponite-water system (inorganic-inorganic)
On the premise of ensuring that the material system can maintain the solid shape without external force, the simulation based on the constitutive model of 20-40% mass fraction laponite-water system is carried out with creep analysis under 20 Pa stress. The comparison of creep with time between the numerical model and experimental results is shown in Fig. 5a. The figure shows that the higher laponite mass fraction increases the system strength of the material and decreases the steady-state strain. However, the numerical fluctuation is more significant than the lower mass fraction at the beginning of creep. This is because laponite is a colloidal structure formed by charge induced lamellar interposition. As a dispersive phase, it is the main factor of creep properties in the material system (water does not possess creep characteristics). Due to the plastic deformation capacity, the steadystate strain is slightly larger than the model value, but the experimental results are still basically in accordance with the simulation.

Alumina NP-pluronic L61 system (inorganic-organic)
Pluronic can form polymer chains with different flow states according to the proportion and length of hydrophilic and hydrophobic blocks. In order to ensure the miscibility of the binary system and keep the solid state of the material system without external force, the numerical simulation based on the constitutive model of 20-40% mass fraction alumina NP-pluronic L61 system is carried out with creep analysis under 20 Pa stress. The comparison of creep with time between the numerical model and experimental results is shown in Fig. 5b. It can be seen from the results that the higher mass fraction of alumina NPs obviously reduces the steady-state strain of the whole system, which is because the strength of aluminum NPs is higher than laponite, and the mechanical strength of the system can be enhanced significantly. But different from laponite-water system, the initial creep rises more dramatic under the lower mass fraction. This is due to the high viscosity of pluronic L61, which is the main factor of creep properties in the system as a continuous phase (the alumina NPs are close to rigid spheres with less creep properties). In terms of plasticity, there is no obvious difference between the system of aluminum NP-pluronic L61 and laponite-water; thus, the experimental results are in good agreement with the model as well.

Gelatin-water system (organic network-inorganic)
Generally speaking, gelatin does not dissolve in water. However, it can absorb 5 ~ 10 times of water to expand and soften in soaking condition. When the gelatin is heated, it will dissolve into colloid, and then form a gel phase after cooling below 35 ~ 40 °C. In addition, the packing density of gelatin is smaller than that of laponite and alumina, so it is difficult to mix at the same mass fraction. In this work, the 15-25% mass fraction gelatin-water system is cooled to room temperature from 45 °C and then stabilized to form a stable gel. The numerical simulation based on the constitutive model was carried out with creep analysis under 20 Pa stress. The comparison of creep with time between the numerical model and experimental results is shown in Fig. 5c. Similar to laponite-water system, the creep behavior of gelatin-water system is also provided by dispersed phase. However, the aqueous solution of gelatin is a gel network formed by hydrogen bonding to restore collagen structure. Its structure intensity has good resilience, and is much lower than that of inorganic material dispersed phase system. Therefore, the strain is obvious at the beginning of creep and quickly reaches into a stable state. The similar trend among different mass fractions means that the elastic force is more dominant than viscous force, which leads to a relatively large deviation between the experimental results and the numerical model in the elastic deformation stage for the gelatinwater system. Thus it can be seen that the description of viscoplasticity in Nishihara model cannot be very effective in this kind of organic network-inorganic material system.

Gelatin NP-water system (organic nanoparticles-inorganic)
Gelatin has excellent biocompatibility and biodegradability. By chemically modifying the structure of gelatin, it can not only retain these properties but also construct materials with stronger mechanical properties. By means of condensed phase separation method, 20% NaSO4 aqueous solution and isopropanol are added into gelatin hot solution to promote the coagulation process with the crosslinking and fixation of molecular chain by glutaraldehyde. Therefore, gelatin NPs has a more loose stacking mode than gelatin (the specific gravity is only about one-sixth of that of water), and it is more difficult to mix with water at the same mass fraction. Due to the better water solubility of gel NPs at room temperature, the system with 5-15% mass friction is analyzed by numerical simulation based on the constitutive model under 20 Pa stress. The comparison of creep with time between the numerical model and experimental results is shown in Fig. 5d. The results reveal that after gelatin is prepared into nanoparticles, the ability of hydrogen bond to restore collagen structure is attenuated to a certain extent. In aqueous solution, the molecular chain of gelatin NPs is partially stretched with the particle as the core and crosslinked to form colloids. Compared with the pure gelatin system, the nano particle system has higher mechanical strength, better plastic deformation ability and weaker creep, and the viscoelastic-viscoplastic characterization is more consistent with the Nishihara model.

Extrusion process -results and verification of plastic deformation stage
In order to explore the critical extrusion shear strength of the material, a system with slower advancing speed, thinner nozzle, and lower modulus/viscosity is selected to verify and control the extrusion behavior in the stage that just can maintain the steady-state process (at this time, if the advancing force is reduced, the material cannot be extruded and return to linear elastic behavior). As shown in Fig. 6a, the relationship between the relative slip velocity and the friction coefficient proves that the critical slip friction can be achieved under the action of lower propulsion velocity or shear force. When the sliding friction is consistent with the maximum static friction, it is the critical extrusion state of 3D printing, and the critical value is the yield stress/strain of the material. In this stable state, the wall shear distribution of the material at the nozzle outlet is calculated by ANSYS fluent module, as shown in Fig. 6b (the material is gelatin NPs-water system of 5% mass fraction). It can be seen from the figure that the wall shear stress rapidly decays from the nozzle axis to the tube wall, showing a nearly concentric distribution. Due to the constraint effect of the material at the outlet, the wall shear stress near the outermost ring is less than that at the bottom of the syringe. This means that although this part of the material can fill the container, it basically does not flow under the shear effect, which is adverse to the extrusion process. It can also be inferred that there is a minimum amount of material in the syringe. If it is lower than this limit, the extrusion process will be much more difficult.
In order to verify the yield stress/strain of the material, the rheological amplitude scanning (1-1000% of strain range, 1 Hz of frequency) is used for measurement, from which the flow yield modulus, yield stress/strain, flow stress/ strain, and shear thinning behavior of the material can be obtained, as shown in Fig. 7.
The figure shows that the yield stress/strain of gelatin NP-water material (5% mass fraction) is 202 Pa/10.25% (linear viscoelastic below this value and viscoplastic beyond this value), which is consistent with the ANSYS fluent simulation results shown in Fig. 6b. The viscosity of the material decreases under the yield state, but it is still not fluidized (G′ is greater than G″). At this time, if the amplitude (shear strength) continues to increase, the composite viscosity of the material decreases rapidly and gradually presents liquid flow characteristics (G′ is less than G″). In order to verify the tendency, the transient measurement of 5% mass fraction gelatin NP-water system under different stress loading is carried out to get the relationship between creep/recovery curve and time. The creep results are shown in Fig. 8 (stress loading is 20 Pa, 50 Pa, 100 Pa, 200 Pa, and 500 Pa).
It can be seen from Fig. 9a-c that with the increase of stress loading, the steady-state strain of the material increases, and the creep behavior tends to be obvious. At the same time, because of the plasticity, the strain cannot completely return to the initial state after removing the stress loading, so that the residual strain also increases. When the loading stress is near the critical state, as shown in Fig. 9d (stress loading is 200 Pa), the material can still maintain its shape in the form of stable creep for a certain period of time. However, when the increasing strain exceeds the critical value (about 10% in the figure, which is consistent with the amplitude scanning results), the curve will go beyond the linear viscoelastic zone and enter the plastic deformation stage. In this case, the creep loses its stable state and continues to increase, and cannot be recovered after removing the loading. Furthermore, when the stress of the material significantly exceeds the critical state, as shown in Fig. 9e (the stress loading is 500 Pa), the curve will immediately enter the plastic deformation stage without behavior of recovery. b) a) Fig. 6 Numerical simulation of critical steady state extrusion (gelatin NPs-water system of 5% mass fraction). a relationship between relative slip velocity and coefficient of friction; b wall shear stress distribution in critical state extrusion

Shear thinning behavior and response characteristics
In order to verify the consilience between model calculation and experimental results, and to evaluate and compare the extrudability in different types of materials, rheological analysis with different proportions is carried out, including frequency scanning (measuring the modulus and composite viscosity of materials), amplitude scanning (measuring the LVE and yield stress/strain of materials), flow scanning (measuring shear thinning characteristics and response degree of materials), and step rate scanning (measuring shear thinning characteristics and response rate of materials).

Laponite-water system (inorganic-inorganic)
The rheological tests and results of laponite-water system with different mass fractions are shown in Fig. 9 and Table 3 (the mass fraction of 20-40% is used to ensure the solid shape and extrudable state of the material).
The results show that the material strength of laponitewater system increases with the increase of mass fraction. However, when the mass fraction reaches 40%, the volume fraction is relatively high (the loose bulk density of laponite is about 1.37 g/cm 3 ), which leads to a rather difficult mixing process due to the dense gel network formed by lamellar interpenetration. Although higher mass fraction can significantly increase the yield stress (271-2182 Pa), it does not significantly change the yield strain (4.034-6.606%), which indicates that the material system has good plasticity. This means that the initial propulsion force before extrusion is required to be high, but it is easier to extrude when reaching steady state.
In terms of extrusion lag, laponite is the main factor affecting the creep of the system. With the increase of mass fraction, both of the creep effect and extrusion delay are more significant. On the other hand, the high mass fraction of the material system needs higher propulsion, which also promotes the behavior of extrusion lag. In the aspect of model error, because the viscous behavior is more obvious at the lower fraction (20%) and the plastic behavior is more obvious at the higher fraction (40%), the deviation from the model is slightly higher at these two mass fractions, while the viscoelastic-viscoplastic properties of the material at the 30% fraction are in the best agreement with the model. Figure 10 and Table 4 show the rheological test and results of different mass fractions of alumina NPs-pluronic L61 system (the mass fraction of 20-40% is used to ensure the solid shape and extrudable state of the material).

Alumina NPs-Pluronic L61 system (inorganicorganic)
The results show a similar trend to lanponite-water system that the strength of alumina Ns-pluronic L61 system also increases with the increase of mass fraction of dispersed phase. Due to the hydrophilic hydrogen bond between alumina NPs and pluronic L61, the interaction between the dispersed phase and the continuous phase is relatively weak at a lower mass fraction (20%), and the whole system tends to Fig. 7 Rheological amplitude scanning of gelatin NP-water system (5% mass fraction) be the suspended phase with lower modulus and liquefaction tendency. At high mass fraction (40%), because of the increasing binding sites and high mechanical strength of alumina NPs itself, the modulus of the whole system increases rapidly, which is 1 order of magnitude higher than that of laponite-water system. However, at the same time, the yield stress/strain is obviously smaller than those of laponite-water system. It is indicated that alumina NP-pluronic L61 system has high strength, poor resilience, moderate viscosity, excellent shear thixotropy, and plastic deformation ability, which is also confirmed by the flow scanning curve.
In addition, the results also show that when the system is affected by shear force, it is easy to extrude because of the obviously changed viscosity and relatively lower yield strain. Thus, there is almost no extrusion lag phenomenon, and the viscosity decay rate does not increase significantly at higher mass fraction (although the yield stress increases, the creep behavior decreases, which counteracts the lag effect). In terms of model error, the system with lower mass fraction has relatively higher error due to poor plasticity. With the increase of mass fraction, the plastic behavior tends to be obvious, and the agreement with the model also increases. However, it can be speculated that this consilience will decrease with the further increase of mass fraction (the bulk density of alumina NPs is slightly lower than that of laponite, and the 40% mass fraction has not reached the critical value of mixing difficulty yet).

Gelatin-water system (organic network-inorganic)
When the dispersed phase is inorganic material, the interaction between particles is relatively weak, while when the dispersed phase is organic material, because of the entanglement effect of molecular chain, it is easy to form stable jelly like structure. Figure 11 and Table 5 show the rheological tests and results of the gelatin-water system (the mass fraction is 15-25% because of the low bulk density of gelatin).
It can be seen from the results that the shear thinning characteristic of gelatin-water system is not obvious under gel state. The mass fraction of organic dispersed phase  system is much lower than that of inorganic dispersed phase system under the similar solid degree (tanθ). The strength of gelatin-water system is also significantly lower than that of laponite-water system and alumina NP-pluronic L61 system. However, through the molecular chain stretching in aqueous solution, gelatin can gradually recover its collagen structure to form dense and firm network after gelation. This kind of structure endows the material system with strong elastic deformation ability, and can withstand quite high shear strength and shear strain. That is why it is difficult to be extruded due to the relatively wide LVE and poor plasticity. In the experiment, gelatin-water system almost has no extrudability, and its elastic behavior is far beyond the viscous and plastic behavior, which largely deviates from the model.

Gelatin NP-water system (organic nanoparticles-inorganic)
Due to the excellent biological characteristics of gelatin, it is hoped that the organic dispersed phase system can also be used in extrusion 3D printing. Figure 12 and Table 6 show the rheological test and results of gelatin NP-water system with different mass fractions after gelatin is prepared into  nanoparticles (in order to ensure the solidity similar to gelatin-water system, the mass fraction is 5-15%).
The results show that the strength of gelatin NP-water system is significantly higher than that of gelatin-water system at the same mass fraction. However, at a similar degree of solidity (tanθ), the strength of the system is lower than that of laponite-water system and alumina NP-pluronic L61 system. Although gelatin NP-water system has good elasticity as well, it can bear much lower shear strain/stress than gelatin-water system. It is proved that the plasticity of nanoparticle system is better than that of fully extended molecular chain system (easier to be extruded).
Although the creep behavior of gelatin NP-water system has been alleviated to a certain extent, it is still more significant than that of inorganic particle dispersed phase system, so it also has obvious extrusion lag phenomenon. With the increase of the mass fraction, the creep behavior changes obscurely, but the viscosity decay time increases rapidly due to significantly increased yield stress/strain. In terms of model error, similar to laponite-water system, when the mass fraction is moderate (10%), the viscoelastic-viscoplastic behavior of the material is most consistent with the model description. When the mass fraction is higher or lower than this value, the model error increases because of the more obvious viscous or plastic behavior.

Extrudability of materials
In the extrusion process of ideal steady-state fluid, the shear stress increases with the increase of material dynamic viscosity, advancing speed, and constraint degree (generally, it is the shape and size of nozzle). The relationship is expressed as follows: where is the shear stress at a point of nozzle section, is the dynamic viscosity of material, u is the advancing speed, and y is the size of nozzle (thus du∕dy is the extrusion velocity gradient).
For ideal Newtonian fluid, the relationship between shear stress and material dynamic viscosity, advancing speed, and nozzle size is shown in Fig. 13. It can be seen from the figure that among the three influencing parameters, the dynamic viscosity of the material changes exponentially on shear stress with the highest weight. The influence of the nozzle size on the extrusion speed is square (inversely proportional to the square of the cross section area) as the second weight on the shear stress. The influence of advancing speed on extrusion speed is only linear, which is the lowest weight of influence on shear stress.
However, for the 3D printing ink with shear thinning characteristics, the viscosity of the system and shear rate shows an approximate linear relationship; that is, the faster the shear rate, the lower the dynamic viscosity of the material. This means that there is a non-monotonic effect between the shear stress and the advancing speed. Taking alumina NP-pluronic L61 system as an example, the relationship between shear stress and advancing speed and nozzle radius is shown in Fig. 14.  It can be seen from the figure that when the advancing speed is relatively low, the shear stress increases with the increase of the advancing speed. Due to the high dynamic viscosity of the material, the extrusion speed gradient is the dominant factor, which leads to a very large shear stress required for the critical steady-state flow. With the increase of advancing speed, the viscosity of the material decreases, resulting in the attenuation of extrusion speed gradient effect. At this time, the viscosity of the material becomes the main cause of the shear stress, and the shear stress decreases gradually. It can also be concluded that if the advancing speed remains constant, the extrusion speed gradient will be significantly reduced by increasing the size of the nozzle (the extrusion speed is inversely proportional to the square of the radius ratio of the nozzle to the syringe), which also tends to reduce the shear stress. Figure 15 shows the relationship between the push force of the syringe (injector thrust) and the advancing speed/nozzle radius. It can be seen that at a low advancing speed or tiny nozzle size, the material needs a considerable amount of thrust to reach the steady-state extrusion, while in general 3D printing device, the thrust setting is mostly below 100 N. At this time, the material can only be maintained  in the elastic creep, but not reach the plastic deformation, so it cannot be extruded. The results show that due to the shear thinning property of the material, the extrusion speed gradient is attenuated by increasing advancing speed, thus gradually reducing the thrust required to maintain the steady-state extrusion. On the other hand, the influence of nozzle size on the required steady-state thrust is more significant. For example, when the current shear viscosity of

Conclusions
In order to quantitatively characterize the differences and characteristics of the extrudability of different material systems in extrusion 3D printing, a Nishihara creep model based on non-Newtonian fluid with shear thinning property is established to simulate the rheological properties of typical dispersed-continuous phase binary mixed material systems, including inorganic-inorganic inorganic-organic, organic network-inorganic, and organic nanoparticle-inorganic types. Through the numerical calculation and experimental verification of the elastic deformation stage and plastic deformation stage in the extrusion process, the extrusion characteristics and difficulty of different material systems are obtained and compared as follow: The creep constitutive equation of the material with shear thinning properties based on the Nishihara model is in good agreement with the experimental verification. The results show that the viscoelastic-viscoplastic characteristics of inorganic-organic system (alumina NP-pluronic L61) are in the best agreement with the constitutive model. Compared with the lamellar interpenetrating stacking structure of laponite, the viscoelastic-viscoplastic properties of the weak network structure composed of partially extended molecular chains of gelatin nanoparticles in aqueous solution are slightly more consistent with the model. Because the viscoelasticity of organic network-inorganic system (gelatin-water) is more significant than that of viscoplasticity, the agreement with the constitutive model is relatively low.
In the extrusion process of 3D printing materials with shear thinning properties, the shear stress is related to the advancing speed of syringe, material viscosity, and nozzle diameter. Because the dynamic viscosity of the material itself has an approximate linear relationship with the shear rate, in terms of the shear stress, the influence of advancing speed and nozzle diameter show antagonistic relationship in a certain range; that is, the velocity gradient is the dominant factor under the condition of lower extrusion speed, and the dynamic viscosity is the dominant factor under the condition of higher extrusion speed.
In terms of extrusion properties of various materials, inorganic-organic system has smaller yield stress/strain, and is easier to extrude under the same strength of other materials. The extrusion process of inorganic-inorganic system is less affected by the creep behavior of the material, but the overall strength of the system is higher, so the extrusion difficulty is inferior to that of inorganic-organic system. The strength of the organic nanoparticle-inorganic system is no less than that of the former two systems due to its partially stretched molecular chain. However, the larger yield strain makes the extrusion delay behavior more obvious, and the extrudability ranks the third. The organic network-inorganic system has good resilience, but poor plasticity, which means that the material can withstand high stress and strain. Due to the serious extrusion delay effect, it is not suitable for extrusion 3D printing process.