Intelligent Technologies for Polymer Injection MoldingView this Special Issue
Research Article | Open Access
Simulation and Optimization of Short Fiber Circumferential Orientation in Short-Fiber-Reinforced Composites Overflow Water-Assisted Injection Molded Tube
The mechanical properties of the water-assisted injection molded tube can be enhanced by the increase in the short fiber circumferential orientation (SFCO). Thus, the numerical method verified by experiments is used to simulate the SFCO distribution in the overflow water-assisted injection molding (OWAIM), with the mechanism of short fiber orientation analyzed as well. The effect of parameters (filling time, melt temperature, mold temperature, delay time, water pressure, and water temperature) on the SFCO is explored by range analysis and variance analysis of the orthogonal experimental scheme. Moreover, both of artificial neural network (ANN) and genetic algorithm (GA) are used to model and optimize process parameters. Results show that the melt temperature, delay time, and water pressure are predominant parameters. The evolution of SFCO increases with the increase of melt temperature and water pressure, whereas the changes in delay time reverse. The value of the maximum SFCO tensor obtained by GA optimization is found to be 0.234.
Injection molding is the main process used to produce plastic products [1, 2]. Water-assisted injection molding (WAIM), similar to the gas-assisted injection molding (GAIM), fills the high-temperature melt followed by the injection of high-pressure water after a short delay time and forms a plastic part with hollow channel finally. Due to the incompressibility, high thermal conductivity, and high heat capacity of the medium water, WAIM has the characteristics of thin residual wall thickness (RWT) and high production efficiency and has obvious technological advantages in the preparation of hollow shaped plastic parts. At present, the research mainly focuses on the influence of factors (such as process parameters and material properties) on primary penetration, secondary penetration, and RWT distribution of plastic parts [3–6].
Short-fiber-reinforced composites (SFRC) can effectively improve the mechanical properties of water-assisted injection molded parts . Liu et al.  have used the short-fiber reinforced polypropylene as raw material in short water-assisted injection molding (SWAIM). Results showed that the short fibers in RWT oriented mostly along the melt flow direction. Meanwhile, the orientation of the poly-butylene-terephthalate composites containing 15% glass fiber has been studied by Liu et al. . They have found that the short fibers near the mold surface are arranged mainly along the melt flow direction, and the orientation gradually decreased with the increase of thickness. In order to further explore the orientation distribution characteristics of short fibers in SWAIM, Huang et al.  have proposed that the short fiber orientation along the melt flow direction in the cross section near the high-pressure water inlet decreases gradually with the increasing thickness.
The mechanical properties of SFRC are closely related to the orientation of short fibers, and the reinforcing effect is mainly controlled by the direction of short fibers arrangement [11–13]. Experimental studies show that the short fibers of RWT in WAIM are mainly oriented along the axial direction and weakly oriented in the circumferential and radial directions. This implies that the mechanical properties of the parts in the circumferential direction are far below the axial direction. In order to improve the pressure resistance, friction resistance, and aging resistance of the water-assisted injection molded pipe, it is desirable to have more short fibers oriented in the circumferential direction. In general, the orientation distribution of the short fibers can be adjusted by altering the melt flow pressure and velocity field distribution in WAIM through appropriate process parameters [14–16].
Up to our knowledge, due to the difficulty in assessment of the performing quantitative experimental characterization, various researches focus on the qualitative observation of the short fiber orientation using scanning electron microscopy (SEM). To quantitatively characterize the short fiber orientation of RWT, it is crucial to carry out effective numerical simulation for WAIM. Thus, based on the improved Anisotropic Rotary Diffusion model and Retarding Principal Rate Model (iARD-RPR) proposed by Tseng [17–19], the numerical simulation of short fiber orientation distribution in overflow water-assisted injection molding (OWAIM) is implemented. The validity of the model and mechanism of short fiber orientation are analyzed. For this purpose, the orthogonal experimental scheme is designed to obtain the sample data. Meanwhile, an artificial neural network (ANN) is used to map the nonlinear relationship between process parameters and short fiber circumferential orientation (SFCO). In addition, the genetic algorithm (GA) uses the ANN model as the fitness function to optimize the process parameters for the maximum circumferential orientation.
2.1. Related Mathematical Models
The numerical simulation of short fiber orientation in OWAIM is based on the fiber orientation tensor evolution equation and the fluid mechanics governing equation. According to the fiber orientation theory, a single fiber is usually regarded as a rigid cylindrical rod. As shown in Figure 1, the orientation of a single fiber can be described by the unit vector .where P is the direction of the fiber; θ and φ are the angles between P and the coordinate axis.
The orientation tensor A is defined to succinctly depict the orientations of a large amount of short fibers .where ψ(P) is the probability density distribution function of the whole orientation space; A is a symmetric matrix, and when A=I/3, it represents the orientation state isotropic, where I is an identity matrix. The diagonal Components A11, A22, and A33 represent orientation tensors in the axial, circumferential, and radial directions, respectively.
Tseng has proposed the iARD-RPR fiber orientation prediction model composed of three parts, such as hydrodynamic model , improved anisotropic rotational diffusion model , and delayed principal rate model .where contains two effective parameters: interfiber interaction factor and fiber-matrix interaction factor ; contains a parameter α (0 < α < 1), used to slow down the response speed of fiber orientation. W is the vortex tensor; D is the deformation rate tensor; ξ is a dimensionless number. The fourth-order orientation tensor is determined by the higher order polynomial approximation of second-order tensor A; is the diagonal tensor material derivative and its superscript is the intrinsic orientation kinetic assumption; is the i-th diagonal component of ; R is a rotation matrix; is the eigenvalue of the matrix A ().
In the injection molding process, the initial condition of the fiber orientation tensor at the entrance was set as the isotropic state. The movement of short fibers in the polymer melt is a transient, non-Newtonian, and nonisothermal process. In the numerical simulation, the melt is regarded as incompressible, laminar, and the inertia term is ignored. The basic governing equations for transient and nonisothermal melt motion in OWAIM include the following :where ρ is the melt density; u is the velocity vector; t is the time; is the total stress tensor; is the gravity acceleration vector; η is the viscosity; P is the pressure; is the specific heat capacity; T is the temperature; is the thermal conductivity; is the shear rate.
The seven-parameter Cross-WLF viscosity model is used as the constitutive equation, which can well describe the nonlinear relationship between polymer melt viscosity, temperature, and shear . where η is the viscosity; is the zero shear viscosity; is the shear rate; is the material constant; is the power law index in high shear rate; T is the melt temperature; is the glass transition temperature; , , , , and are the Constants associated with the selected material.
2.2. Geometric Model and Process Parameters
The geometric model (Figure 2) consists of two parts: an overflow cavity and a plastic part with length 280 mm and diameter 20 mm. The 3D solid combination model is built using Pro/E and meshed in Moldex Designer. The total number of mesh nodes and mesh elements is 205,056 and 2,208,359, respectively. The 3D solid meshed model is imported into the commercial software Moldex3D for material selection and process parameter settings. The PP (Fiberfil J-68/30/E) containing 30% short glass fiber with the fiber aspect ratio 20 is used as material in simulation of OWAIM. First, the melt was injected into the mold cavity of the plastic part. Second, after a short delay time, the high-pressure water was injected into the melt and penetrates along the position with the least flow resistance and pushed the melt into the overflow cavity to form a plastic part with a hollow cross-section. The process parameters considered in this study are shown in Table 1.
2.3. Experimental Verification
The experiment of OWAIM was carried out in the laboratory using the same process parameters as the simulation to verify the prediction result of short fiber orientation distribution. The platform includes an injection molding machine, a high pressure water injection system, a mold temperature control system, and a water temperature control system. The injection molding machine is the TTI-250FT automatic injection molding machine produced by Donghua Machinery Co., Ltd. The high-pressure water injection system is composed of a control system, a high pressure pump, and a water injection nozzle. The maximum water injection pressure is 20 MPa. The raw material is short-fiber-reinforced polypropylene (PP/GH43 with a short glass fiber mass fraction of 30%) produced by South Korea's Samsung Total Co., Ltd.
As shown in Figure 3, a section was taken in the middle of the plastic part along the axial direction. The orientation distribution of the short fibers was observed by SEM (Nova NanoSEM 450) with an acceleration voltage of 5kV. The sample was first immersed in liquid nitrogen and cryofractured after two hours. The prepared specimen was gold-sputtered before the observation.
2.4. Orthogonal Experimental Design
During the injection molding process, the melt is sheared and stretched by the surrounding melt due to the difference in viscosity gradient, pressure, and velocity field distribution. The shearing action leads the short fibers to be aligned along the melt flow direction, and the stretching effect induces the short fibers to orient in the stretching direction. The process parameters of OWAIM considered in this research include filling time, melt temperature, mold temperature, delay time, water pressure, and water temperature. Set the range of values of each process parameter based on the software's recommendation for the selected material. In order to reduce the number of experiments and comprehensive investigation of the effect of process parameters on the SFCO, orthogonal experimental design (L25 (56)) with six factors and five levels is arranged, as shown in Tables 2-3. The objective A22 represents the fiber orientation in the circumferential directions of the pipe. The values of circumferential component of orientation tensor were taken at eighteen different points along the thickness direction of RWT, and A22 is the average value of these points.
2.5. Artificial Neural Network (ANN)
ANN inspired by biologic neural system is a computing model used to map linear or nonlinear relationships between factors and responses. ANN can work as a human brain to establish a sample model from the perspective of information processing, without prior information or heuristic assumptions. An ANN model comprises three parts: one input layer, one or more hidden layers, and one output layer. The numbers of neurons in the input and output layer are determined by the numbers of factors and responses, respectively. In general, the number of neurons in hidden layer is determined by trial and error. The relationship between factors and responses can be depicted as follows:where Ii denotes the i-th factor in input layer, wji denotes the weight between the i-th neuron in input layer and the j-th neuron in hidden layer, wkj denotes the weight between the j-th hidden layer and the k-th output layer, bj and bk denote the bias assigned to j-th neuron in hidden layer and k-th neuron in output layer, respectively, fh denotes the transfer function employed in the hidden layer, fo denotes the transfer function employed in the output layer, and Ok denotes the k-th response in output layer.
According to Kolmogorov’s theorem , an ANN model with a single hidden layer has the ability to map any complex nonlinear relationship between the factors and responses. In this study, an ANN model with one hidden layer was used for modeling circumferential component of orientation tensor. The transfer functions used in the hidden and output layer are “Tansig” and “Purelin,” respectively. The Levenberg-Marquardt algorithm was used to train the ANN model, and ANN can converge with much fewer iterations. By changing the number of neurons in the hidden layer from 5 to 15, the ANN topology of 6-13-1 was determined according to the minimum mean square error between the targets and the outputs, indicating that there were six neurons in the input layer, thirteen neurons in the hidden layer, and one neuron in the output layer (Figure 4).
2.6. Genetic Algorithm (GA)
Genetic algorithm based on natural selection and survival of fitness is a global searching algorithm and is widely used in the fields of optimization, pattern recognition, robots, and prediction . Compared with other optimization methods, it has many advantages including being not easy to be trapped into the local minima, requiring little prior information about the searched objectives, and easy identifying of the optima in a complex search space.
The flow chart of genetic algorithm is shown in Figure 5. The major operations of GA are summarized as follows. (1) Selection: individuals are selected based on their fitness so that better individuals are given a higher chance of being chosen. (2) Crossover: exchange the information of the two parents to generate a new individual according to the crossover probability. (3) Mutation: randomly alter the information of each chromosome according to the mutation probability. After the predefined evolution generations or the resulting solution is satisfied, the genetic algorithm is stopped.
3. Results and Discussion
3.1. Short Fiber Circumferential Orientation Distribution in RWT
After the simulation is completed, a cross-section extracted from the middle of the model is used to observe the short fiber orientation distribution. As shown in Figure 6, the SFCO distribution in the RWT of the tube has a distinct layered structure. The outer layer has a small SFCO distribution with a component of orientation tensor about 0.15. The SFCO distribution in the inner layer near the water channel increases significantly. The component of orientation tensor in the region near the water channel is about 0.33, indicating that the short fibers tend to be freely oriented.
Figure 7 is a typical short fiber orientation distribution state of RWT in OWAIM. The RWT can be divided into two regions according to the characteristics of the short fiber orientation distribution. The short fibers in the region near the mold wall are arranged mainly along the axial direction, while these near the water channel are irregularly arranged. Figure 7(a) is a magnified view of the area selected from the inner layer, in which there are dense holes caused by the cryofractured process, indicating many short fibers oriented in the circumferential direction. Figure 7(b) is a magnified view of the area selected from the outer layer of the RWT, in which the existing short fibers have a regular arrangement in the axial direction. The number of holes in outer layer is significantly reduced, and many elongated vacancies parallel to the axial direction are left in the photograph, revealing the short fibers in outer layer mainly arranged in the axial direction and very few short fibers oriented in the circumferential direction. A conclusion can be inferred from Figure 7 that the value of A22 in the inner layer is larger than that in the outer layer. Compared with the result obtained in Figure 6, it can be found that the iARD-RPR model is well suited for the short fiber orientation prediction in OWAIM.
3.2. Orientation Mechanism Analysis of the Short Fibers
The process of OWAIM includes a melt filling stage and a high-pressure water filling stage. Figure 8 shows the components of orientation tensor in the position where the RWT is located during two filling stages. The melt filling stage of OWAIM is similar to that of conventional injection molding, in which the large shear stress and velocity gradient changes exist between the thin layers near the mold cavity. This makes the short fibers highly orient along the axial direction due to the shear action. There are very few short fibers aligned in the circumferential direction and A22 is small in the whole RWT after melt filling stage. After the delay time, influenced by the low temperature of the mold cavity, the temperature of the surface melt is lowered to form a high-viscosity layer, wherein the short fiber orientation is solidified. During the penetration process of high-pressure water, the water column is in a turbulent state because of the short injection time. The residual melt is squeezed and rubbed by the water column at the interface between the melt and the water column. In addition, the flowing melt driven by the water column stretches the residual melt. The short fibers in the layer near the water channel, which are originally oriented along the axial direction, readjust the orientation posture, so that A11 decreases and A22 increases.
3.3. Sensitivity Analysis of Process Parameters
Table 3 shows the values of A22 in the orthogonal experimental scheme. The results of range analysis are listed in Table 4, in which is the sum value of A22 of all level in each factor. According to the magnitude of the range R, the order of sensitivity of A22 to the six process parameters is determined as follows: water pressure, delay time, melt temperature, filling time, water temperature, and mold temperature. The optimal process parameters combination recommended by range analysis is A4B3C5D1E4F3, that is, melt filling time 3s, melt temperature 240°C, mold temperature 80°C, delay time 2s, water pressure 7.5MPa, and water temperature 30°C. From the analysis of variance (Table 5), it is concluded that melt temperature, delay time, and water pressure are the predominant factors with the F values of 7.782, 12.875, and 21.231, respectively, while filling time, mold temperature, and water temperature have no significant effects on the A22 with F values of 1.918, 0.938, and 1.061, respectively.
indicates being significant.
indicates being highly significant.
3.4. Modeling the Process Parameters Using ANN
Experimental data (Table 3), obtained from orthogonal experimental design, were divided randomly into three data sets. 19 of overall 25 points were used to train the ANN model, and the other points (3+3) were used to validate and test this model, respectively. The trained neural network was used to predict the values of A22 in the orthogonal experimental scheme. The predicted values shown in Figure 9, which are very consistent with the expected value, indicate that the neural network has established a nonlinear relationship between the process parameters and the A22.
3.5. Effect of Significant Factors on the SFCO
Figure 10 shows the effect of a single significant parameter on the A22 predicted by the ANN model. It indicates that the values of A22 increase as the melt temperature and water pressure increase, while A22 has a reverse relation with the delay time. As mentioned above, the interferences of high-pressure water column with the inner layer of RWT result in an increase in A22. The greater the melt temperature and water pressure are, the higher the disturbance to the inner layer is, causing an increase in A22. The longer delay time is, the larger thickness of solidified layer is, and the thinner thickness of melt is disturbed by the high-pressure water, leading to the decrease of the A22.
Figures 11–13 are the response surface contours indicating the interaction effects of two significant process parameters on A22. Figure 11 shows the interaction effect of melt temperature and delay time. The value of A22 varies within the range of , and maximum of A22 is obtained with delay time 2s and melt temperature 260°C. Figure 12 shows the interaction effect of delay time and water pressure. The value of A22 varies within the range of [0.185, 0.215], and maximum A22 appears with delay time 2s and water pressure 8MPa. As shown in Figure 13, the value of A22 varies within the range of [0.182, 0.204] for the interaction effect of melt temperature and water pressure, and the maximum A22 is got with melt temperature 260°C and water pressure 8MPa.
3.6. Optimization and Verification of the Circumferential Orientation
Taking the trained ANN model as the fitness function, the optimization of A22 was carried out using the solver of “ga, genetic algorithm” in the MATLAB R2015b optimization toolbox, with the principle “the larger, the better.” The parameters of GA were set as the follows: population size 50, initial range , crossover Fraction 0.8, mutation function 0.01, and generation 100. The evolution of optimization process is recorded in Figure 14. After 100 generations, the optimized A22 is found to be 0.234. The corresponding process parameters are filling time 3.5s, melt temperature 260°C, mold temperature 80°C, delay time 2s, water pressure 8MPa, and water temperature 36°C, respectively. This value of A22 is greater than any set of experimental results in the orthogonal experimental scheme. The simulation experiments using the optimal parameters combination recommended by GA and range analysis are conducted by Moldex3D, and the obtained A22 are 0.233 and 0.228, respectively. Results show that the SFCO distribution is improved by two methods, and the combination of ANN and GA is better than that of range analysis in orthogonal experimental scheme.
In this research, the short fibers orientation distribution of RWT in OWAIM was simulated using the iARD-RPR model. Compared with the SEM micrographs, the simulation results indicated that this model was suitable for OWAIM. The simulation results showed that the value of A22 in the outer layer is small and that in the inner layer increases significantly. In addition, the simulations indicated that the penetration of water column readjusted the orientation state of short fibers, increasing the circumferential orientation of short fibers in the inner layer of RWT. Through the range analysis and variances analysis, it was found that melt temperature, delay time, and water pressure were significant factors. The values of A22 increase as the melt temperature and water pressure increase, while decreasing with the delay time rise. The maximum A22 was 0.234 in the optimization with the combination of ANN and GA. Compared with range analysis, the combination of ANN and GA resulted in a better optimization result. This study will help to further understand the orientation mechanism of short fibers in OWAIM.
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
The work in this paper is funded by the National Natural Science Foundation of China (Grant No. 21664002, Grant No. 51563010) and Jiangxi Provincial Key Technology R&D Program (Grant No. 20161BBE50073). The authors would like to express their sincere gratitude to those who made comments on the changes proposed in this article.
- J. Zhang, P. Zhao, Y. Zhao, J. Huang, N. Xia, and J. Fu, “On-line measurement of cavity pressure during injection molding via ultrasonic investigation of tie bar,” Sensors and Actuators A: Physical, vol. 285, pp. 118–126, 2019.
- P. Zhao, W. M. Yang, X. M. Wang, J. G. Li, B. Yan, and J. Z. Fu, “A novel method for predicting degrees of crystallinity in injection molding during packing stage,” Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture, vol. 233, no. 1, pp. 204–214, 2019.
- T. Q. Kuang, B. P. Xu, G. F. Zhou, and L. S. Turng, “Numerical simulation on residual thickness of pipes with curved sections in water-assisted co-injection molding,” Journal of Applied Polymer Science, vol. 132, no. 34, pp. 1–10, 2015.
- N. Xia, P. Zhao, T. Kuang, Y. Zhao, J. Zhang, and J. Fu, “Nondestructive measurement of layer thickness in water-assisted coinjection-molded product by ultrasonic technology,” Journal of Applied Polymer Science, vol. 135, no. 33, pp. 2–9, 2018.
- T.-Q. Kuang, K. Zhou, L.-X. Wu, G.-F. Zhou, and L.-S. Turng, “Experimental study on the penetration interfaces of pipes with different cross-sections in overflow water-assisted coinjection molding,” Journal of Applied Polymer Science, vol. 133, no. 1, pp. 1–8, 2016.
- S.-J. Liu and S.-P. Lin, “Factors affecting the formation of fingering in water-assisted injection-molded thermoplastics,” Advances in Polymer Technology, vol. 25, no. 2, pp. 98–108, 2006.
- Q. Jiang, H. Liu, Q. Xiao, S. Chou, A. Xiong, and H. Nie, “Three-dimensional numerical simulation of total warpage deformation for short-glass-fiber-reinforced polypropylene composite injection-molded parts using coupled FEM,” Journal of Polymer Engineering, vol. 38, no. 5, pp. 493–502, 2017.
- S. Liu and Y. Chen, “The manufacturing of thermoplastic composite parts by water-assisted injection-molding technology,” Composites Part A: Applied Science and Manufacturing, vol. 35, no. 2, pp. 171–180, 2004.
- S. J. Liu, M. J. Lin, and Y. C. Wu, “An experimental study of the water-assisted injection molding of glass fiber filled poly-butylene-terephthalate (PBT) composites,” Composites Science and Technology, vol. 67, no. 7-8, pp. 1415–1424, 2007.
- H. Huang, R. Zhou, and C. Yang, “Fiber orientation propelled by high-pressure water penetration in water-assisted injection molded fiber-reinforced thermoplastics part,” Journal of Composite Materials, vol. 47, no. 2, pp. 183–190, 2012.
- I. Fiebig and V. Schoeppner, “Factors influencing the fiber orientation in welding of fiber-reinforced thermoplastics,” Welding in the World, vol. 62, no. 5, pp. 997–1012, 2018.
- R. Taktak, N. Guermazi, and T. K. Kallel, “Effect of E-Glass fibre and ply orientation on the mechanical behaviour of FRP composites used for pressure pipe,” The International Journal of Advanced Manufacturing Technology, vol. 92, no. 5-8, pp. 1741–1749, 2017.
- S. Mortazavian and A. Fatemi, “Effects of fiber orientation and anisotropy on tensile strength and elastic modulus of short fiber reinforced polymer composites,” Composites Part B: Engineering, vol. 72, pp. 116–129, 2015.
- U. Delpy and G. Fischer, “Effect of mold‐filling conditions on fiber distribution in injection‐molded disks and on the mechanical properties of such disks,” Advances in Polymer Technology, vol. 5, no. 1, pp. 19–26, 1985.
- S. Lee, D. Yang, J. Ko, and J. Youn, “Effect of compressibility on flow field and fiber orientation during the filling stage of injection molding,” Journal of Materials Processing Technology, vol. 70, no. 1-3, pp. 83–92, 1997.
- K. Li, S. Yan, W. Pan, and G. Zhao, “Optimization of fiber-orientation distribution in fiber-reinforced composite injection molding by Taguchi, back propagation neural network, and genetic algorithm–particle swarm optimization,” Advances in Mechanical Engineering, vol. 9, no. 9, pp. 1–11, 2017.
- H. Tseng, R. Chang, and C. Hsu, “Numerical prediction of fiber orientation and mechanical performance for short/long glass and carbon fiber-reinforced composites,” Composites Science and Technology, vol. 144, pp. 51–56, 2017.
- H. Tseng, R. Chang, and C. Hsu, “Effect of the packing stage on fiber orientation for injection molding simulation of fiber-reinforced composites,” Journal of Thermoplastic Composite Materials, vol. 31, no. 9, pp. 1204–1218, 2017.
- H. Tseng, R. Chang, and C. Hsu, “An objective tensor to predict anisotropic fiber orientation in concentrated suspensions,” Journal of Rheology, vol. 60, no. 2, pp. 215–224, 2016.
- S. G. Advani and C. L. Tucker, “The use of tensors to describe and predict fiber orientation in short fiber composites,” Journal of Rheology, vol. 31, no. 8, pp. 751–784, 1987.
- B. O. Agboola, D. A. Jack, and S. Montgomery-Smith, “Effectiveness of recent fiber-interaction diffusion models for orientation and the part stiffness predictions in injection molded short-fiber reinforced composites,” Composites Part A: Applied Science and Manufacturing, vol. 43, no. 11, pp. 1959–1970, 2012.
- V. Kůrková, “Kolmogorov's theorem and multilayer neural networks,” Neural Networks, vol. 5, no. 3, pp. 501–506, 1992.
- K. Li, S. Yan, W. Pan, and G. Zhao, “Warpage optimization of fiber-reinforced composite injection molding by combining back propagation neural network and genetic algorithm,” The International Journal of Advanced Manufacturing Technology, vol. 90, no. 1-4, pp. 963–970, 2017.
Copyright © 2019 Haiying Zhou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.