Numerical Infeasibilities of Nanofibrous Mats Process Design

Abstract

A new computer-aided method to design electrospun, nanofibrous mats was implemented and tested. In this work, the standard nonlinear algebraic model led to the terminal fiber diameter ${\mathrm{F}}_{\mathrm{D}}$ being examined in detail. The analysis was performed in terms of numerical feasibility. The study specified the limit value of the axial length scale, parameter $\chi$, that determined valid solutions. The presented approach has vast practical potential (i.e., biomedical applications, air/water purification systems, fire protection and solar industries).

1. Introduction

Electrospinning (ES) has been one of the most dynamically developing technologies over past decades [1]. A wide range of natural and synthetic nanofibrous mats have been successfully electrospun [2,3]. ES contributes to the formation of polymeric fibers with a nano-scale diameter, that can mimic extracellular matrix [4]. Thanks to this feature, ES may be considered as a technology suitable for biomedical applications. Dettin and co-authors [5] presented in their work a problem associated with the mechanism of cell–biomaterial adhesion. For this purpose, they characterized electrospun poly-ε-caprolactone (PCL) scaffolds with an increasing concentration of self-assembling peptides (SAPs). Bazzolo and co-authors proposed the use of polycaprolactone electrospun structures in preclinical studies in relation to effective breast cancer therapy [6]. Bombin et al. drew attention to the latest trends in using ES to develop structures applicable in wound healing [7]. Some recent studies have focused on the use of electrospun materials in the filtration/purification processes [8], the flame retardant sector [9], and the solar industry [10].

Most studies have been content with an experimental method of material design. Celebioglu and Uyar conducted the experimental optimization of the electrospinning process of gamma-cyclodextrin/poly(ethylene oxide) aqueous solutions [11]. As a result, they obtained uniform fibers capable of removing aniline. Extensive experimental optimization has also been explored in the work [12]. The authors established the optimum conditions in which to fabricate hydrolyzed Polymer of Intrinsic Microporosity (PIM-1) structures of a high affinity in relation to cationic dyes (Methylene Blue). Wang and Ziegler [13] proposed optimized processing parameters to produce fibers from starch-pullulan aqueous dispersions.

Traditional methods used by engineers and scientists have been changed due to the appearance of powerful, high-performance computers. This general trend is reflected in numerical simulations, that together with experiment and theory, have become the central aspects of scientific discovery [14]. Today, a high-quality system design proceeds with the use of a computer. This approach may contribute to the growth of innovative solutions, new materials, and advanced technologies.

In one study [15], the Artificial Neural Network model was used to predict the terminal diameter of ferrofluid-polyvinyl alcohol fibers. The model took into account variables such as voltage, flow rate, spinning distance and the (collector) rotation speed. Šimko and Lukáš [16] presented the mathematical modeling strategy based on works [17,18]; the simulation parameters were employed for a 6 wt % aqueous solution of poly (ethylene oxide).

Akkoyun and Öktem [19] determined the fiber diameter of polyamide 6 based on numerical and experimental investigations. The proposed approach utilized Finite Extensible Non-linear Elastic–Chilcott and Rallison (FENE-CR) model; the authors examined the influence of initial polymer contribution. It is worth noting that a considered numerical approach allowed both the viscoelastic material behavior and fiber solidification to be captured.

Mathematical models may significantly accelerate the design cycle and allow us to investigate the system’s behavior. It is important to construct a model that reflects the considered tasks with appropriate accuracy. Some questions related to the solution procedure arise at this point.

Therefore, to obtain numerical data that represent system behavior in different conditions and for a wide range of selected model parameters, it is desirable to perform numerous computer simulations. Then, based on the results obtained it is possible to focus on the detailed optimization of the received approximated solution. Nevertheless, in some instances for the chosen set of parameters the model cannot be solved. In this situation, we usually obtain the information about a calculation error. It is worth considering whether computer simulations, which in engineering practice may terminate with a numerical error (especially when we consider new and difficult issues), are in fact worthless, or whether they might be a source of useful information.

The aim of this work is to discuss the computational possibilities of the commonly known and used nonlinear algebraic model [20,21,22,23]. Most studies have relied on terminal fiber diameter prediction or model validation. The current work examines the validity of the solution returned by the model in terms of numerical feasibility. So far, no studies have been presented in this field.

2. Theoretical Support

The control of the fiber diameter is an important element of the electrospinning process [24]. There are several attempts to model the electrically charged jet behavior described in the literature [25,26]. A basic electrospinning setup consists of an infusion pump with a polymer solution, a high voltage power supply, and a collector [27]. The process involves applying a high voltage to electrodes, which are a syringe needle and the collector [28]. In other words, a strong electrical field is applied to the polymer solution. This operation induces a charge on the surface of the fluid. As the electrical field strength increases, a conical shape called the Taylor cone begins to be created. A straight jet is then formed, which undergoes a bending instability phenomenon, as it travels towards the collector [29,30]. This process causes the fiber diameter to decrease (by orders of magnitude) relative to the diameter of the syringe needle [31].

Basically, there are four zones in the ES: the Taylor cone formation (T); the stable part of the jet (J); the unstable, whipping part (W); and the area where fibers are stopped and collected (C) [32]. A schematic view of the basic ES equipment with the four separate states is depicted in Figure 1.

Fridrikh and co-authors [20] in 2003 presented the model of a charged Newtonian fluid jet applicable to the whipping instability; the model allowed the prediction of a terminal fiber diameter with the known values of an electric current, flow rate and surface tension. The model neglected elastic effects and solvent evaporation. The model of Fridrikh et al. [20] has been widely discussed and modified [33,34,35,36,37,38,39,40,41,42].

Gadkari, in [40], has noticed that the model presented in [20], due to the omission of important dependencies (i.e., viscosity, solution evaporation, polymer contribution), leads to a significant exaggeration in the importance of the surface tension. Ismail et al. [21,22] proposed the combined version of the stable and unstable jet based on works [43,44] and [20]. The model enabled the prediction of the final fiber diameter regarding the volumetric flow rate, the applied electrical field, and the polymer concentration. Moreover, the authors of the work [38] described the approach where the exponential model was proposed. The model related the stable jet length and the terminal fiber radius. The authors showed that for fully charged fibers, the scaling exponent was “−2/3” and that this becomes smaller for partially charged fibers.

The parameter $\chi$ is a very interesting part of the work described above. The literature definitions of $\chi$ can be presented as follows: “${\chi }^{-1}$ is the local aspect ratio, which is assumed to be small” [43], “$\chi ~R/h$ is the dimensionless wavelength of the instability responsible for the normal displacements” [20], $\chi$ is “ a characteristic length parameter (…) influenced by the first stage of electrospinning, which is the stable jet” [21], and ”$\chi$ is the aspect ratio of the jet” [38].

Thus, we can also state:

$$\chi ~\frac{\lambda }{r_{sj,end}},$$

(1)

where $\lambda$ denotes the wavelength and $r_{sj,end}$ is the critical radius of the stable jet (at the end of the straight jet region), and

$$\lambda =\frac{2\pi }{k},$$

(2)

that can be estimated by the radial instability wavelength from the fastest growing mode [38].

The definition of ${\chi }_{bi}$ described as “the ratio of the characteristic axial length scale” can also be found in the work of [21]. It is also worth exploring the $\chi$ range and determining which values lead to nanofibers at certain diameters. Fridrikh et al. [20] in all computations used $\chi =100$; the authors assumed also that the $\chi$ value between $10 \mathrm{and} 1000$ was reasonable. The estimation method of the $\chi$ parameter and the ${\mathrm{F}}_{\mathrm{D}}\left(\chi \right)$ relation was presented, analyzed and discussed in detail in [23].

3. Problem Statement

The aim of this study is a further analysis of the model presented in the work [23]; the model is a classic example, commonly known and widely used in the literature [20,21,22]. The conventional formulation makes it possible to obtain an explicit ${\mathrm{F}}_{\mathrm{D}}$ for given model parameters. However, it is quite surprising that not every set of such parameters results in an expected value. The relationships between the values of ${\mathrm{F}}_{\mathrm{D}}$ and $\chi$ presented in the prior research [23] encourage the further examination of these connections in terms of numerical simulations. It is worth emphasizing that the goal is not to correct or improve the model, but to thoroughly investigate the situations in which its application may result in an acceptable solution.

The relevant model is formulated as follows:

$${\mathrm{F}}_{\mathrm{D}}\left(\mathsf{\chi }\right)={\mathrm{C}}^{0.5}{\left(\mathsf{\gamma }\mathsf{\epsilon }\frac{{\mathrm{Q}}^2}{{\mathrm{I}}^2}\frac{2}{\mathsf{\pi }\left(2\ln \mathsf{\chi }-3\right)}\right)}^{1/3},$$

(3)

subject to

$${\mathrm{F}}_{\mathrm{D}}>0$$

(4)

$$\mathsf{\gamma },\mathsf{\epsilon },\mathrm{C},\mathrm{Q},\mathrm{I},\mathrm{E},\mathrm{K}>0$$

(5)

$$\chi >0$$

(6)

$$2\ln \mathsf{\chi }-3>0$$

(7)

and

$$\mathrm{I}~{\mathrm{EQ}}^{0.5}{\mathrm{K}}^{0.4}$$

(8)

with $\mathsf{\gamma },\mathsf{\epsilon },\mathrm{C},\mathrm{Q},\mathrm{I},\mathrm{E},\mathrm{K}\in ℛ$ and scalar real-valued function ${\mathrm{F}}_{\mathrm{D}}$. Moreover, from the previous work [23] we know that

$$4.5\le \chi \le 6.70\times {10}^5$$

(9)

where ${\mathrm{F}}_{\mathrm{D}}$ determines the final diameter of the fibers collected from the collector; $\mathrm{C},\mathrm{Q},\mathsf{\Delta }\mathrm{V},\mathrm{z}$ are the concentration of the polymer solution, flow rate, applied voltage and spinning distance; $\mathsf{\gamma },\mathsf{\epsilon },\mathsf{\chi }$ are surface tension, permittivity of the outside medium and the axial length scale; $\mathrm{I}$ determines the electric current, $\mathrm{K}$ denotes solution conductivity and $\mathrm{E}$ determines the electrical field.

The considered limitations are dictated by the process nature, Equations (4) and (5), mathematical principles (7), or both (6).

4. Design Challenge

The design of nanofibrous mats can be achieved by using an experimental method. This approach is concerned with the literature studies, detailed statements about research objectives, data collection, and data analysis. Numerical simulations enable a reduction in the cost of experiments, as well as design time. Unfortunately, even simple mathematical expressions may lead to undesired results.

Let us consider the motivating example given by Beykal et al. [45]

$$y\left(t\right)=\frac{y_0}{1-y_0·t}$$

(10)

subject to

$$t>0$$

(11)

$$y_0>0$$

(12)

$$y\left(t\right)>0$$

(13)

where $y_0=y\left(0\right)$ denotes the initial value, $t$ is the simulation time and $t\in \left[0,3\right]; y_0\in \left[0,3\right]$.

The overall problem (10) revolves around the validity of the solution; the solution is valid when all constraints (11)–(13) are satisfied. The sample set of solutions found for this case is presented in Figure 2.

The authors of the motivating example [45] noted that although it is not difficult to derive limitations for the Equation (10), the valid solution may not be trivial (or may not even exist) for various complex engineering issues. By following the same procedure, we can find the sample set of solutions and designate feasibility and infeasibility regions for $FD\left(\chi \right)$ dependency, presented by Equations (3)–(9); see Figure 3.

5. Numerical Procedure

Figure 4 outlines the general procedure of the numerical investigations executed in this work. The analysis was based on the nonlinear algebraic model of the ES (Equation (3)). The first step was to assume preliminary suppositions and constraints (Equations (4)−(9)). Then, feasible and infeasible regions were determined; sampling points that satisfied these constraints were depicted as blue circles (Figure 3). Then, the detailed analysis imposed on the $\mathsf{\chi }$ parameter was performed. The obtained effects were presented in Section 6 and Section 7. Concluding remarks are given in Section 8. All calculations were performed in the MATLAB Environment. Simulation parameters are listed in

Table 1. Simulation parameters.

Solution Concentration (wt %)	Flow Rate (mL/h)	Tip-to-Collector Distance (cm)	Applied Voltage (kV)
8	1.6	15	16
8	1.6	15	25
8	1.6	15	27

.

6. Simulation Results

Five graphs, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, were obtained as a result of numerical investigations. Figure 5 presents the detailed values of the $\chi$ parameter that lead to infeasible and feasible regions. Figure 6, Figure 7, Figure 8 and Figure 9 present the model results for the wide (Figure 6), narrow (Figure 7 and Figure 8), and very tight (Figure 9) ranges of $\chi$.

7. Discussion

The target issue is the distinction between failed and successful simulations. The simulation is considered to have failed when the examined set of points do not satisfy at least one of the imposed constraints or if they do not belong to the function domain.

Let us reduce Equation (3) to the form

$${\mathrm{F}}_{\mathrm{D}}=\mathrm{A}\sqrt[ 3]{\frac{1}{2\ln \mathsf{\chi }-3}},$$

(14)

where $\mathrm{A}={\mathrm{C}}^{0.5}{\left(\mathsf{\gamma }\mathsf{\epsilon }\frac{{\mathrm{Q}}^2}{{\mathrm{I}}^2}\frac{2}{\mathsf{\pi }}\right)}^{1/3} \in ℛ$, for a given set of process parameters.

Figure 5 depicts the detailed value of the $\chi$ parameter ($\chi *\approx 4.48$) that leads to the feasible region. More accurate and time-consuming numerical computations resulted with ${\chi }^{*}\approx 4481.69\times {10}^{-3}$. The infeasible region was marked in grey; these results are beyond expectations.

Therefore, in order to obtain successful simulations, we should take $\chi \ge \chi *$. Otherwise, our simulations will have failed. The $\chi <\chi *$causes a denominator (Equation (14)) of less than zero. As a result we have a cube root of a negative number. This leads to

• Negative values of the ${\mathrm{F}}_{\mathrm{D}}$ (real cube root), where $\chi <\chi *$violates the Equation (4);
• To the situation, where the ${\mathrm{F}}_{\mathrm{D}}$ value is expressed as the complex number (complex cube root), where $\chi <\chi *$violates the Equation (7).

In both cases the obtained result is undesirable. In practice, the obtained solutions were expressed in the form of complex numbers due to the power function used in all calculations (see [46,47]).

Figure 6 shows the model results when the fiber diameter approaches 500 nm and $\chi \in \left[\chi *, 7.00\times {10}^5\right]$. Figure 7 depicts the area where $\chi \in \left[\chi *,50\right]$. Here, one can observe the decreasing function ${\mathrm{F}}_{\mathrm{D}}\left(\chi \right)$. Let us find

$$\underset{\chi \to +\infty }{\lim }{\mathrm{F}}_{\mathrm{D}}\left(\chi \right)=0$$

(15)

Thus, we can see that as the $\chi$ increases the ${\mathrm{F}}_{\mathrm{D}}\left(\chi \right)$ goes to zero. This situation is more clearly illustrated in Figure 8. Here, we can see the initial significant change in the ${\mathrm{F}}_{\mathrm{D}}$ value (sharp slope area). ${\mathrm{F}}_{\mathrm{D}}$ varies at a slower rate. It is worth noting that in the sharp slope area even a small change in the $\chi$ value may lead to incorrect, distorted results, whereas approaching the desired solution may be achieved by inconsiderable change of the $\chi$ parameter, especially for the smaller ${\mathrm{F}}_{\mathrm{D}}$ (i.e., when ${\mathrm{F}}_{\mathrm{D}}=100$ nm or less).

It should be also added that

$$\underset{\chi \to {\chi }^{*+}}{\lim }{\mathrm{F}}_{\mathrm{D}}\left(\chi \right)=+\infty$$

(16)

This is particularly well demonstrated by Figure 9, where ${\mathrm{F}}_{\mathrm{D}}$approaches $3.25$ mm ($3.25\times {10}^6$) nm. According to Equation (16), the ${\mathrm{F}}_{\mathrm{D}}$ goes to $+\infty$, hence ${\mathrm{F}}_{\mathrm{D}}$ may reach high values. However, ${\mathrm{F}}_{\mathrm{D}}=3.25$ mm is against the process nature.

Moreover, we can also write, that ${\mathrm{F}}_{\mathrm{D}}\left(\chi \right)=\mathrm{A}\sqrt[3]{\frac{1}{2\ln \mathsf{\chi }-3}}$ and

$$\widetilde{{\mathrm{F}}_{\mathrm{D}}}\left(\chi \right)=\frac{a}{\mathsf{\chi }-\chi *}$$

(17)

with $A,a\in ℛ$, characterized by similar asymptotical behavior.

One should also emphasize that the sampling frequency affects the results. The sampling should be appropriately adjusted to the process nature and to the model that describes it. The sampling frequency for the simulations depicted above was: $1.00$ (Figure 6); $1.00\times {10}^{-5}$ (Figure 5, Figure 7 and Figure 8); $1.00\times {10}^{-11}$ (Figure 9). Dense sampling needs to be carried out for uncertain, suspicious regions. Nevertheless, dense sampling may significantly affect the calculation time—for example, the simulation time took about 21 minutes for results illustrated in Figure 9 for Intel(R) Core(TM) i7-8665U CPU 1.90 GHz with 16.0 GB RAM and the Windows 10 Pro operating system.

8. Concluding Remarks

This study has investigated the nonlinear algebraic model in terms of numerical feasibility. The model is commonly known and used by numerous research groups; the model relates the ES process parameters with the final fiber diameter. The identification of the feasible and infeasible regions led to the determination of the proper set of points that satisfied all imposed constraints; the work was particularly focused on constraints imposed on the axial length scale: χ parameter. The detailed value of the χ parameter, that ensured failure-free computations, was $\chi *\approx 4.48$. A minimal change in the $\chi *$ value ($\mathsf{\chi }$ less than $\chi *$) went to negative or complex numbers. In both cases, the results were undesirable. The analysis enabled the achievement of valid solutions and assessment of the mathematical nature of the model. Former studies have almost exclusively focused on the selected areas of the model applicability without discussing the entire spectrum in which they can be used. The objective was not to improve the model, but to analyze, in-depth, various results of numerical computations. It is worth noting that even simulations that resulted in unexpected values were effectively used to estimate the feasible region of the $\mathsf{\chi }$ parameter.