Experimental & Computational Fluid Dynamics Study of the Suitability of Different Solid Feed Pellets for Aquaculture Systems

: Fish feed delivery is one of the challenges which ﬁsh farmers encounter daily. The main aim of the feeding process is to ensure that every ﬁsh is provided with sufﬁcient feed to maintain desired growth rates. The properties of ﬁsh feed pellet, such as water stability, degree of swelling or ﬂoating time, are critical traits impacting feed delivery. Some considerable effort is currently being made with regard to the replacement of ﬁsh meal and ﬁsh oil with other sustainable alternative raw materials (i.e., plant or insect-based) with different properties. The main aim of this study is to investigate the motion and residence time distribution (RTD) of two types of solid feed pellets with different properties in a cylindrical ﬁsh tank. After experimental identiﬁcation of material and geometrical properties of both types of pellets, a detailed 3D computational ﬂuid dynamics (CFD) study for each type of pellets is performed. The mean residence time of pellets injected at the surface of the ﬁsh tank can differ by up to 75% depending on the position of the injection. The smallest residence time is when the position is located at the center of the liquid surface (17 s); the largest is near the edge of the tank (75 s). The maximum difference between the two studied types of pellets is 25% and it increases with positions closer to the center of the tank. The maximum difference for positions along the perimeter at 3/4 tank radius is 8%; the largest residence times are observed at the opposite side of the water inlet. Based on this study, we argue that the suitability of different solid feed pellets for aquaculture systems with speciﬁc ﬁsh can be determined, and eventually the pellet composition (formula) as well as the injection position can be optimized.


Introduction
Modern aquaculture starts to be more reliable upon using recirculating aquaculture systems (RAS) as an effective way to produce large quantities of fish using small footprints for building [1]. The rule of thumb in such intensive systems is to remove solid particles or uneaten food as soon as possible.
In some system designs, in-tank solid particle separation by devices as a dual drain and particle traps,

Experimental
Two different pellets, from the point of view of composition, were used in the study. The first one, TM0, contained fish meal, plant based ingredients and fish oil as protein and lipid source, respectively. In the second one, TM75, 75% of fish meal was replacement by insect meal (see Figure 4, right, for illustration of these pellets).
Feed formulation and proximate composition is presented in Table 1. Feed mixer HLJ-700/C (Saibainuo, China), dual-screw extruded SLG II 70 (Saibainuo, China) and seven-layers air dryer KX-7-8D (Saibainuo, China) were used for the preparation of experimental feeds. Temperature and pressure ranged from 96 to 106 • C and from 19 to 22 atm during the course of feed production, respectively. Maximal temperature of 138 • C was used during drying process, which lasted 25-30 min. Oil was applied after extrusion. The final dry matter content was 7.8-8%.

Pellet's Diameter Measurement Using Image Processing
The pellets' diameters were measured using image processing to acquire the appropriate size distribution of the pellets accurately. Pellets from each group were scattered on a white paper and placed in a single-layer white fabric tent to ensure uniform light intensity over the sample. They were then photographed with a 24-megapixel Nikon D3300 digital camera (Nikon Corp., Tokyo, Japan) under a fluorescent lighting system. All settings of the camera were on manual, with the following parameters: exposure mode manual, shutter speed 1/40 s, aperture f/4.0, and ISO sensitivity 350.
Within the tent, the digital camera was positioned vertically 46 cm from the sample. The images were recorded in Nikon raw format (NEF) and transferred to a computer for further processing. Each pellet was segmented from the background based on histogram thresholding. The threshold, T = 60, was found to be the most suitable one for the segmentation of pellets. The binary image [BW(x, y)] of the pellets could be defined as: where f is the original image, x, and y are coordinates of threshold value points, and T is the threshold value. Subsequently, any image pixel ≥ 40 was labeled as 1 (white), while all other pixels were classified as a background and were labeled 0 (black). Afterward, pellet diameters were measured ( Figure 2). The pixel size of the projected scene was 0.0461mm. It was calculated from the camera resolution and the distance between the camera and the sample. The accuracy of the diameter measurement is mainly influenced by the pellet edge localization (±0.5 pixel) and the image distortion close to the border of the image. The accuracy was estimated as ±0.0461 mm. All images were processed using MATLAB image processing toolbox (MathWorks Inc., Natick, MA, USA).

Pellet's Size Distribution
Commonly used representation of droplets or particles size distribution is the Rosin-Rammler distribution, which describes the mass fraction of particles with a diameter greater than d, according to the following equation Parameters d and n in this distribution were fitted on our experimental data obtained in the image processing experiment (see the previous section) as for TM0 pellets (see the previous subsection for more details about pellets TM0 and TM75). Pellets with diameters below 0.5 mm were excluded from the data set, giving 173 data points in the fitting procedure. Figure 3 shows the graphical representation of the Rosin-Rammler distribution. For the different pellet type, TM75, the following parameters of Equation (2) were obtained (number of data points was 349).   (2)).

Effective Density of Pellets
Because the pellets are porous, it is crucial to determine their effective density in water, which could then be used in a CFD analysis. It can be determined in an experiment based on measuring the settling velocity. The settling velocity of particles in a fluid is derived from the balance of gravitational, drag and lift (buoyancy) forces. For spherical particles, the drag coefficient is determined by the following equation [15] where d p is particle diameter, u the characteristic velocity, p and f are effective density of the particle and liquid phase, respectively. The drag coefficient C d is usually expressed in terms of Reynolds number, Re = ud p /ν, where ν = µ/ is kinematic viscosity. An example of a correlation describing this dependency could be where constants a 1−3 are chosen according to corresponding range of Reynolds number [12]. Haider and Levenspiel [13] presented the following equation, which takes into account the non-spherical shape of particles The parameters b 1−4 in (7) are defined as follows where shape factor ϕ (sphericity) represents the ratio of the surface of a sphere with the same volume as the real particle (s) and surface S of the real particle, i.e., ϕ = s/S. The shape factor cannot exceed a value of 1, and it is (2/3) 1/3 = 0.874 for a cylinder with height equal to diameter. Equation (7) is used in ANSYS Fluent solver [16], which was used to perform CFD simulations in this work. This correlation is reported to be valid for Re < 2.5 × 10 4 .
Settling velocities of several particles in a water column (see Figure 4, left) were measured [17]. An HD camera with 29.879 fps recorded the particle trajectory within a water column of total height 300 mm. The measured settling trajectory was then within the bottom part of it (224 mm), where constant (terminal) velocity could be expected. The settling trajectory measurement accuracy was estimated as ±0.14 mm (based on the camera resolution). The time necessary to evaluate the settling velocity was determined from the recorded movie with accuracy limited to one frame. The frame rate accuracy of the camera was estimated as 0.06%; therefore, the maximum time measurement inaccuracy is 0.004 s (for max. time 7.1 s, see Table 2), and the inaccuracy for the determined settling velocity should be below 1%. Assuming that the characteristic size (diameter) of the particle is known, the effective density of the particle could be expressed from Equation (5) as where u is the determined settling velocity. The characteristic size d p was based on a static camera picture with estimated accuracy 0.06 mm. (right) Photo of feed pellets TM75 (75% insect + 25% fish-based source of lipid and protein). Table 2. Experimental results of settling velocity for TM0 pellets [17]. The resulting effective density of particles expressed from Equation (9) Table 2 summarizes the measured settling velocities of 10 selected TM0 feed pellets. Using these values, corresponding Reynolds numbers and drag coefficients according to Equation (7) have been calculated. Shape factor ϕ = 0.874 (cylinder with equal height and diameter) was used in the correlation. Then, average value of the effective density can be evaluated from Equation (9), in this case it was 1054.4 kg/m 3 with confidence interval ±1.1%. The fluid (water) used in our calculations corresponded to temperature 20 • C (density f =998.2 kg/m 3 and kinematic viscosity ν = 1.004 × 10 −6 m 2 /s).

Settling Velocity and Effective Density Results
The same experiments were done using TM75 feed pellets; see Table 3 for resulting settling velocities. According to Equation (9), the resulting effective density of these pellets was determined as 1072.9 kg/m 3 ± 1.9%. Table 3. Experimental results of settling velocity for TM75 pellets [17]. The resulting effective density of particles expressed from Equation (9)

CFD Simulations
We used ANSYS Fluent software to perform CFD simulations of fluid flow in a fish tank, including particles representing the feed pellets. Due to the better use of the RAS rearing volume, see [10,18] and references within there, the cylindrical geometry of the fish tank (see Figure 1) was used in this study. Particles properties were determined in the previous experimental part of this work. We used the classic DPM (discrete phase model) approach, where the particles and their trajectories are modeled in the Lagrangian reference frame, whereas the continuous phase (water, in our case) is modeled in the Eulerian reference frame. We assumed that the particles representing the feed pellets do not have substantial impact on the velocity field, that is, we assumed one-way coupling between the continuous phase (water) and the discrete phase (pellets).

CFD Setup
The geometry used in CFD simulations corresponded to 1.4 m 3 cylindrical fish tank with inlet at the top close to the wall, and outlet at the bottom, operated at the Research Institute of Fish Culture and Hydrobiology, University of South Bohemia (see Figure 1). The same fish tank served for the validation of CFD analysis of the velocity flow field in the MSc thesis [19]. The mesh created in ANSYS Fluent Meshing had 420,000 mesh elements in total, using hexahedral mesh elements in the core of the geometry and polyhedral mesh elements between boundary layer mesh elements and the core of the geometry, see Figure 5 for illustration of the mesh in 3D fluid domain. This is a relatively new feature in ANSYS Fluent software package, which should deliver better efficiency and accuracy compared to purely tetrahedral or polyhedral meshes (see Mosaic technology for more details [20]). Grid convergency index based on average velocities in three horizontal planes was evaluated as 11.6% for this size of the mesh.
The inlet flow rate was 0.5 L/s, representing 1.27 of the tank volume per hour and residence time 0.787 h. Corresponding mass flow rate boundary condition was set at the inlet pipe. Properties of the water inside the fish tank used in our simulations corresponded to 20 • C (density 998.2 kg/m 3 and kinematic viscosity 1.004 × 10 −6 m 2 /s). Turbulence modeling using ANSYS Fluent is rather technical; on the basis of our previous experimental measurements and subsequent studies [19,21], the SST k − ω intermittency turbulence model was used for the description of the turbulent flow. This model belongs to the family of transition RANS models, combining k − ω model in the near-wall region and k − ε further from the wall, along with another equation for intermittency variable, covering the laminar-turbulent transition process [22,23]. This model should give better predictions of flow field in cases where the fully developed turbulent flow is not established in the whole fluid domain. After the developed steady-state velocity field in the tank was obtained (it took approx. 70 min on 3 cores of Intel i7-8700 processor), particle tracking followed. As one-way coupling between the solid and fluid phase was assumed, the particle tracking was part of post-processing procedure based on the steady-state solution. It took up to 20 min for the largest amount of feed pellets (42 thousand) injected at the liquid surface. Corresponding characteristic size (diameter), density and sphericity for the discrete phase (feed pellets) based on the experiments described in the previous section were used. One kilogram of particles were injected at the surface of the liquid (water) in the tank at the beginning of the simulation.
The validation of the flow field obtained in CFD simulations with experimental data was done according to [19], where flow velocities in different positions of the tank were measured using acoustic doppler velocimeter (FlowTracker). The differences varied between 0.2% up to 37%, depending on the position where the velocity was measured. The average difference between the measured data and simulation data in 24 positions of the tank is 11%. A comparison for 6 positions located on 2 vertical lines (at 3 different distances from the bottom) is illustrated on Figure 6. sim. data exp. data

CFD Results
CFD simulations were performed to obtain a description of particle distribution in the tank. Flow velocities in the tank are relatively small (maximum is around 0.8 m/s, but the average is 0.08 m/s), see Figure 7 illustrating contours of velocity magnitude in three horizontal planes (on the left). On the right, the vertical velocity components are depicted, and comparing with Table 2, we can see that their values are smaller than the measured settling velocities; therefore, they cannot simply prevent the pellets from settling down to the bottom relatively quickly. In some locations (near the wall and at the bottom central part near the outlet), the vertical velocity components are even negative, so they actually accelerate the pellets on their way down to the bottom there. So, if we wanted to increase the residence time of pellets in the tank for some reason, we should probably use a different configuration, for example a different position of the inlet or different inlet flow rate. Or, different properties of pellets (size and effective density) could affect the residence time as well, and they could be a part of the parameters (variables) to be optimized.   The mean residence time for pellets according to the identified Rosin-Rammler distribution parameters (see previous section) is 63.95 s ± 0.21% for TM0, and 65.35 s ± 0.25% for TM75, in the case of inlet flow rate 0.5 L/s. Figure 9 illustrates the normalized residence time distribution for TM0 pellets described by the identified parameters of the Rosin-Rammler distribution, Equation (3), with minimum and maximum diameters 1.9 and 4.74 mm, respectively. Table 4 presents the summary of mean residence times for three different inlet flow rates.   Table 5 and Figure 10 describe the impact of the injection position on the mean residence time. Samples of 500 particles reflecting the Rosin-Rammler distribution parameters determined in the experimental part were injected at the liquid surface in a circular area with radius 7.5 cm around the corresponding points 0-8 depicted on the left of Figure 11. Points 1-8 were placed at the circle with 3/4 of tank radius R. This dependency could help us to find out an optimum position with, for example, largest residence time; see Figure 10 describing the mean residence time dependency with respect to the angular coordinate. TM0 pellets have maximum mean residence time at point 5 (angle 180 • ), and TM75 pellets at point 6 (angle 225 • ), that is, they are situated on the opposite side of the inlet. Table 5. Impact of the injection position on the mean residence time. Positions of injections 0-8 are depicted on Figure 11. In all cases, samples of 500 particles reflecting the Rosin-Rammler distribution parameters determined in the experimental part were injected at the liquid surface in a circular area around the corresponding points. TM0 TM75 Figure 10. Dependency of the mean residence time with respect to the position angle (see Figure 11 on the left) along the perimeter of circle with 3/4 of the vessel radius. The inlet flow rate used to obtain these data in CFD simulations was 0.5 L/s.  Table 5).
Notably, 500 particles were injected in a circular area of radius 75 mm around these points (group injection type + staggered positions in ANSYS Fluent). (right) Illustration of TM0 pellets trajectories injected at position 5. Different colors represent particle diameters (based on the Rosin-Rammler distribution parameters).
When searching for an optimum (maximum) mean residence time, we might be interested in the radial coordinate as well. Figure 12 shows such dependency of the mean residence time at angular position, corresponding to point 5 in Figure 11. The mean residence time is monotonically decreasing towards the center of the fish tank. In addition, Figure 13 illustrates CFD results describing the impact of the pellet's density on the mean residence time of the pellets injected at position 5 (see Figure 11). It is clear that, for densities closer to the density of water (998.2 kg/m 3 in our simulations), the mean residence time increases substantially. Therefore, if there is a need for a particular residence time, such dependency could provide necessary information. For different geometrical and operating parameters, we can expect similar results, even though quantitatively different.  Figure 11), that is on the opposite side of the water inlet. The inlet flow rate used to obtain these data in CFD simulations was 0.5 L/s.

Discussion
The CFD simulations done in this work illustrate the impact of various parameters on the mean residence time of feed pellets in a cylindrical fish tank (Figure 1), which might be quite significant concerning the sustainability of such aquaculture systems. For example, the impact of the pellet characteristic size (diameter) is illustrated on Figure 8, and it is clear that with a smaller size, the mean residence time increases. A similar trend can be observed for lower flow rates, see Table 4.
Another important parameter affecting the mean residence time of the feed pellets in a fish tank is their injection position. It is usually located at the liquid surface, and the exact position could be a part of an optimization procedure in CFD simulations. In this particular case of the cylindrical fish tank with an outlet at the bottom center, the smallest residence time is at the center (it is 17.0 s for TM0 pellets, and 21.3 s for TM75), and it increases when moving towards the outer edge of the tank in the radial direction (75.5 s for TM0, 75.3 s for TM75). The largest difference between the two types of pellets is at the center, up to 25%. However, there might exist an optimum (maximum) residence time with respect to the angular coordinate, see Figure 11. The maximum difference along the perimeter at 3/4 tank radius is 8%. Of course, the results could be different for different geometrical and operating parameters. For example, if we switched the inlet and outlet, that is, if the water inlet was at the bottom of the tank, we could substantially change the flow pattern and affect the residence time of pellets as well as their trajectories. Once the methodology described in this paper is well established, this could be a relatively easy task.
Concerning the feed pellet properties which were determined in the experimental part of this work, they can change substantially with the time spent in water (see [11] describing the expansion rate of feed pellets, for example). Implementing such dependency might not be an easy task in ANSYS Fluent solver; a user-defined function (UDF) would have to be created based on some model describing experimental data. Moreover, other parameters of the pellets could change; for example, their resistance to higher shear stresses will probably decrease with larger residence time in the water so that they could disintegrate into smaller particles of quite different properties than at the beginning. Because the amount of pellets used for fish feeding in RAS systems is relatively very small, they were modeled as completely passive in our CFD simulations, that is, we assumed the so called one-way coupling between the fluid and solid (pellets) phases. The feeding is usually done several times per day with a small amount of pellets. Therefore, the effect of pellets settling to the flow field is very low. With a larger amount of pellets, two-way coupling could be used, and consequently, the impact of particles on the flow field, turbulence quantities and other variables could be investigated.
This article can be served as a basis for the next simulating of feed particle transport and modeling of zones of sedimentation and zones of uniform distribution. This can be utilized in the feed production industry, as this sector is constantly testing novel feed ingredients, which influence feed pellets' properties. The CFD-based methodology developed in this study can also be used in the optimization of tank design and recirculating aquaculture technology. Future studies should aim to effect inlet modification, outlet modifications, and aeration equipment to optimize feed particle distribution within rearing tanks.

Conclusions
In this work, the motion (flow pattern) of solid feed pellets in a fish tank was simulated (using ANSYS Fluent software), in order to discern pellet suitability for aquaculture systems. First, we designed the experiments to determine the effective density of feed pellets in water (based on the settling of pellets in a water column and determining the settling velocity). Having determined the properties of the pellets, 3D CFD simulations of fluid flow and the distribution of feed pellets in a fish tank were performed. The mean residence times which might be one of the important parameters in the design of fish tanks were evaluated for various operating and geometrical parameters, namely pellet characteristic size (diameter), volumetric flow rate, density, and the location of the pellet injection. Knowing that the quantitative results depend on both (i) the specific design and operating conditions of the fish tank and (ii) feed pellets properties, we present some of them. The mean residence time of pellets injected at the surface of the fish tank can differ by up to 75% depending on the injection position. The shortest residence time is when the position is located at the center of the liquid surface (17 s); the largest is near the edge of the tank (75 s): The maximum difference between the studied two types of pellets is 25%, and it increases with position closer to the center of the tank. The maximum difference for positions along the perimeter at 3/4 tank radius is 8%; the most extended residence times are observed at the opposite side of the water inlet.
As far as we know, this is the first study dealing with both the identification of properties of the pellets and correct settings of the model in CFD solver as well. Such a CFD model might be beneficial for researchers and aquaculturists who are concerned in fish farming and looking for optimal properties of the feed pellets, as well as optimal design and operating parameters of fish tanks.
Concerning our future goals in this research field, we want to decrease the discrepancy between CFD simulations and a real system. One thing is reflecting the change of pellet properties with time, and the other one, quite challenging, is the effect of fish on the flow field and possibly the particle (pellet) distribution in the tank. This might be neglected for small fish and a small number of fish, but the larger fish amount can substantially impact the system's hydrodynamics.