Mesoscale Simulations of Structure Formation in Polyacrylonitrile Nascent Fibers Induced by Binary Solvent Mixture

Polyacrylonitrile (PAN) is widely used as a raw material for the production of high-modulus carbon fibers, the internal structure of which is directly affected by the spinning of the precursor. Although PAN fibers have been studied for a long time, the formation of their internal structure has not been sufficiently investigated theoretically. This is due to the large number of stages in the process and the parameters controlling them. In this study, we present a mesoscale model describing the evolution of nascent PAN fibers during the coagulation. It is constructed within the framework of a mesoscale dynamic density functional theory. We use the model to study the influence of a combined solvent of dimethyl sulfoxide (DMSO, a good solvent) and water (a non-solvent) on the microstructure of the fibers. A porous structure of PAN is formed as a result of the microphase separation of the polymer and the residual combined solvent at a high water content in the system. The model shows that one of the possible ways to obtain the homogeneous fiber structure is to slow down the coagulation by increasing the amount of good solvent in the system. This result is in agreement with the existing experimental data and confirms the efficiency of the presented model.


Introduction
Polyacrylonitrile (PAN) has a low density, thermal stability, high strength, elastic modulus, and chemical stability. After functionalization with different chemical groups, it can form stable materials with different structures. These properties have established PAN as essential material in chemical and textile industry, for foods packing materials, electronics, and cosmetics [1,2]. PAN fibers are also used for outdoor awnings and clothing, and they play a key role in high-tech applications such as air/water filters [3,4]. They are also utilized as a precursor for the manufacturing of carbon fibers (CFs) due to its good mechanical properties, high carbon content, and ability to undergo cyclization reactions at high temperatures [2,[5][6][7][8][9][10][11][12][13][14]. CFs are widely used in high-tech applications where heat resistance, chemical inertness, and high tensile strength are required [5][6][7][8][15][16][17]. Their formation is accompanied by complex physical and chemical transformations in the precursor, triggered by the influence of solvent quality, mechanical stress, and temperature. The supramolecular structures (nanofibrils, networks, crystals, etc.) formed during the multiphase transition have a decisive influence on the properties of the final product. It is believed that the properties of CFs can be controlled by varying the formation conditions of precursor fibers [6,9,14,16].
One of the most effective ways to form the desired structure of a polymeric material is to perform a controlled phase transformation in polymer solutions [18][19][20]. Wet and dry-wet spinning processes for the production of PAN fibers are well suited for this purpose [5,7,12,15,16,21]. In these technologies, the spinning solution (polymer solution on the methodology of simulations. Finally, we summarize the results of our work in the Section 5. Additional information can be found in the Supplementary Materials file.

Results
First, we performed calculations for systems with different PAN volume fractions C P of 20 vol%, 60 vol%, and 80 vol% and different water contents in the composition of the combined solvent f w = 0-1 (f w = C W /(C D + C W ), where C D and C W are DMSO and water volume fractions). Note that the following abbreviations are used in the text and Figures: P-PAN, D-dimethyl sulfoxide (DMSO), W-water. The rate of microphase separation increases with increasing water content in the combined solvent (see Figure 1A). Due to the large difference in the values of the solubility parameters δ P and δ W at a low percentage of PAN (C P = 20 vol%), the polymer forms isolated domains in the final state (see Figure 1B, snapshot 1). As C P increases, the observed picture is reversed. In this case, water begins to form isolated domains embedded in the polymer matrix. At their boundary, a transition layer of DMSO is formed (see Figure 1B, snapshot 3). At the same time, DMSO (in all cases) remains in the domains formed by the polymer, because its solubility parameter is close to that of PAN.
The rest of the article is organized as follows. The results obtained are discussed in Section 2. The limits of our model and the portability of the results to other systems are discussed in the Section 3. Section 4 gives brief information on the modeling method used and a description of the model construction, its parameterization and provides information on the methodology of simulations. Finally, we summarize the results of our work in the Section 5. Additional information can be found in the Supplementary Materials file.

Results
First, we performed calculations for systems with different PAN volume fractions СP of 20 vol%, 60 vol%, and 80 vol% and different water contents in the composition of the combined solvent fw = 0-1 (fw = СW/(СD + СW), where СD and СW are DMSO and water volume fractions). Note that the following abbreviations are used in the text and Figures: P-PAN, D-dimethyl sulfoxide (DMSO), W-water. The rate of microphase separation increases with increasing water content in the combined solvent (see Figure 1A). Due to the large difference in the values of the solubility parameters δP and δW at a low percentage of PAN (СP = 20 vol%), the polymer forms isolated domains in the final state (see Figure 1B, snapshot 1). As СP increases, the observed picture is reversed. In this case, water begins to form isolated domains embedded in the polymer matrix. At their boundary, a transition layer of DMSO is formed (see Figure 1B, snapshot 3). At the same time, DMSO (in all cases) remains in the domains formed by the polymer, because its solubility parameter is close to that of PAN.   (1) C P = 20 vol%, f w = 0.77; (2) C P = 60 vol%; f w = 0.77; (3) C P = 80 vol%, f w = 0.7. Red color corresponds to the PAN matrix (ρ P > 0.8), green color corresponds to DMSO (ρ D > 0.02), and blue color corresponds to water (ρ W > 0.8). The following abbreviations are used here and in the Figures below: P-PAN, D-DMSO, W-water. Figure 2 shows an instance of the distribution profiles N(ρ α )/N total characterizing the volume fraction of the simulation cell regions where the density of the components of the system (P, D, and W) is equal to certain values of densities fields ρ α (α = P, D, W) for the system with C P = 80 vol%, f w = 0.5. See Section 4 for a description of ρ α . As can be seen, the distribution N(ρ α ,t)/N total for PAN has a clear maximum at ρ P ≈ 1.02. The figure also shows that DMSO does not form regions of high density. system with CP = 80 vol%, fw = 0.5. See Section 4 for a description of ρα. As can be seen, the distribution N(ρα,t)/Ntotal for PAN has a clear maximum at ρP ≈ 1.02. The figure also shows that DMSO does not form regions of high density.  Figure 3A shows that the DMSO:water ratio has a strong influence on the rate (tps decrease with increasing water content) and the result of microphase separation in the system. By visualizing the distribution of DMSO and water in the volume of the simulation cell, it is possible to extract information about the dependence of the morphology of the water domains on the change in the water fraction in the system. Figure 3B (snapshot 1) shows that PAN and solvent form a homogeneous structure at a low water content (region I in Figure 3A). As the water content increases (fw > 0.15), PAN forms a separate phase and water initially forms spherical domains (region II in Figure 3A, snapshot 2 in Figure  3B). Then, with a further increase in water content, the water domains start to coalesce and form elongated dumbbell-shaped structures (region III in Figure 3A, snapshot 3 in Figure 3B). Finally, at a certain threshold (fw ≈ 0.8), the domains can coalesce to form a porous structure that percolates throughout the simulated cell (region IV in Figure 3A, snapshot 4 in Figure 3B).
We have also already discussed above that the DDFT model shows the formation of a "blue cheese" structure that occurs during the coagulation of the spinning solution at a high water content starting from about fw ≈ 0.5 (see Figures 1B, 3B and 4C). It is well known that the solvent-filled voids, embedded in the polymer matrix, can form pores during the fiber-drying step. At the same time, regions with the highest PAN density can act as a nucleus for the growth of the crystalline phase in the next stages of the process of drawing and heat treatment of the fiber.
Since the high strength of carbon fibers is limited by internal morphological defects, it is desirable that the precursor PAN fibers also have a homogeneous structure [61]. In order to assess the possible influence of the conditions of the coagulation process on the size of the high-density PAN regions, we have analyzed the final states of the simulated systems at different water contents and fixed volume fraction of PAN (СP = 80 vol%). It should be noted that in our DDFT model, we do not deal directly with the packing of polymer chains, but only with the densities of their spatial distribution. By high-density domains, we mean domains where the field density ρP > 1.03. An example of the visualization of domains with the highest PAN density and their fraction in the volume of the simulation cell is shown in Figures S2 and S3.  Figure 3A shows that the DMSO:water ratio has a strong influence on the rate (t ps decrease with increasing water content) and the result of microphase separation in the system. By visualizing the distribution of DMSO and water in the volume of the simulation cell, it is possible to extract information about the dependence of the morphology of the water domains on the change in the water fraction in the system. Figure 3B (snapshot 1) shows that PAN and solvent form a homogeneous structure at a low water content (region I in Figure 3A). As the water content increases (f w > 0.15), PAN forms a separate phase and water initially forms spherical domains (region II in Figure 3A, snapshot 2 in Figure 3B). Then, with a further increase in water content, the water domains start to coalesce and form elongated dumbbell-shaped structures (region III in Figure 3A, snapshot 3 in Figure 3B). Finally, at a certain threshold (f w ≈ 0.8), the domains can coalesce to form a porous structure that percolates throughout the simulated cell (region IV in Figure 3A, snapshot 4 in Figure 3B).  The pores present in the precursor fibers are usually considered defects that promote fiber destruction under high loads in the next drawing stages of PAN production [6,[48][49][50]. Thus, the results obtained confirm that one of the possible ways to reduce porosity is to slow coagulation by increasing the amount of good solvent in the system [5,[49][50][51][52][53][54][55][56][57][58][59][60].
We have also already discussed above that the DDFT model shows the formation of a "blue cheese" structure that occurs during the coagulation of the spinning solution at a high water content starting from about f w ≈ 0.5 (see Figures 1B, 3B and 4C). It is well known that the solvent-filled voids, embedded in the polymer matrix, can form pores during the fiber-drying step. At the same time, regions with the highest PAN density can act as a nucleus for the growth of the crystalline phase in the next stages of the process of drawing and heat treatment of the fiber.   Figure 5 shows the dependence of the average radius of the high-density domains in PAN on the water content of the system, on which four regions can be identified. At a low water content (region I in Figure 5), there are no dense domains. As the water content increases, PAN begins to form the dense domains, but their characteristics are not stable (region II). This leads to a large error in estimating the size of such regions in samples of the simulated system. With a further increase in the water content, the dense domains stabilize and their average radius increases rapidly (region III). The largest domains of Since the high strength of carbon fibers is limited by internal morphological defects, it is desirable that the precursor PAN fibers also have a homogeneous structure [61]. In order to assess the possible influence of the conditions of the coagulation process on the size of the high-density PAN regions, we have analyzed the final states of the simulated systems at different water contents and fixed volume fraction of PAN (C P = 80 vol%). It should be noted that in our DDFT model, we do not deal directly with the packing of polymer chains, but only with the densities of their spatial distribution. By high-density domains, we mean domains where the field density ρ P > 1.03. An example of the visualization of domains with the highest PAN density and their fraction in the volume of the simulation cell is shown in Figures S2 and S3. Figure 5 shows the dependence of the average radius of the high-density domains in PAN on the water content of the system, on which four regions can be identified. At a low water content (region I in Figure 5), there are no dense domains. As the water content increases, PAN begins to form the dense domains, but their characteristics are not stable (region II). This leads to a large error in estimating the size of such regions in samples of the simulated system. With a further increase in the water content, the dense domains stabilize and their average radius increases rapidly (region III). The largest domains of high density are formed in the polymer matrix at the highest water content in the system (region IV). The positions of the regions I-IV in Figure 5 correlate with the positions of the corresponding regions in Figure 3. At the same time, region I (low water content) corresponds to a homogeneous distribution of PAN, DMSO, and water density fields, while region IV corresponds to states where the pores formed by the water form a percolating structure.
bath. As a result, a core-shell structure is formed (a loose core with capillaries and voids filled with a solvent and a dense surface). The presence of a dense surface layer with numerous nanopores can also cause the fiber cross-section to deviate from a round shape. This can occur as a result of a non-isotropic distribution of nanopores relative to the symmetry axis of the fiber. Thus, our model clearly shows that the formation of precursor fibers with uniform structure is achieved at the high content of the good solvent in the coagulation bath [62,63]. In this way, the behavior of our model agrees with the conclusions of experimental studies and can be used to follow the processes occurring in the fiber both during immersion in the coagulation bath and during the mass transfer process.

Discussion
Let us discuss how else we can interpret the results of our model. If we arrange the visualizations of the final states of the cells of the system with different PAN contents in the system as a sequence ( Figure 1B), a resulting sketch (see Figure S4) of A-C snapshots (I) can be correlated with the change in fiber structure from the surface to the center. At the same time, the set of the same snapshots 1-3 (II), arranged in a different direction in Figure S4, can be identified with the change in a fiber fragment structure over time. If we consider the simulation cell as part of the fiber near-surface layer, the formation of a large number of high-density domains ( Figure 5, region IV) can also be interpreted as the development of a rigid shell on the surface of the nascent fibers. Such a shell can significantly slow the mass transfer between the inner regions of the fiber and the coagulation bath. As a result, a core-shell structure is formed (a loose core with capillaries and voids filled with a solvent and a dense surface). The presence of a dense surface layer with numerous nanopores can also cause the fiber cross-section to deviate from a round shape. This can occur as a result of a non-isotropic distribution of nanopores relative to the symmetry axis of the fiber. Thus, our model clearly shows that the formation of precursor fibers with uniform structure is achieved at the high content of the good solvent in the coagulation bath [62,63]. In this way, the behavior of our model agrees with the conclusions of experimental studies and can be used to follow the processes occurring in the fiber both during immersion in the coagulation bath and during the mass transfer process.

Discussion
Let us discuss how else we can interpret the results of our model. If we arrange the visualizations of the final states of the cells of the system with different PAN contents in the system as a sequence ( Figure 1B), a resulting sketch (see Figure S4) of A-C snapshots (I) can be correlated with the change in fiber structure from the surface to the center. At the same time, the set of the same snapshots 1-3 (II), arranged in a different direction in Figure S4, can be identified with the change in a fiber fragment structure over time.
We have already written that the results obtained indicate the formation of a porous structure at the high proportion of non-solvent in the simulation cell ( Figure 3B), which is usually considered a defect that weakens the mechanical properties of fibers. On the other hand, high surface porosity may be preferred for applications that require the deposition of a stable color pigment layer on the fiber surface, as well as for the production of filters and sorbents for air/water purification [2][3][4]64]. Additionally, based on the chemistry of the thermal stabilization and the carbonization process [5], it can be assumed that high porosity can be beneficial to better oxygen penetration and the removal of evolved gasses from the fiber during the final stages of CFs production.
It should also be noted that because our model does not work with the specific chemical structures but with the density distributions of their coarse-grained equivalent representations, the results of its predictions are also applicable to other chemical fibers where the Hildebrand solubility parameters of the polymer/combined solvent are close to the values chosen in this study. This can be seen as a disadvantage and an advantage of using the mesoscale simulation technique.

Mesoscale Dynamic Density Functional Theory
We used Dynamic Density Functional Theory to describe the evolution of the PAN/ DMSO/water mixture. Details of this method can be found in this list of articles [44][45][46][47]65].
Here, we will only briefly present the main points of this simulation technique.
In DDFT, a model of a molecular system is built in the form of a cubic volume V, which is called a simulation cell. Periodic boundary conditions are imposed on the cell faces, allowing the continuum nature of the system to be taken into account. Chemical structures of polymer chains, solvent, and other components of the system are mapped onto a sequence of beads of the same volume ν. To do this, we should identify key molecular structures (objects) that play the main role in the behavior of the system. They may be fragments of polymer chains, structural elements of the filler volume, several solvent molecules, etc. As a result, each selected key molecular object in the initial chemical structure is associated with a spherical particle of its own type α in the mesoscale model, α = {1, 2, 3, . . . }. The value of ν is usually evaluated as the average volume of molecular structures corresponding to the beads of the model. Thus, σ = (6ν/π) 1/3 determines the characteristic scale of the model. Thus, the original molecular structures are replaced by a simplified (coarsened) equivalent representation, which drastically reduces the number of degrees of freedom of the modeled system. The use of coarsening makes it possible to access large-length scales and to run simulations over long time intervals, wherein information about the particular chemical composition of the system is not completely lost. The interaction between the molecular subsystems is controlled by the Flory-Huggins parameters [66,67]. They can be measured experimentally or obtained by classical molecular dynamics (MD) simulations [45].
The distribution of beads of each type is described by time-dependent density fields ρ α (r,t), α = {1, 2, 3, . . . }. Their evolution obeys the Langevin diffusion equations [44,65] where F[ρ α (r,t)] is the free energy functional, D α are bead mobility parameters, and η α (r,t) is stochastic noise with a distribution obeying the fluctuation-dissipation theorem [68]. The functional derivative δF[ρ α (r,t)]/δρ α (r,t) is the partial chemical potential µ α (r,t) of the α subsystem at point r. Thus, the evolution of the density fields is determined by random "thermal" noise and gradients of chemical potentials arising from the features of the interaction between different molecular subsystems.
In the general case, the system described by Equation (1) is in a non-equilibrium state. Therefore, additional external fields, ω α (r), are introduced to act on each subsystem. They are chosen so that at each moment of time, t, the distribution of densities ρ α (r,t) corresponds to the equilibrium state. For a system consisting of n M molecules, the relationship between ρ α (r,t) and ω α (r) at a fixed time is given by the following partial functional derivative: where k B is the Boltzmann constant, T is the absolute temperature, and Z M is the intramolecular statistical integral. The form of this integral depends on the model representation of the polymer chain. We use the Z M from [65], which was obtained using the Gaussian model of polymer chains. The expression for the free energy F[ρ,ω] can be written as the following sum: where <ρ α > is the average density of the α subsystem, and κ is the compressibility parameter (κ = 10 k B T [44]). Note that external fields ω α (r) do not affect intermolecular interactions and are only present in the part of the free energy that relates to the individual molecular subsystems of the model. The force constants ε αβ in Equation (3) define the intermolecular interactions. They can be calculated from the Flory-Huggins parameter χ αβ using the known relationship [41,67]: This equation relates the chemical properties of a real system to the interaction energy of coarse particles.
It is should be noted that the problem of calculating chemical potentials is reduced to the problem of calculating external fields ω α (r) based on the standard self-consistent mean-field method. The implementation of the DDFT method and the corresponding calculation scheme are described in detail in our previous publications [46,47].

Model and Parameterization
To start simulations, we must first: (i) determine the external conditions affecting the system under study, (ii) select the composition of the system, (iii) map all participating molecular objects onto equivalent coarse-grained representations, (iv) determine a characteristic scale of the model, and (v) calculate the parameters of interactions between all molecular subsystems.
Let us describe the principles behind our model construction. When the spinning solution comes into contact with the coagulation bath (the mixture of good solvent and nonsolvent in different ratios), the rapid diffusion exchange between the contents of the bath and the spinning solution begins. During this process, the good solvent leaves the nascent fiber while the non-solvent enters it. Phase separation of the combined solvent and polymer then occurs as their miscibility deteriorates. After some time, when the polymer reaches a high density, the mass transfer slows down, and equilibrium conditions are reached between the precipitated polymer and the combined solvent phases. Therefore, the state of the nascent fiber can be considered as a mixture of polymer, solvent, and non-solvent [22]. As a result of the nature of polymers, this is not a true thermodynamic equilibrium, and the final state depends on the conditions of the process used. The specific composition of the system depends on the residence time in the coagulation bath, the chemistry of the spinning solution, the contents of the bath, and its temperature.
Thus, in our model, we assumed that the mass transfer that occurs between the spinning solution and the coagulation bath forms the PAN fiber, and the final state of each fiber fragment modeled is a quasi-equilibrium state. We chose dimethyl sulfoxide and water as good solvents and non-solvents because they are commonly used for the production of PAN fibers [5,12,16,[27][28][29]51,[69][70][71][72][73][74]. Since the real sizes of nascent fibers of 60 µm are three orders of magnitude larger than the spatial scale available in DDFT, our studies focus on the processes occurring within a relatively small volume corresponding to different parts of a real fiber. Based on data from the literature [27,29,48,75,76] and the speculations above, we also assume that the nascent fiber consists of a mixture of polymer and residual solvent composition with varying ratios of dimethyl sulfoxide and water.
Since the internal state of our model faciliates multiple interpretations (see Figure S1 ("S" denotes a reference to the Supplementary Materials file)) in our study, we consider a modeling fiber piece to be located near a nascent fiber surface. This means that the simulation cells with different DMSO:water ratios can be understood as states of the system that occur at different compositions of the coagulation bath. This allows us to simulate the effect of different coagulation bath compositions on the structure of the near-surface layer of the fiber by varying the DMSO:water ratio in the cells.
Thus, the modeled system contains three main components: PAN, DMSO, and water. As mentioned above, the DDFT approach assumes that all molecular components of the coarse-grained model consist of spherical particles of the same volume ν, and each of the subsystems is assigned its own type of coarse-grained particles. Figure 6 shows the chosen principle of mapping the atomistic structure of the model components to an equivalent coarse-grained representation. P-type beads are used to model PAN chains, D-type beads correspond to DMSO molecules, and W-type beads correspond to water. In this case, a P-type bead is mapped onto a segment of the polymer chain, and D-and W-type beads correspond to multiple DMSO and water molecules, respectively. The choice of volume ν determines the characteristic scale of the system.

Parameters and Characteristic Scale
The intermolecular interaction force constants εαβ in Equations (3) and (4) are calculated using the Hildebrand solubility parameters δα [79][80][81], which are related to the Flory-Huggins' interaction parameters [65,67,79]: In our model, we assume that the spinning solution at the inlet of the coagulation bath has the following composition: PAN-20 vol%; DMSO-80 vol% (volume fractions are used for convenience in setting the composition of the coarse-grained model) [8,16,17,71,72]. We also believe that in the late stages of coagulation, the fiber consists of 80 vol% PAN, 20 vol% DMSO + water [30,72,77,78], because the non-solvent has not completely displaced the good solvent in the freshly coagulated fiber.
It should be noted that the choice of δ P is based on the fact that although in our model, for simplicity, we consider PAN as a homopolymer, we nevertheless mean that itaconic acid (IA) is present in its composition. The IA comonomers are introduced to PAN to weaken its intermolecular interaction because PAN contains highly polar nitrile groups (-CN), which results in a stiffer chain and an increase in the viscosity of the spinning solution. In addition, the introduction of comonomers affects the cyclization process of PAN. The percentage of comonomers is usually between 2% and 15% [5,15,16]. In our case, we assume that PAN contains 5% of itaconic acid, which determines the correction of δ P relative to the value obtained for PAN homopolymer chains. Thus, using the correction (see Section S2) of the solubility parameters for P-type beads, our model can take into account the multifunctional nature of the polymer chain structure (i.e., the presence of comonomers of several types).
The DDFT model also allows us to take into account the real molecular weights of the polymer chains. A chain length, N, of the polymer model is determined by the degree of polymerization and the characteristic ratio C ∞ (the measure of the stiffness of a polymer chain) [79,84]: where M p is the molecular weight of the chain, and M mon is the molecular weight of the monomer. Using the Bicerano model [84], we obtain for PAN that C ∞ = 7.33. At the same time, the experimental value of C ∞ is in the range of 10.1-18.5 [97]. Thus, assuming C ∞ = 10, for PAN with a molecular weight of 72,000 Da, we have N ≈ 136 beads (M mon = 53.1 g/M). Since such a model length requires sufficiently large simulation cells, to save computational cost, we reduced the chain length by a factor of 4 to N = 34 (the size of an equivalent Gaussian chain is 5.8 σ) and then rescaled all interaction parameters by multiplying them by a factor of 4 [79]. With this assumption, one bead corresponds to~40 PAN monomers, giving us a characteristic length scale of~1.82 nm. At the end of this section, it should be noted that the polymer chain in the DDFT model could also be constructed with explicit consideration of the presence of several types of alternating comonomers. For example, this has been realized in our work [45][46][47], where the size of a monomer was used as the scale unit. In this study, we use large molecular weights of the polymer chains and, at the same time, implicitly take into account the multifunctional structure of the polymer chain, as it has a relatively low percentage of randomly distributed itaconic acid.

Methodology of Simulations
In the course of the simulation, the effect of the ratio of the solvents (DMSO:water) on the structure of the nascent PAN fibers is investigated. All calculations are performed at a fixed temperature of 278 K. The set of Langevin Equation (1) are solved on a regular cubic grid with a total number of nodes N total equal to 32 × 32 × 32, with a time step ∆t = 0.5τ (τ-time scale). Since the characteristic length scale is σ = 1.82 nm, the step between the nodes for convenience is chosen to be 2 nm. Thus, the selected parameters allow us to study structure formation in a cubic cell with an edge length of 64 nm.
To set up the coefficients D α in Equation (1), we used the averaged value of the diffusion coefficients for DMSO and water equal to 3.3 × 10 −11 m 2 /s [98], which gives the time scale τ of ≈ 121 ns. For all calculations, we use the homogeneous distribution of densities ρ α (r,0), α = P, D, W as the initial state. The duration of calculations, t max , is controlled by the free energy convergence condition (|F(t 2 ) − F(t 1 )| < ε where ε is a small value) (for example, see Figure 7A). The total calculation time, t max , for each set of parameters, varies from 100 ∆t to 10,000 ∆t, which is about ≈6 µs to ≈600 µs. It should be noted that the diffusion coefficients values for water and DMSO found in the literature range from 0.22 × 10 −11 m 2 /s to 15 × 10 −11 m 2 /s [27,30,76,98], which allows us to interpret the simulation time range from 1.3 µs to 9091 µs.  Using the chosen values of the solubility parameters, the following values were obtained for the Flory-Huggins interaction parameters: χ PD = χ DW = 0 and χ PW = 9.5. The choice of χ DW should be clarified. If we use the chosen δ D and δ W to calculate χ DW , we obtain the value 9.7. This contradicts the known experimental fact about the good miscibility of DMSO and water. Therefore, to eliminate this discrepancy, we chose a zero value for χ DW . An estimate of the effect of error in determining the parameters of the interactions shows that small variations in their values have little effect on the radius of the emerging domains.
The volume fraction of PAN C P varies in the range from 20 to 80 vol%, and the total volume fraction of the combined solvent DMSO + water C D + C W is varied from 80 to 20 vol%. The ratio of solvents DMSO:water is controlled by a parameter f w = C W /(C D + C W ). Its value is changed from 0 to 1. To avoid the influence of the initial density distribution in the simulation cells, we have performed a series of three calculations using independently generated initial states for each set of parameters. The results obtained are used for averaging.
To control the multiphase separation occurring in the PAN/DMSO/water mixture, we calculate the order parameter Λ α . This property is defined as the volume-averaged difference between the square of the local time-dependent density field ρ α (r,t) and the square of the average density ρ α of particles of type α, respectively: Values of the order parameter close to zero for all subsystems correspond to a system with a uniform distribution of all system components, and Λ α > 0 indicates a phase separation in the system. Figure 7A shows typical time dependences of the free energy F and the order parameters Λ α using the example of the system with C P = 80 vol%, f w = 0.5. Their shape indicates phase separation induced by non-solvent, the initial phase of which (at t < 25 τ) is accompanied by a sharp decrease in the free energy and a rapid increase in the order parameter. To identify the nature of structural changes, we visualize the distribution of specific densities ρ α (r,t) by constructing isosurfaces. An example of such visualization for the PAN-specific density ρ P > 0.8 over time in Figure 7B clearly shows that the system undergoes microphase separation.
The intense structural changes occur at characteristic times close to the time t ps of phase separation (coagulation) ( Figure 7B, snapshots 1, 2, and 3). The time of phase separation t ps is estimated by linear extrapolation of the region of rapid change in the free energy of the system along the time axis (see Figure 7A). After t > 50 τ, the values of Λ α (α = P, D, W) slowly tend towards saturation ( Figure 7A), the structure is still evolving, but the changes that occur are (essentially) small ( Figure 7B, snapshots 3 and 4). This can also be seen in Figure 8, which shows the values of a certain fraction of PAN regions with different densities of ρ P , N(ρ α )/N total , where N(ρ) is the number of nodes with a density equal to ρ, and N total is the total number of grid nodes in the simulation cell, at different simulation times. Therefore, we stop our simulation after the end of the microphase separation (t max ), when the values of F and Λ α reach saturation ( Figure 7A). The time of this event depends on the composition of the system and varies in our calculations from 100 to 10,000 integration steps of Equation (1).
At the final state, t >> t PS , a dense porous structure of water (the "blue cheese" type) with an average pore cross-section size of ≈12 nm is present in the PAN matrix ( Figure 7B). The pores are filled with a mixture of DMSO and water. This is confirmed by the combined plot of the maximum values of the densities ρ α (α = P, D, W), shown in Figure 4. It can be seen that PAN and water form an unoriented fibrillar 3D network structure with regions of high-density PAN matrix and water domains. densities of ρP, N(ρα)/Ntotal, where N(ρ) is the number of nodes with a density equal t and Ntotal is the total number of grid nodes in the simulation cell, at different simula times. Therefore, we stop our simulation after the end of the microphase separation (t when the values of F and Λα reach saturation ( Figure 7A). The time of this event depe on the composition of the system and varies in our calculations from 100 to 10,000 i gration steps of Equation (1). At the final state, t >> tPS, a dense porous structure of water (the "blue cheese" ty with an average pore cross-section size of ≈12 nm is present in the PAN matrix ( Fig  7B). The pores are filled with a mixture of DMSO and water. This is confirmed by combined plot of the maximum values of the densities ρα (α = P, D, W), shown in Fig  4. It can be seen that PAN and water form an unoriented fibrillar 3D network struc with regions of high-density PAN matrix and water domains.

Conclusions
The article presents a coarse-grained model based on the method of dynamic den functional theory. It aims to study the structure formation caused by the occurrenc phase transitions in solutions based on polyacrylonitrile. This model allows us to obse the processes of phase separation in the nascent fiber as a result of the coagulation of spinning solution in the coagulation bath. The calculations performed allow us to fol the changes in the internal structure of the system and estimate the size of the struct features. Figure 8. Distribution profiles of the fraction of domains with density ρ P in the simulation cell at different simulation times (N(ρ P ) is the number of nodes with a density equal to ρ P = 0-1.4, and N total is the total number of nodes in the grid of the simulation cell). C P = 80 vol%, f w = 0.5.

Conclusions
The article presents a coarse-grained model based on the method of dynamic density functional theory. It aims to study the structure formation caused by the occurrence of phase transitions in solutions based on polyacrylonitrile. This model allows us to observe the processes of phase separation in the nascent fiber as a result of the coagulation of the spinning solution in the coagulation bath. The calculations performed allow us to follow the changes in the internal structure of the system and estimate the size of the structural features.
In general, the model adequately describes the phase separation processes occurring in the system as a function of the solvent quality, which is controlled by the water content in the combined solvent. Several conclusions can be drawn from the results obtained. By changing the composition of the coagulation bath, it is possible to control the fiber morphology. The greatest inhomogeneity of the PAN fiber structure occurs at a high water content in the combined solvent. In this case, a large number of high-density domains are also formed. Conversely, under conditions where a good solvent is present, a homogeneous spatial structure of the polymer matrix is formed. The latter can be favorable for multi-stage fiber drawing, helping to achieve high orientation with a high degree of crystallinity.
The model developed in this work can be used to optimize the structure of PAN nascent fibers and to better understand the kinetics of the processes involved in the formation of chemical fiber structures based on other polymers. Although our model has limitations, taking into account the complexity of the phenomena occurring during the formation of the structure of the nascent fiber (from the point of view of the theoretical description), it nevertheless gives an adequate physical picture with the assumptions made. Therefore, the results obtained give reason to expect that we have obtained an additional method to analyze the formation of chemical fibers.