Molding and Structure Co-Analysis of Injection Molding Composites Considering Heterogeneity in the Thickness Direction

Abstract

During automobile crashes, the main behavior of plastics in the vehicle is bending deformation. Because of this, it is important to accurately calculate the edge stress and allowable stress on structural surfaces. For that purpose, it is necessary to consider the distribution of the fiber orientation angle and degree in the plate thickness direction. In a general coupled analysis for injection-molded materials, a material model is created based on the physical properties obtained from an arbitrary test piece. Therefore, even if the orientation angle in each part to be analyzed is taken into consideration, the degree of orientation is constant in each part, and the degree of orientation distribution in the plate thickness direction is also treated as constant. We tried to improve the accuracy of crash analysis by calculating the edge stress and allowable stress with high accuracy by considering the distribution of orientation angle and degree in the plate thickness direction. For that purpose, we created a model that calculates properties of composite materials from resin and fiber properties, and a process that maps the distribution of orientation angle and degree in the plate thickness direction to each integration point of the shell element in structural molding analysis.

 

This Paper was Originally Published in Japanese in J. Soc. Mater. Sci., Japan 72 (2023) 675–682.

1. Introduction

Injection-molded resins such as polypropylene are widely used for automobile interior materials because of the high degree of freedom of shape and molding cycles required.¹⁾ In recent years, there has been a need to reduce vehicle weight to reduce CO₂ emissions during driving, and fiber-reinforced resins with high specific stiffness and specific strength are being applied.²⁾ In addition to considering the stress triaxiality of yield and fracture properties, models such as the anisotropic yield functions of Hill (1984) and Barlat (2003),^3,4) which change yield properties depending on the fiber orientation angle, are now used for the simulation-based design of these resin parts. According to this trend, coupled analysis techniques have been developed to include the fiber orientation angle in structural analysis based on injection molding analysis.⁵⁾ Coupled analysis is also used in automotive development, and the parameters for the materials model are based on materials properties obtained from test specimens. Therefore, the analysis is performed based on the assumption that the fiber orientation is constant and the degree of fiber orientation is the same in all regions of the specimen and parts. However, as shown in Fig. 1, the fiber orientation is generally different in each region of the part. It is also known that even the same part has a distribution of orientation angle and degree in the direction of the plate thickness, such as skin and core layers, due to fountain flow during forming.⁶⁾ There are few verified cases of the effects of these orientation angles and distributions within a part and in the thickness direction on the accuracy of crash performance prediction.

Fig. 1

Distribution of fiber orientation in in-plane and out-of-plane.

Because many parts of an automobile are subjected to bending deformation in vehicle collisions, it is important to accurately predict the stiffness and strength of the parts during bending deformation. Accurate prediction of strength under bending deformation requires accurate prediction of the edge stresses at the surface of the member, which is the starting point for fractures. In contrast, conventional methods treat the material as homogeneous in the thickness direction, so its properties are often determined based on in-plane properties such as those obtained by tensile tests. In other words, the physical properties of the skin layer and core layer are assumed to be homogeneous as a whole. However, as shown in Fig. 2, the skin layer on the surface of an actual part has a high degree of orientation, and the conventional method determines the onset of fracturing at a lower edge stress than that for the actual phenomenon.

Fig. 2

Difference between modeling the characteristics of thickness direction as a homogeneous body and heterogeneous body.

In this study, the influence of the degree of orientation on material properties was clarified experimentally, and a coupled analysis that accounts for the distribution of the degree of orientation and orientation angle in the in-plane and plate thickness directions was performed. The accuracy of the analysis was verified for application to the development of vehicle. A new mapping system was developed and a coupled analysis process was examined in combination with a material model that can account for the orientation degree.

2. Overview of Coupled Analysis

To predict the strength and stiffness of injection-molded composites with high accuracy, a coupled process of forming and structural analysis was developed that can account for heterogeneity in the thickness direction. Because a crash analysis for a single vehicle requires several million elements, an explicit process using first-order shell elements was used to reduce the computational cost. An overview of the coupled analysis developed in this study is shown in Fig. 3.

Fig. 3

Process of co-analysis for injection molded composite.

Moldex-3D was used for the molding analysis and LS-Dyna for the structural analysis. In the molding analysis, the Folgar-Tucker fiber orientation equation is used to estimate fiber orientation from the resin flow and convective flow. For each integration point of a solid element in the molding analysis, the fiber orientation is obtained from the second-order tensor.⁷⁾ From this orientation tensor, eigenvectors are obtained as the principal direction of fiber alignment, and eigenvalues are obtained as the statistical ratio of fibers aligned in that direction. The oriented ellipsoid defined by the eigenvalues and eigenvectors is projected onto the plane of the shell element in the structural analysis to obtain a plane ellipsoid showing the fiber orientation in 2D. From this plane ellipsoid, the eigenvalues of the vectors are obtained as the orientation degree, and the angle between the vectors and the reference angle of the shell element is obtained as the orientation angle. By doing this for each integration point for all shell elements, we obtain a structural analysis model that reproduces the distribution of orientation angle and orientation degree in the in-plane and thickness directions for the part being investigated.

2.1 Stiffness representation of composite materials

For composite materials, it is necessary to determine not only the anisotropy but also the triaxial stress distribution of the yield and rupture properties. Therefore, it is common to use models such as MAT54, 58, 157, 215, 261, and 262, when analyzing with LS-Dyna.^8,9) In this study, we used MAT215, which can predict composite properties based on the properties of the fibers and base resin, to represent the dependence of the properties on the degree of orientation. This model uses the Mori-Tanaka model to predict composite properties.¹⁰⁾ The average stress and strain for the composite are given in eqs. (1) and (2), and are decomposed into fiber and resin stress and strain based on the fiber content φ. By introducing the concentration tensor B, the relationship between the average stress and the average strain for the fiber and the resin is given by eqs. (3) and (4).

The concentration tensor B is defined using the stiffness S and the compliant tensor C, as shown in eqs. (5), (6), and (7).

2.2 Representation of composite damage and fracture

The MAT215 model used in this study represents resin yield and damage, as well as fiber tensile rupture. The Damage Initiation and Evolution Model¹¹⁾ was used to consider resin damage. A damage parameter D is introduced and a linear relationship is defined as shown in eq. (8), where the stress decreases as damage progresses. The evolution of D introduces a plastic displacement U_p and is given in eq. (9). When U_p reaches the critical plastic displacement U^p_f, D becomes 1, rupture occurs, and the element is deleted. U_p is given in eq. (10) in terms of the element representative length l and the plastic strain ε^P when the damage initiation parameter ω_D is introduced and its value is greater than 1. The parameter ω_D, which determines the onset of damage, is given in eq. (11). When the sum of the plastic strain ε_p in each stress state reaches the critical plastic strain ε^p_f, ω_D is set to 1 and damage is initiated.

The critical plastic strain ε^p_D is expressed in terms of the triaxial stress η as shown in eqs. (12) and (13).

3. Verification of Accuracy of Molding Analysis

Table 1 shows the forming conditions and specifications for the flat plates. Forming analysis was performed to obtain fiber orientation information at the cutout position in each specimen. Because errors in orientation information obtained from the forming analysis affect the prediction errors in the structural analysis, nondestructive testing was used to verify the accuracy of the analysis. Talbot-Lau imaging was used to determine fiber orientation over the entire plate area, and X-ray micro-computed tomography (micro-CT) was used to determine the local fiber orientation.

Table 1 Test piece spec and molding condition.

3.1 Evaluation of fiber orientation by Talbot-Lau imaging

For Talbot-Lau imaging of injection-molded composites, we used the technique described by Yoshimura et al.¹²⁾ The principle is shown in Fig. 4. In this technique, a diffraction grating is installed to obtain equivalent X-ray images at 0°, 45°, and 90° to obtain fiber orientation information. The relative fiber content is indicated by brightness in Fig. 5, and the orientation averaged in the thickness direction is shown as brightness and the orientation angle as color contours in Fig. 6. The flow ends were excluded from the evaluation of physical properties due to their high relative fiber content, which was expected to differ from the fiber content at the center. The degree of orientation at the center and the side of the flat plate was significantly different, indicating that specimens with different degrees of orientation can be obtained depending on the cutout position of the specimen in the mold direction (MD). To evaluate the orientation angles more quantitatively, image analysis was conducted on the same area as the specimen’s gauge area, from Area 1 to Area 5. The frequency of fiber orientation angles averaged in the thickness direction in the area of interest is shown in Fig. 7. The black solid line shows the frequency of orientation angles for the entire flat plate area as a reference. In Area 1 and Area 5, most of the fibers were oriented at 0 degree, whereas in Area 3, the fiber orientation angle was distributed in the range up to ±45 degrees.

Fig. 4

Principle of Talvo-Low interference measurement.

Fig. 5

Relative fiber contents in a plate.

Fig. 6

Distribution of fiber orientation angle and strength of orientation averaged in the plate thickness direction.

Fig. 7

Frequency distribution of averaged fiber orientation angle in each evaluation area.

3.2 Evaluation of local fiber orientation by X-ray micro-CT

Talbot-Lau imaging cannot distinguish skin and core layers because the information in the thickness direction is superimposed. Therefore, X-ray micro-CT was used to obtain fiber information. Table 2 shows the evaluation conditions, and Fig. 8 shows the results for the fiber orientation tensor A_ij with A₁₁ as the degree of orientation in the MD direction. The evaluation range was the same as that in Fig. 7, with three evaluation points (outside, intermediate, and center, corresponding to Areas 1, 2, and 3, respectively). As expected, the degree of orientation varied in the thickness direction, with the skin layer having a higher degree of orientation than the core layer, and the outside layer having a higher degree of orientation than the center layer over the entire thickness range. The latter result was consistent with the Talbot-Lau imaging results and is considered to be reasonable. The length-weighted average fiber length at each cutout position was only 600 to 700 µm, and the fiber length distribution in the in-plane direction could not be determined, and was therefore not considered in this model.

Table 2 Measurement conditions of micro x-ray CT.

Fig. 8

Orientation tensor from image analysis of micro x-ray CT.

3.3 Prediction of orientation information of base plate

The viscosity model for the forming analysis was the general Cross-WLF model.¹³⁾ The defining equations are shown in eq. (14), where η is the viscosity, T is the temperature, $\dot{\gamma }$ is the shear rate, p is the pressure, η₀ is the zero-shear viscosity, n is the flow index, T₀ is the reference temperature, D₃ is a pressure-dependent coefficient, and τ^*, D₁, A₁, and $\tilde{A}_{2}$ are fitting coefficients.

The viscosity characteristics set by the model are shown in Fig. 9, and the analysis conditions are shown in Table 3. To verify the prediction accuracy for the fiber orientation, the orientation tensor information in the thickness direction in the center, intermediate, and outside ranges was compared with the results calculated by the X-ray micro-CT described in Section 3.2 and Fig. 10. In the experiment, the most highly oriented area was 0.3 mm from the surface layer to the center of the plate thickness, and the analytical results are consistent with this trend. However, the analysis tended to underestimate the degree of orientation in the transverse direction (TD).

Fig. 9

Pressure dependence of viscosity.

Table 3 Molding analysis condition.

Fig. 10

Comparison of experimental of micro x-ray CT and molding analytical results.

4. Obtaining Physical Properties of Injection Molded Composites

Using the above-mentioned plates, physical properties were obtained at different orientation angles and degrees of orientation by tensile and compression tests at different cutout angles and positions. Three-point bending tests were also conducted to verify the accuracy.

4.1 Tensile test

In general, composite materials are tested using a specimen with a large area of the gauge area to reduce variability in the material properties. However, in this study, because it was necessary to model local changes in fiber orientation, tensile tests were conducted in accordance with ASTM D 1822 Type L, which requires a small area of the score section. In addition, to determine the effect of the size of the score section, an original specimen with a large score section was also used for evaluation. The number of tests in both cases was n = 5. A hydro shot tester was used to obtain the velocity dependence. Table 4 shows the test conditions, Fig. 11 shows the cutout positions, and Fig. 12 shows the test results.

Table 4 Tensile test condition.

Fig. 11

Conditions for cutting out test pieces from a flat plate.

Fig. 12

Experimental results of tensile strength in ASTM and Original specimen at 0° and 90° of cut angle.

As shown in Fig. 12(a), the tensile properties in the 0° direction clearly differ depending on the cutout position. The specimen cut from outside near the edge of the plate had higher modulus and strength, while the stiffness and strength decreased towards the center of the plate. This is considered to be due to the higher orientation at the edge of the specimen based on the results of Talbot-Lau imaging and X-ray micro-CT. In contrast, as shown in Fig. 12(b), there was little dependence of the elastic modulus in the 90° direction on the cutout position. The Talbot-Lau imaging results showed that there is little change in the degree of orientation in the flow direction, even if the position in the specimen is different, which is thought to be due to the lower difference in the degree of orientation of the 90° specimens depending on the cutout position. As shown in Fig. 12(c), the ASTM specimens with a smaller grading section showed a larger variation in physical properties, and the static strength tended to be slightly higher. This indicates that the model, which is intended to obtain local properties, should be evaluated using specimens with a small grade section. The results for the strain rate dependence are shown in Fig. 13. The strain-rate dependence was low in both the 0° and 90° directions, and the material model did not consider the strain-rate dependence.

Fig. 13

Confirmation of strain rate dependence.

4.2 Compression test

A tensile specimen based on ASTM D 1822 Type L was gripped using a highly rigid chuck without a universal joint to avoid buckling, and a compression test was performed. The evaluation conditions are shown in Table 5.

Table 5 Conditions for compression test.

The relationship between the measured compressive properties in the flow direction and the specimen cutout position is shown in Fig. 14. As in the tensile test, the elastic modulus and compressive strength were high for the outside position, while the elastic modulus and compressive strength were low for the center position. This confirms that there are differences in compressive properties depending on the degree of orientation.

Fig. 14

Relationship between cut out position and properties on 0deg compression test.

4.3 Bending test

To verify the accuracy of the analysis, a three-point bending test was conducted in accordance with JIS K7171. The pitch of the supports was 48 mm, where L/h = 16. Table 6 shows the test conditions, and Fig. 15 shows the load-displacement characteristics of the loader. The bending stiffness and strength in the 90° direction were higher than those in the 0° direction, confirming anisotropy. As in the tensile test, there was a difference in stiffness and strength in the 0° direction depending on the cutout position, whereas there was no significant difference in the 90° direction.

Table 6 Conditions for three-point bending test.

Fig. 15

Force-displacement curve of three point bending on 0° test.

5. Parameter Identification of Material Model

The parameters for MAT215 were identified based on the results of the tensile and compression tests described in Section 4. For comparison with the conventional method, which considers the thickness direction as being homogeneous, test results at the center and intermediate cutout positions were used to identify the parameters for MAT157, a general orthotropic anisotropic elastic-plastic model.⁸⁾ The identification procedure is shown in Fig. 16. The analytical model used in the identification is shown in Fig. 17, which simulates a tensile test as a representative example.

Fig. 16

Material model correlation procedure considering the fiber orientation strength of the material test pieces.

Fig. 17

Structural analysis model assuming tensile test.

Each model simulated the specimen geometry, and stresses and strains were obtained from the reaction forces generated in the rigid body corresponding to the specimen chuck and the displacements of the nodes at the evaluation points. For the compression test, the rigid body range was modified to match the chuck position. The distance between nodes was set to 2.1 mm, which matched the length of the strain gage used in the experiment. Ten integral points in the thickness direction of the shell element were used to match the orientation measurement results of the X-ray micro-CT test. For material parameters, the manufacturer’s nominal value¹⁴⁾ was used as the initial value for fiber properties and the value of PA6-GF45 by Reithofer et al.¹⁵⁾ was used as the initial value for resin properties. Fiber orientation information was obtained from the results of the flow analysis of the base plate, and material model parameters were adjusted by performing an analysis simulating a material test at each cutout position. The parameters for the material model are shown in Figs. 18 and 19. Parameters that could be unambiguously determined from physical properties tests or actual measurements are shaded in gray, and the parameters adjusted as described above are shaded in blue. The definitions of the parameters for MAT215 are as follows: W_f is the fiber content, length is the typical fiber length, φ is the fiber diameter, ρ_m is the density of the resin, E_m is the elastic modulus of the resin, ν_m is the Poisson ratio of the resin, ρ_f is the fiber density, E_xx is the axial modulus of the fiber, E_yy is the transverse modulus of the fiber, G_xy is the shear modulus of the fiber, ν_xx is the longitudinal Poisson ratio for the fiber, ν_yy is the orthogonal Poisson ratio for the fiber, C_ij is the stiffness matrix of the composite, and R₀, R₄₅, and R₉₀ are the Hill constants³⁾ of the Hill 1948 yield function. The GISSMO model,¹⁶⁾ which is commonly used when considering triaxial stress, was used to define the fracture properties.

Fig. 18

Material parameter of MAT215.

Fig. 19

Material parameter of MAT157[4].

The reproducibility of each material test mode is shown in Fig. 20. It was found that not only anisotropy but also differences in properties depending on the cutout position were reproduced. However, to reproduce the stiffness in the 90° direction, the orthogonal modulus of the fibers had to be higher than the original value. This is considered to be one of the reasons for the low prediction accuracy for the orientation in the TD direction.

Fig. 20

Results of confirmation of the material test.

6. Verification of Analysis Accuracy

To verify the accuracy of the newly developed coupled analysis process, a three-point bending test was performed. The analysis results using MAT157, which was assumed to be homogeneous in the thickness direction, are shown in Fig. 21(a), and the analysis results using MAT215 are shown in Fig. 21(b). Table 7 shows the analytical error. The new method using the MAT215 model was able to predict changes in stiffness and strength depending on the cutout position and angle, with a maximum analytical error of 5.4% for the maximum load, and 3.4% for the rupture strain. In contrast, MAT157 had a maximum analysis error of 33.2% for the maximum load and 20.8% for the rupture displacement, confirming the superiority of the new method.

Fig. 21

Verification of analysis accuracy.

Table 7 Analysis error at 3 points bending.

To determine the reason for the improved accuracy of the analysis, the difference in stress in the thickness direction is shown in Fig. 22. At the center position, the maximum stress for the conventional method was 292 MPa, while the maximum stress for the new method was 429 MPa, indicating that the new method can express the increase in strength due to the high orientation of the skin layer. In the outside position for the new method, the maximum stress occurred in layers 2 and 9, which is consistent with the tendency of the orientation distribution in the thickness direction shown in Fig. 8, where the degree of orientation was highest in a layer with a depth of about 0.3 mm from the top surface. This point also confirms that consideration of thickness directional heterogeneity contributes to the accuracy of the analysis.

Fig. 22

Comparison results of edge stress on 0deg three point bending at generate maximum force.

7. Conclusion

We showed that conventional coupled analysis methods for injection-molded composites have low prediction accuracy in principle for edge stresses during bending deformation. To improve the prediction accuracy, it is necessary to reflect changes in physical properties due to the distribution of orientation angle and orientation degree in the thickness direction in the structural analysis, and we proposed a coupled analysis procedure for this purpose. To validate the proposed method, the accuracy of the fiber orientation tensor obtained from the injection-molding analysis was confirmed using Talbot-Lau imaging and X-ray micro-CT, and the material model was identified by obtaining the relationship between physical properties and orientation angle and orientation degree from material tests. Based on these results, we confirmed that the proposed method improves the accuracy of the analysis of the three-point bending test, with an analysis error of approximately 5% or less, compared with an error of 20–30% range for the conventional method. To further improve the accuracy, it is going to be necessary to improve the prediction accuracy of the fiber orientation in the forming analysis.

Acknowledgments

This work was performed as a commissioned project (JPNP14014) of the New Energy and Industrial Technology Development Organization. We thank JSOL Corporation for the mapping system and Mechanical Design Inc. and Nagoya Municipal Industrial Research Institute for their assistance in obtaining the physical properties, Dr. Tsuyoshi Matsuo of the University of Tokyo for analysis and discussion of the obtained physical properties, and Konica Minolta Inc. for Talbot-Lau imaging.

REFERENCES

• 1)   J.  Agarwal,  S.  Sahoo,  S.  Mohanty and  S.  Nayak: J. Thermoplast. Compos. Mater. 33 (2020) 978–1013. doi:10.1177/0892705718815530
• 2)  K. Lisk, C. Soyarslan, H. Richter and A. Erman: “Analysis of formability of advanced high strength steel sheets with phenomenologically based failure criteria with separate treatment of instability, shear and normal fracture”, 8th European LS-DYNA Users Conference, Strasbourg, (2011).
• 3)   R.  Hill: Proc. R. Soc. London, Ser. A 193 (1948) 281–297. doi:10.1098/rspa.1948.0045
• 4)   F.  Barlat and  J.  Lian: Int. J. Plast. 5 (1989) 51–66. doi:10.1016/0749-6419(89)90019-3
• 5)  K. Yunoki, K. Yamamoto, N. Hirayama and K. Terada: Proceeding of the Conference on Computational Engineering and Science 22, (2017) p. 4.
• 6)   R.  Nakano,  T.  Shiota,  A.  Hirakimoto,  K.  Taki,  H.  Takashio and  C.  Yong-Sung: Seikei-Kakou 20 (2008) 237–241. doi:10.4325/seikeikakou.20.237
• 7)   S.  Yashiro,  T.  Okabe and  K.  Matsushima: Journal of the Japan Society for Composite Materials 36 (2010) 11–18. doi:10.6089/jscm.36.11
• 8)  W. Korte, S. Pazour and M. Stojek: LS-Dyna Forum 2016, https://www.dynamore.de/de/download/papers/2016-ls-dyna-forum/Papers 2016/Dienstag-11.10.16/.
• 9)  A. Cherniaev, J. Montesano and C. Butcher: 15th International LS-DYNA Users Conference, 03-A, (2018) pp. 1–17.
• 10)  P. Reithofer, A. Fertschej, B. Hirschmann, B. Jilka, A. Erhart and S. Hartmann: 15th International LS-DYNA Users Conference, 45-B, (2018) pp. 1–19.
• 11)  T. Borrvall, T. Johansson and M. Schill: “A General Damage Initiation and Evolution Model (DIEM) in LS-DYNA”, 9th European LS-DYNA Conference, (2013).
• 12)  A. Yoshimura, Y. Kasai, M. Arai, K. Goto, M. Ichiki, A. Yamanaka and H. Shindo: Proceedings of the 2019 Annual Meeting of the Japan Society of Mechanical Engineers, No. 19-1, 2019.9.8–11.
• 13)   S.  Goto: Seikei-Kakou 31 (2019) 297–301. doi:10.4325/seikeikakou.31.297
• 14)  Toray carbon fiber ameria, inc, Technical data sheet, No. CFA-005.
• 15)  P. Reithofer, A. Erhart, A. Fertschej, S. Hartmann and B. Jilka: “MAT_4A_MICROMEC micro mechanic based material model”, Automotive CAE Grand challenge, (2017).
• 16)  D. Horling: “Parameter identification of GISSMO damage model for DOCOL 1200M”, Faculty of science and technology.

Appendix

  

\begin{equation} \bar{\sigma}^{\textit{Composite}} = \varphi \bar{\sigma}^{\textit{Fiber}} + (1 - \varphi)\bar{\sigma}^{\textit{Matrix}} \end{equation} (1)

  

\begin{equation} \bar{\varepsilon}^{\text{Composite}} = \varphi \bar{\varepsilon}^{\textit{Fiber}} + (1 - \varphi)\bar{\varepsilon}^{\textit{Matrix}} \end{equation} (2)

  

\begin{equation} \underline{\bar{\sigma}}^{\text{Fiber}} = \underline{\underline{\text{B}^{\sigma}}}\, \underline{\bar{\sigma}}^{\textit{Matrix}} \end{equation} (3)

  

\begin{equation} \underline{\bar{\varepsilon}}^{\text{Fiber}} = \underline{\underline{\text{B}^{\varepsilon}}}\, \underline{\bar{\varepsilon}}^{\textit{Matrix}} \end{equation} (4)

  

\begin{equation} \underline{\bar{\sigma}}^{\text{Fiber}} = \underline{\underline{S^{\textit{Fiber}}}}\, \underline{\bar{\varepsilon}}^{\textit{Fiber}} \end{equation} (5)

  

\begin{equation} \underline{\underline{\text{B}^{\sigma}}} = \underline{\underline{S^{\textit{Fiber}}}}\, \underline{\underline{\text{B}^{\varepsilon}}}\, \underline{\underline{C^{\textit{Matrix}}}} \end{equation} (6)

  

\begin{equation} \underline{\underline{C^{\textit{Matrix}}}} = (\underline{\underline{S^{\textit{Matrix}}}})^{-1} \end{equation} (7)

  

\begin{equation*} 0 \leqq \text{D} < 1 \end{equation*}

  

\begin{equation} \sigma = (1 - \text{D})\tilde{\sigma} \end{equation} (8)

  

\begin{equation} \dot{D} = \frac{\dot{u^{p}}}{u_{f}^{p}} \end{equation} (9)

  

\begin{equation} \dot{u^{p}} = \begin{cases} 0 & \text{$\omega_{D} < 1$}\\ l\dot{\varepsilon^{p}} & \text{$\omega_{D} \geqq 1$} \end{cases} \end{equation} (10)

  

\begin{equation} \omega_{D} = \int_{0}^{\varepsilon^{p}}\frac{d\varepsilon^{p}}{\varepsilon_{D}^{p}} \end{equation} (11)

  

\begin{equation} \varepsilon_{D}^{p} = \varepsilon_{D}^{p}(\eta,\dot{\varepsilon^{p}}) \end{equation} (12)

  

\begin{equation} \eta = -\text{p}/\text{q} \end{equation} (13)

  

\begin{equation} \eta(\text{T},\dot{\gamma},\text{p}) = \cfrac{\eta_{0}(T,p)}{1 + \biggl(\cfrac{\eta_{0}(T,p)}{\tau^{*}}\dot{\gamma}\biggr)1 - n} \end{equation} (14)

  

\begin{align*} &\eta_{0}(T,p) = D_{1}\,\mathit{exp} \left[-\frac{A_{1}(T - T_{0})}{A_{2} + (T - T_{0})}\right]\\ &T_{0} = D_{2} + D_{3}\text{p}\\ &A_{2} = \tilde{A}_{2} + D_{3}\text{p} \end{align*}