MODELING OF PROTEIN AND SALT REDISTRIBUTION DURING DRYING OF A SOLUTION FROM A SQUARE CELL

Background: Drying of biological fluids and saline biopolymer solutions is an actively researched topic, since the textures of the films remaining on the substrate after evaporation carry information about the state of the organism or biopolymer. Earlier, we have shown that the texture area and the amount of zigzag crystallization patterns depend on the structural state of the biopolymer. Models of flows and particle redistribution in a round sessile drop have been described by multiple authors, but in our experiments the solution completely fills a square cell to its walls, which results in different drying conditions and dynamics. Objectives: Conducting a numerical simulation of liquid evaporation and the corresponding redistribution of BSA and NaCl particles for a square 20×1×20 mm 3 cell. Materials and methods: Liquid evaporation was simulated in the OpenFOAM continuous media modeling package using the interThermalPhaseChangeFoam module, and the redistribution of the particle concentrations (BSA, Na + and Cl – ) in the solution was simulated using the method of biased random walk on a discrete Cartesian lattice. Results: According to the obtained results, for most of the drying time, the flows in the liquid are from the corners to the center of the cell and from the diagonals to the walls of the cell, which leads to the accumulation of a significant fraction of the particles near the walls. When the liquid-air interface in the central part of the cell reaches the bottom, the surface tension forces quickly withdraw the solution to the cell walls, although a small amount of liquid can remain as a drop in the center. After complete drying, the majority of BSA and NaCl particles that did not settle at the cell walls were found to be concentrated at a distance of 1-3 mm from the cell walls, in the form of a rounded band. A small amount of BSA is also present in the central part of the cell, while the remainder of the salt is uniformly distributed over the entire area of the cell due to its greater diffusivity. Conclusions: Taking into account the experimental distributions of textures, these results support the hypotheses that textures are not formed in the absence of a biopolymer, and that zigzag patterns are formed in locations with high drying rate.

Evaporation of colloidal solutions has promising applications in a number of industries, such as medicine [1][2][3], agriculture [4,5] and microelectronics [6,7]. Films obtained as a result of the evaporation of biological fluids or saline solutions of biopolymers were investigated by various authors [8][9][10][11][12][13][14]; in particular, our experiments have shown that the structural state of the biopolymer affects the area of visible textures [15] and the number of zigzag crystallization patterns on the film [16][17][18].
Mechanisms of drying of colloid-salt solution droplets, as well as mechanisms of pattern formation on the resulting films, are an actively studied subject [19][20][21][22][23][24]. However, instead of the standard method of sessile droplet evaporation, we use a flat cell filled with solution to Modeling of protein and salt redistribution during drying of a solution from a square cell 53 improve reproducibility of the contours and the areas of the obtained films. This changes the evaporation conditions and the course of drying compared to the free drop, which does not let us draw direct analogies with the existing models of particle redistribution in the drying droplet. Due to the fact that the geometry of our cell (square 20x20 mm 2 cell with 1 mm high walls) cannot be reduced to an axially symmetric 2D case, the theoretical analysis of the problem is considerably complicated; therefore, in this study we conducted a numerical simulation of fluid evaporation and of the corresponding particle motion for a square cell.

MODEL
As a first approximation, the following assumptions were made during the simulation: 1. Available concentrations of protein and salt do not affect the overall course of hydrodynamics and evaporation: for most of the drying time, the mass fraction of BSA (bovine serum albumin) and NaCl in the solution is very low (about 0.2-2%), so their effect on viscosity, thermal conductivity, surface tension, and other physical and chemical properties of the solvent can be considered negligible; accordingly, liquid evaporation can be simulated separately from particle redistribution.
2. Inhibition of the Marangoni effect: although in theory the contribution of the Marangoni flows must be noticeable in evaporating droplets of pure water, in practice their effect is very weak due to residual contamination by surface-active substances [25].
3. The effect of contact line pinning is ignored (due to the complexity of modeling this phenomenon in the OpenFOAM package). 4. Absence of particle aggregation: the present concentration of NaCl salt is not large enough to completely compensate the charge of BSA, so during the simulation it was assumed that protein molecules effectively do not aggregate.

Simulation of fluid evaporation in a square cell
The simulation of fluid evaporation was performed in the OpenFOAM continuous media finite element modeling package. For hydrodynamic problems with two phases this package uses a volume-of-fluid method, where the amount of liquid in each spatial element of the computational domain is represented by the "volume-of-fluid" quantity  . To model the phenomena of mass and heat transport, surface tension and phase transition, the interThermalPhaseChangeFoam module [26,27] was used, which solves the conservation equations for mass (using the Poisson equation), momentum (taking into account viscous forces, pressure, surface tension and gravity), thermal energy (taking into account the convection and heat conduction phenomena) and volume-of-fluid (with corrections to compress the interface and counteract numerical diffusion). The cells we use for obtaining films have a volume of 20x1x20 mm 3 . Since the initial volume of the solution is 0.5 ml, the height of the liquid can reach 1.3 mm, so an additional 1 mm of air above the cell was included in the computational domain. To accelerate the simulation, the domain was limited to a rectangular volume of 10x2x10 mm 3 along the axes X, Y and Z, respectively ( Fig. 1), which corresponds to a quarter of the cell. This volume was divided into cubic 0.1x0.1x0.1 mm 3 finite elements in the form of a three-dimensional grid (100x20x100 elements along the corresponding axes).
To assign the boundary conditions, the surface of the computational domain was divided into 4 areas: the XY side, the YZ side, the surface below 1 mm ("cell") and the surface above 1 mm ("air"). The boundary conditions for these areas are given in the Table 1.  Table 1. Boundary conditions for the values used in modeling of evaporation. N   denotes the gradient in the direction normal to the boundary surface. The value of the wetting angle was taken from [28], the temperature corresponds to the experimental values [16,17], and the rest of the boundary conditions were chosen in accordance with the conventional practices [29] "Cell" area "Air" area XY and YZ sides  In our experiments, it was observed that after exp t = 20 min in the drying chamber, 0.060±0.003 g of a solution evaporates from the cell, and 0.224±0.003 g after exp t = 80 min. Correspondingly, 0.5 ml of a solution completely evaporates in 3 hours on average.
Due to the simulation taking a lot of time to compute, it was decided to start with 0.333 ml of water (which corresponds to exp t = 1 hour). Assuming that the nature of the currents within the liquid stays the same for most of the drying time, this shouldn't significantly affect the final deposition distributions. To account for the liquid meniscus, which forms near the cell walls due to surface tension, the field  was initialized using the formula (1) where E L = 0.1 mm is the size of the finite element, C h = 1 mm is the height of the cell walls, 2 / C w = 10 mm is the half-width of the cell; the rest of the parameters were chosen empirically to approximate the surface of 0.333 ml of liquid: W h = 0.75 mm (height of the liquid at the cell center), M w = 3.5 mm (width of the meniscus), P = 1/3 (empirical coefficient determining the shape of the meniscus). The min() and max() denote the minimum and maximum of a set of values, respectively (the same notation is used throughout the paper). Modeling of protein and salt redistribution during drying of a solution from a square cell 55 The simulation was performed using the incompressible Newtonian fluid and laminar flow approximations, the Brackbill surface tension model [30] and the "Interface Equilibrium -No Dilatation" phase change model [26]. The properties of the liquid and gaseous phases that were used in the simulation are given in the Table 2. Thermal conductivity of water, W/(m·K) 0.631 [31] Specific heat of water, J/(kg·K) 4178 [31] Kinematic viscosity of water vapor, m 2 /s 20.5·10 -6 [31] Density of water vapor, kg/m 3 0.598 [31] Thermal conductivity of water vapor, W/(m·K) 0.025 [31] Specific heat of water vapor, J/(kg·K) 2077 [31] Surface tension, kg/s 2 0.0695 [32] Phase change enthalpy, J/kg 2264705 [33] Since the direction of the phase change (i.e., whether evaporation or condensation takes place) in interThermalPhaseChangeFoam [26,27] is determined by the difference between the temperature of the finite element T and the saturation temperature sat T of the liquid, we had to set the saturation temperature lower than the temperature of the cell ( sat T = 10 °C to reduce the computational time) and to use a model which did not incorporate formation of vapor bubbles. As a result, the rate of evaporation in the simulation was 2-3 orders of magnitude higher than in the experiment; however, the qualitative nature of the obtained currents can still be useful as a first approximation in particle redistribution modeling. In order to correlate the simulation results with the experiment, simulation time sim t was matched to the "experimental time" exp t based on the volume of water remaining in the cell (2): where the total drying time total t = 10800 s, the volume of the finite element E V = 10 -12 m 3 , the initial volume of the solution 0 is the three-dimensional index of the corresponding finite element. In the simulation of particle redistribution, velocity values U  were also scaled accordingly.

Modeling of redistribution of colloidal particles and dissolved salt
To simulate the particle redistribution in the solution, biased random walk on a discrete Cartesian lattice (with cubical elements) was used, where the probability of particle movement between adjacent lattice cells depends on the vector of fluid velocity (the details are laid out in the "Calculation of the probability of particle motion" section). The same number of elements (100x20x100) was used as the size of the discrete lattice. BSA and salt concentrations were initialized as evenly distributed over the volume of the liquid; assuming that the nature of the currents within the liquid stays the same for most of the drying time, this shouldn't significantly affect the final deposition distributions. The algorithm for simulating particle motion was based on the ideas described in [34,35]. The overall structure of the algorithm is depicted in Fig. 2. During the simulation, BSA molecules and dissolved salt ions (Na + , Cl -) were considered separate types of particles. Initial concentrations and physical parameters of these particles are shown in Table. 3   Table 3. Parameters of particles used in calculations. For BSA, p r was taken from [36], p V from [37], p m from [38]. For Na + and Cl -, p r was taken from [39], p V was calculated as the volume of spheres of Van der Waals radius [40], and p m was taken from [41] Parameter BSA Na + Cl - Interpolation of the U, h, α fields at time t The fields of velocities t U  and volume-of-fluid t  at a time t were obtained by linear interpolation of snapshots (3)(4)(5), which were recorded at regular intervals during the simulation of fluid evaporation. Since water-air interface can move more than one lattice cell per snapshot, interpolation of the field  was not performed directly. Instead, a height-map n h was constructed for each n  snapshot, and the field t  at time t was calculated from the interpolated height-map t h (5): Interpolation of the U, h, α fields at time t

Calculation of the simulation step Δt
Calculation of volume fraction and diffusion coefficient of the particles Calculation of the probability of particle motion

Redistribution of particle concentrations
Next simulation step t = t + Δt If t < t end Initialization of particle concentration distribution Plotting particle concentration distributions where n U  , n h are snapshots of the velocity field and the height-map at the time n t , the , the size of the finite element E L = 0.1 mm.

Calculation of the simulation step Δt
The simulation of motion of colloidal particles was performed using the exp t time scale (i.e. it started at exp t = 1 hour and ended at exp t = 3 hours) and a variable time step.
To ensure the correctness of the calculations, the duration of the simulation step t  must be such that the distance traveled by any particle during the step is not larger than the linear size of the lattice element E L : where max U is the maximum projection of velocity in this simulation step It was decided to limit the maximal particle distance per step to 1/3 of lattice element size, as a compromise between the simulation fidelity and the time spent on computation. Correspondingly, Calculation of volume fraction and diffusion coefficient of the particles Each element of the lattice may contain colloidal particles (BSA) and ions of dissolved salt (Na + , Cl -). In contrast to the papers [34,35], the number of particles in this case is too large for them to be modeled individually. Thus particle movement was simulated as the redistribution of the amount of substance between the adjacent lattice elements.
If we consider colloidal (BSA) and dissolved (Na + , Cl -) particles to be solid spheres, their amount in the lattice element is limited by the condition where   is the amount of substance of particles of type p in the corresponding lattice element, A N is the Avogadro number, p V is the volume of a particle of type p , E V is the volume of the lattice element.
The value of max  is determined by the packing geometry, which depends on the shape and size distribution of the particles, but not on the particle size. For solid spheres, the theoretical value is usually taken max  = 0.64 [42], which corresponds to random close packing. However, according to the experimental studies of solid sphere dispersions [43,44], molecular mobility disappears when the volume fraction reaches the value of colloidal glass transition g  [43,44]. For monodisperse solid spheres g  = 0.58 [45,46].
Our system is polydisperse (BSA, Na + and Clparticles have different size), so max  of each lattice element will depend on its content. Since max  estimation formula appears to only exist for the bi-disperse case in the published literature [47], here we use it under the assumption that Na + and Clcan be reasonably treated as one species (compared to BSA): where mono max  is the maximum packing fraction of a monodisperse system ( g  in this case), large  is the volume fraction of large particles, small  is the volume fraction of small particles.
Freedom of particle movement in a lattice element ) , , ( k j i can be characterized as This value modulates the number of particles that can move between the lattice elements in a single simulation step, and is also used in the Krieger and Dougherty's semi-empirical formula [48] for calculating local viscosity of the dispersion: where F  is the dynamic viscosity of the liquid. In addition to fluid advection, the particles can also move between lattice elements due to diffusion. The diffusion coefficient p D for a spherical particle p can be estimated from the Stokes-Einstein equation: where B k is the Boltzmann constant, T is the temperature, p r is the Stokes particle radius. Calculation of the probability of particle motion When simulating random walk on a grid, during one simulation step a single particle may, with a certain probability P , move into one of the neighboring lattice elements or stay in place (mutually exclusive events). In the case of a large number of particles, these probabilities can be considered proportional to the amount of substance n  that has moved in the corresponding direction:             ,  ,  ,  ,  ,  ,  ,   ,  ,  ,  ,  , , where   k j i n , , is the amount of substance in a lattice element with a three-dimensional index ) , , ( k j i , and } 1 , denote the direction of particle movement. For biased random walk, particle movement along the direction of the flow should be more likely than the movement in the opposite direction. Three methods for calculating the corresponding probabilities were considered.

Particle motion probability calculation: method 1
The first method was described in [34], where the probabilities of particle movement depended only on the velocity field: where x P  , x P  , y P  , y P  , z P  , z P  are the probabilities of movement in the direction of the corresponding axes, w P is the probability of the particle staying in place, N = 3 is the dimensionality of space, x U , y U , z U are the components of the velocity field vector at a given point, max U is the maximum projection of velocity in this simulation step. To take into account the duration of the simulation step, the probabilities of transitions were calculated in the following way: where t  is the simulation step duration, E L is the size of the lattice element.

Particle motion probability calculation: method 2
The idea of the second method was taken from the paper [35], according to which the probability of particle movement between lattice elements is proportional to the Boltzmann distribution: where S P is the probability of finding a particle in the state S , S E is the energy of state S , B k is the Boltzmann constant, T is the temperature, exp() is the exponential function. In this case, the S E corresponds to the work that a particle must perform against the Stokes' drag to move to a certain distance in a given direction: where U  is the vector of liquid velocity at a given point, v  is the direction vector along which the particle may move during the simulation step t    where ijk  is the unit vector in the direction of the neighboring lattice element (or a zero vector in the case of the particle remaining in place).
Combining the equations (20)(21)(22), we obtain where x P  , x P  , y P  , y P  , z P  , z P  are the probabilities of movement in the direction of the corresponding axes, w P is the probability of the particle staying in place, x U , y U , z U are the components of the velocity field vector at a given point, p D is the diffusion coefficient at a given point, C is an arbitrary constant. Accordingly, Particle motion probability calculation: method 3 In the third method of probability calculation, a simplified geometric interpretation of particle motion was used (Fig. 3). Assuming that at a time t the particles uniformly fill the volume which will partially intersect with several lattice elements. Accordingly, the probability of particles to remain in the element ) , , ( k j i or to move to a neighboring lattice element will be proportional to the volume of the intersection of the cube ' V with the volume of the corresponding cubic element:

Additional rules for calculating probabilities and redistribution of concentrations
The following additional rules were also used in the modeling of particle redistribution.
In the case when a lattice element contains no water ( ) , , ( k j i  = 0), the particles of this element are only allowed to move vertically down (or remain in place, if j = 0). It simulates the effect of gravity and prevents the movement of particles in completely dried regions.
On the boundaries of the computational domain, particle motion is limited by the corresponding boundary conditions (motion beyond the cell walls and the cell bottom is prohibited, and is reflected at the symmetry planes). Other types of particle interaction with the cell walls and bottom (e.g. sticking to the surface) were not modeled in this study.
Since particles can not move beyond the liquid-air boundary and should be less likely to move to densely packed lattice elements, the motion in the direction of the neighboring elements is also limited as follows:  , where  is the amount of fluid in the element,  is the degree of mobility (14).
After applying these rules, the probability sum can be less than one, so the probability of particles remaining in the element ) , , ( k j i is recalculated as 0  ,  0  ,  0  ,  ,  ,  ,  ,  ,  ,  ,  '  1  0  ,  0 In the absence of an external electric field, the concentrations of Na + and Clions, averaged over the volume of the cubic lattice element, can always be considered the same, so during the simulation, the redistribution of Na + and Clparticles between the lattice elements was synchronized.
To prevent accumulation of rounding errors, at the end of each simulation step the amount of substance of BSA, Na + and Clin the lattice elements was normalized to maintain the total number of the corresponding particles.

RESULTS AND DISCUSSION
The course of liquid evaporation, obtained as a result of numerical simulation, is shown in Fig. 4-7. Until the surface of the liquid reached the bottom of the cell, the convective flows in the solution exhibit the character illustrated in Fig. 4: a diagonal flow from the corner to the center of the cell, as well as the circulation between the diagonal and the cell walls (spiraling toward the cell corner). As the evaporation progresses, the flow becomes more chaotic (Fig. 5a), but its general character is maintained almost up to 2/3 of the total drying time (Fig. 5b). When the liquid-air interface at a certain point reaches the cell bottom (Fig. 6a), the surface tension forces quickly withdraw the solution to the cell walls, although a small amount of liquid can remain as a drop in the center (Fig. 6b). After this, the flows stabilize, resulting in a circulation between the cell bottom, the cell walls and the surface of the liquid (Fig. 7a), and the rest of the time the liquid dries toward the cell walls (Fig. 7b). The distribution of protein and salt concentrations for the first method of probability calculation is shown in Fig. 8 (a, b), respectively. Since the probability of particle movement in this method depends only on the fluid velocity field, the resulting BSA and NaCl distributions are the same. At distances < 1 mm from the cell walls, there is an increased concentration of particles; at a distance of 1-3 mm from the walls (where the solution recedes after the "hole" is formed in the liquid surface), the distribution is slightly smaller and uneven; at a distance of 3-7 mm (the central part that quickly lost all the liquid) the distribution is relatively homogenous; finally, there is a diagonal strip of slightly lower concentration, which likely reflects the nature of the diagonal currents. The second (Fig. 9) and the third (Fig. 10) probability calculation methods take into account particle diffusivity, therefore they give slightly different and substantially less uniform distributions for protein and salt. Except for a few points at the cell walls with excessive concentration of particles, BSA and NaCl are mainly concentrated in a rounded band that is located at a distance of 1-3 mm from the cell walls. Unlike protein (Fig. 9a, 10a), a notable amount of salt is also uniformly present outside of this band (Fig. 9b, 10b), due to greater diffusivity of NaCl particles.
(a) (b) Fig. 9. Distribution of BSA (a) and NaCl (b) for the 2nd method of probability calculation.
Cell center is in the bottom left corner of the image.
(a) (b) Fig. 10. Distribution of BSA (a) and NaCl (b) for the 3rd method of probability calculation.
Cell center is in the bottom left corner of the image.

1-3 mm band
These results indicate that the weakness of method 1 (not accounting for particle diffusivity) is much more severe than the weakness of method 2 (exponential dependence of particle motion probability on the fluid velocity). Since method 3 lacks these disadvantages, we expect the results obtained with it to be the most accurate among the three.
Finally, comparing the results of the simulation with the experimental distributions of textures on films (Fig. 11a, 11b), and zigzag patterns (Z-patterns) in particular (Fig. 12b), the following observations can be made: 1. Textures stop forming at a distance of 1-1.5 mm from the cell edges, which approximately corresponds to the boundary of the 1-3 mm band in Fig. 9 and Fig. 10. Since pure salt solutions (without biopolymer) do not create textures [16], the simulation results suggest that very low BSA concentration between the 1-3 mm band and the cell walls is a likely reason for why textures do not cover the whole area of the cell.
2. Z-patterns (Fig. 12a) are concentrated approximately at the same distance from the cell edges as the distance to which the simulated liquid retreats after the "hole in water" is formed. Z-patterns are also practically absent in a small central part of the cell, which approximately corresponds to the droplet that separated from the main body of liquid (Fig. 12c). This suggests that rapid drying is one of the main prerequisites for the formation of Z-patterns.
Cell center is in the bottom left corner of the image.

CONCLUSIONS
Thus, the distributions of BSA and NaCl, obtained via a simulation of solution drying in a square cell, approximately correspond to the experimental texture distributions on the cell surface. The results support the hypotheses that the area of the textures depends on the local biopolymer concentration, and that zigzag patterns are formed in locations with high drying rate. In the future, we intend to combine fluid dynamics and particle motion and account for additional phenomena to improve the accuracy of the model.