Turbidity currents are important carriers for transporting terrestrial sediment into the deep sea, facilitating the transfer of matter and energy between land and the deep sea. Previous studies have suggested that turbidity currents can exhibit high velocities during their movement in submarine canyons. However, the maximum vertical descent velocity of high-concentration turbid water simulating turbidity currents does not exceed 1 m/s, which does not support the understanding that turbidity currents can reach speeds of over twenty meters per second in submarine canyons. During their movement, turbidity currents can compress and push the water ahead, generating propagating waves. These waves, known as excitation waves, exert a force on the seafloor, resuspending bottom sediments and potentially leading to the generation of secondary turbidity currents downstream. Therefore, the propagation distance of excitation waves is not the same as the initial journey of the turbidity currents, and the velocity of excitation waves within this journey has been mistakenly regarded as the velocity of the turbidity currents. Research on the propagation velocity of excitation waves is of great significance for understanding the sediment supply patterns of turbidity currents and the transport patterns of deep-sea sediments. In this study, numerical simulations were conducted to investigate the velocity of excitation waves induced by turbidity currents and to explore the factors that can affect their propagation velocity and amplitude. The relationship between the velocity and amplitude of excitation waves and different influencing factors was determined. The results indicate that the propagation velocity of excitation waves induced by turbidity currents is primarily determined by the water depth, and an expression (v2 = 0.63gh) for the propagation velocity of excitation waves is provided.
Submarine turbidity currents, often referred to as underwater rivers, are important carriers that transport terrestrial sediments to the deep sea [1,2,3,4,5,6,7]. These turbidity currents, carrying a large amount of silt and sand, not only have strong erosive capabilities on the seabed [8,9,10], but also pose a threat to underwater communication cables, resulting in significant economic losses [11,12,13]. For example, the 2006 Pingdong earthquake in Taiwan caused the rupture of 11 submarine cables within the Kaoping Canyon, resulting in a slowdown in network speed in Southeast Asia for 49 days and requiring the deployment of 11 cable ships for repairs [13,14,15]. Investigating the velocity and patterns of turbidity currents in submarine canyons is of great significance for the protection of infrastructure such as pipelines and cables in these canyons. One of the main methods for quantitatively studying the velocity of turbidity currents in submarine canyons is to infer their speed through cable ruptures. The first confirmed occurrence of cable rupture caused by a turbidity current was in 1929, when the Grand Banks earthquake triggered the continuous rupture of 12 submarine cables. Inferred maximum turbidity current velocities reached 28 m/s [16,17,18]. Subsequently, multiple cable rupture incidents caused by turbidity currents have occurred worldwide. Table 1 summarizes the inferred maximum turbidity current velocities from these cable rupture incidents.
Event
Maximum Turbidity Velocity
References
18 November 1929 Grand Banks earthquake
28 m/s
[16,19,20,21]
1953 Suva earthquake in the Fiji Islands
5.1 m/s
[22]
The Orleansville earthquake of 9 September 1954, Algeria
20.6 m/s
[23]
Earthquake, Solomon Islands, Western Pacific, 23 December 1966
10.3 m/s
[24]
Incident at Nice airport, France, 16 October 1979
7 m/s
[25]
Taitung earthquake, 22 August 2002
9.8 m/s
[26]
21 May 2003 earthquake in Algeria
15.8 m/s
[27]
The Taitung earthquake of 10 December 2003
16.5 m/s
[26]
The Taitung earthquake of 18 December 2003
18.6 m/s
[26]
Pingtung earthquake on 26 December 2006
20 m/s
[28]
Typhoon Morakot on 7–9 August 2009
16.6 m/s
[29]
The 15 January 2022 eruption of Hunga volcano
33.9 m/s
[30]
Table 1. Cable breakage events caused by turbidity currents worldwide.
Previous studies have shown that the maximum vertical velocity of high-concentration turbidity currents in water does not exceed 1 m/s, and the maximum downward velocity of spherical particles in water does not exceed 10 m/s [31]. The maximum velocity of professional athlete Usain Bolt in the 100 m sprint on land is 9.58 m/s, while dolphins in the ocean can reach speeds of up to 20 m/s. Deep-sea turbidity currents, characterized by a small density difference compared to water, are primarily driven by the gravitational component along the direction of flow. However, factors such as bed friction also need to be considered. The driving force behind turbidity currents is primarily the density difference between the turbulent flow and the surrounding water, as well as the gravitational downslope component. Previous studies have detected a maximum sediment concentration of 12% in the basal layer of turbidity currents [32]. However, even high concentrations of suspended sediment, such as 1720 g/L, in seawater with a density of 1020 g/L, do not exceed a maximum vertical velocity of 1 m/s [33]. Similarly, spherical particles also have a maximum settling velocity in water of less than 10 m/s [33]. Turbidity currents, being density-driven flows, have relatively low density differences compared to water, and the gentle slope of submarine canyons also contributes to a smaller gravitational downslope force. Additionally, the influence of bed friction and other factors related to sediment deposition needs to be considered. It is incredible to think that turbidity currents can achieve flow velocities as high as 28 m/s [16,18,28,34,35]. When submarine landslides occur on continental slopes, the sliding mass entering the bottom of submarine canyons can cause the destruction of soft sediment beds. The mixing of sliding or flowing sediment with water forms turbidity currents. Turbidity currents exert pressure and propel the water ahead, forming an excitation wave. This aligns with Paull’s hypothesis that in the course of turbidity currents, a high-pressure zone is formed ahead, capable of causing an increase in pore water pressure in the sediment ahead [36]. Similar to surging waves, the excitation waves generated can propagate downstream along the submarine canyon, with a propagation velocity much greater than the velocity of turbidity currents [31]. The rapid propagation of excitation waves can exert a force on the seafloor of the submarine canyon, causing the resuspension of sediment in front of the head of the turbidity currents, which may lead to the formation of secondary turbidity currents at some downstream locations. The distance between the secondary and initial turbidity currents is actually the propagation distance of the excitation waves, rather than the journey of the initial turbidity currents. Therefore, the speed of the excitation waves within this distance is mistakenly considered as the velocity of the turbidity currents (see Figure 1). This may explain why the velocity of the turbidity currents as deduced from cable breakages is so high.
Figure 1. Diagram of excitation wave propagation due to turbidity current (v1 is the velocity of turbidity current. This refers to the ratio of distance to time experienced by a turbidity current mass moving underwater. v2 is the velocity of secondary turbidity current: the rapidly propagating excitation wave applies a force on the submarine canyon floor, leading to the destruction of the soft sediment floor and the secondary turbidity current. v is the propagation velocity of the excitation wave; this refers to the propagation velocity of the turbidity current excitation wave. This speed is not the velocity of the motion of the water mass. At time t0, the initial turbidity current moves underwater, pushing the stationary water in front to generate an excitation wave. At time t1, the excitation wave is propagating. At time t2, the rapidly propagating excitation wave exerts pressure on the soft bottom bed, resulting in the destruction of the bottom bed and secondary turbidity current).
Turbidity currents are mass movements composed of sediment particles, with a high concentration of the dense basal layer near the seabed. Depending on their density and granulometric composition, turbidity currents can move along submarine canyons through mechanisms such as diffusion, collapse, and flow [37], which differ from the downward movement as a single entity of landslide bodies after slope failure (this distinguishes them from surges). Additionally, during the long-distance movement of turbidity currents in canyons, the completion of subsequent water replenishment may generate multiple excitation waves. Furthermore, secondary excitation waves may also occur during the movement of secondary turbidity currents triggered by the initial turbidity current, which differs significantly from the surges caused by submarine landslides. Furthermore, previous studies [38,39,40,41] on sediment supply during turbidity current movements have mostly focused on the scouring action on the seabed, whereas the resuspension of sedimentary deposits in front of the initial turbidity current caused by excitation waves may serve as an effective mode of sediment supply during the long-distance transport of turbidity currents. In 2023, Ren et al. proposed that the cause of the long-distance high-speed motion of turbidity currents is due to the excitation waves caused by the primary turbidity currents. However, only preliminary research has been conducted on the comparison of excitation wave velocity and solitary wave velocity, and there has been no specific discussion on the reasons for the excitation wave velocity being much greater than that of the turbidity current. In an experiment conducted using an indoor flume, it was observed that the wavelength of the excitation waves was much larger than the water depth, similar to shallow water waves [33]. The amplitude of excitation waves in proportion to their wavelength was small, consistent with the theory of small-amplitude waves. Similar to the velocity model of shallow water waves, it is expected that the propagation speed of excitation waves is also influenced by the water depth. However, since excitation waves are triggered by sediment-laden turbidity currents, the velocity model may differ from that of surface waves induced by gravitational flows. The purpose of this study is to simulate and investigate the effects of different factors on the propagation velocity and amplitude of excitation waves through a validated numerical model based on laboratory experiments. The study aims to determine the maximum propagation velocity of excitation waves at a field scale and whether there is attenuation in the long-distance propagation after their formation. In recent studies, seafloor sediment flows have been collectively referred to as turbidity currents [42]. Therefore, we simulated the movement of turbidity currents by sediment flow. This study uses the CFD-based fluid computation software FLOW-3D to simulate the underwater movement process of turbidity currents. The numerical model is validated against indoor experimental results. During the simulation process, a velocity model for surging wave generation triggered by submarine landslides is used as a reference, and multiple factors that may affect the propagation velocity of the excitation wave are considered. By controlling a single variable, the main factors influencing the excitation wave propagation velocity are determined, and the corresponding expression for excitation wave propagation velocity is provided. The results indicate that the propagation velocity of the excitation wave induced by turbidity currents is primarily determined by the water depth. This research provides a new perspective for understanding the high-speed movement of turbidity currents in submarine canyons and enriches the understanding of the movement patterns of turbidity currents in submarine canyons. In addition, studying the propagation speed of excitation waves is highly significant for the resuspension of underwater sediments, as well as the re-circulation of carbon sequestration, nutrients, heavy metals, and microplastics.
2. Experimental Study on Excitation Waves Induced by Turbidity Currents
2.1. Experimental Design and Apparatus
The experimental apparatus used for the turbidity current-induced excitation wave tests is a straight water tank [33]. The water tank is 12.5 m long, 0.5 m wide, and 0.7 m high. A turbidity source area is located on the right side of the tank to generate turbidity currents. The tank is equipped with a terrain with a certain slope. Turbidity currents are generated underwater using a weir. The mass ratio of silt and clay used in the experimental turbid water solution was 8:2, with a density of 1600 kg/m3. Previous experiments have shown that this turbid mixture can reach a maximum flow velocity of 18.7 cm/s [31]. Three pressure sensors are placed along the straight section of the tank at intervals of 0.4 m. These sensors continuously monitor the bottom shear stress caused by the turbidity current-induced excitation wave, as well as the force exerted by the turbidity current itself on the bed. The monitoring frequency is set at 100 Hz.
2.2. Experimental Phenomenon and Results
In the laboratory water tank experiments, it was observed that as the turbidity current propagates, a wave is generated ahead of the turbidity front, moving in the same direction as the current and with a velocity greater than the turbidity current velocity [33]. By monitoring the pressure changes on the bed during the turbidity current motion [33], the propagation velocity of the excitation wave, the head movement velocity of the turbidity current, and the amplitude of the excitation wave (obtained from the measured surface elevation changes caused by the wave) can be estimated based on the distances between the sensors and the time when the pressure change peaks occur. The results of indoor experiments on turbidity currents indicate that they can compress and propel the water ahead of them, generating excitation waves similar to pulses. The propagation speed of these excitation waves caused by turbidity currents is found to be much greater than the velocity of the turbidity current movement at its head, as determined by pressure sensors installed on the seabed.
3. Numerical Simulation of Excitation Waves Induced by Turbidity Currents
FLOW-3D is a powerful computational fluid dynamics (CFD) software that excels in making accurate calculations of free surface and six-degrees-of-freedom motions of objects. Similar to other CFD software, FLOW-3D consists of three modules: pre-processing, solver, and post-processing. In recent years, there have been many simulations of turbidity currents using FLOW-3D due to its superior capabilities. For example, Heimsund (2007) simulated turbidity currents in the Monterey Canyon system using FLOW-3D based on high-resolution bathymetry and flow data [43]. Zhou et al. (2017) used FLOW-3D software to simulate turbidity currents in a flume with obstacles, analyzing the impact of the proportion between obstacle height and flume height on the movement of turbidity currents, including their velocity, flow state, and morphological evolution [44]. In this study, using the CFD software FLOW-3D, the underwater motion process of turbidity currents is simulated. The model is validated by comparing it with experimental results, and the motion of the waves induced by turbidity currents is simulated based on this validation.
3.1. Control Equations
FLOW-3D, a mature three-dimensional fluid simulation software, is used in this study. It employs the RNG turbulence model, which is capable of handling high strain rate flows and is suitable for simulating excitation waves. The research focus of this paper is on sediment gravity flows (turbulent flows), and the control equations used in the calculations include the basic continuity equation, the momentum equation, the turbulent kinetic energy k equation, and the turbulent kinetic energy dissipation rate ε equation.
The continuity equation:
The momentum equation:
The turbulence model:
k equation:
ε equation:
where u, v and w is the flow velocity component in x, y and z directions; Ax, Ay and Az represent the area fraction that can flow in x, y and z directions; Gx, Gy and Gz are the gravitational acceleration in x, y and z directions; fx, fy and fz are the viscous forces in the three directions; VF is the fraction of the volume that can flow; ρ is the fluid density; p is the pressure acting on the fluid element; k is the turbulence energy; ε is the turbulence kinetic energy dissipation rate; μ is turbulence viscosity coefficient
where u, v and w is the flow velocity component in x, y and z directions; Ax, Ay and Az represent the area fraction that can flow in x, y and z directions; Gx, Gy and Gz are the gravitational acceleration in x, y and z directions; fx, fy and fz are the viscous forces in the three directions; VF is the fraction of the volume that can flow; ρ is the fluid density; p is the pressure acting on the fluid element; k is the turbulence energy; ε is the turbulence kinetic energy dissipation rate;
μ is turbulence viscosity coefficient where Cμ = 0.0845;
Gk is the turbulent kinetic energy generation term, expressed as
and σk and σε are the Prandtl numbers corresponding to the turbulent kinetic energy and dissipation rate, respectively, both of which are 1.39.
In addition,where Cε1 and Cε2 are the empirical constants, 1.42 and 1.68, respectively.
Furthermore,
where Eij=12∂ui∂xj+∂uj∂xi, η0 = 4.377, β = 0.012.
The general mass continuity equation is as follows:
where VF is the fractional volume open to flow, ρ is the fluid density, RDIF is a turbulent diffusion term, and RSOR is the mass source.
3.2. Model Validation
To determine the factors affecting the velocity of the turbidity-induced excitation wave and its velocity expression, first, the indoor flume test was taken as the prototype. Then, a 1:1 geometric solid model was established, and the simulation parameters were set to be consistent with the flume test parameters [33]. Finally, the simulation results were compared with the laboratory test results.
The computational domain employs the method of unstructured grid and is entirely divided into structured orthogonal grids. Nested grids are used for local refinement at the interfaces of straight sections, resulting in a total of 800,000 grid cells after refinement.
The simulation results were compared with the indoor experimental results, with the velocity of the excitation wave and the turbidity current head being represented by changes in surface elevation and water density. The experimental and simulation results are shown in Table 2, and the calculation formula for the error is |Calculated value−Test value|Test value×100%Calculated value-Test valueTest value×100%.
Result
Propagation Velocity of Excitation Wave (m/s)
Velocity of Turbidity Current (m/s)
Excitation Wave Amplitude (m)
Sensor 1 to 2
Sensor 2 to 3
Sensor 1 to 2
Sensor 2 to 3
Sensor 1 to 2
Sensor 2 to 3
Test results
1.54
1.48
0.24
0.23
0.029
0.03
Computed results
1.55
1.52
0.25
0.23
0.03
0.03
Error range
0.6%
2.7%
4.2%
0%
3.4%
0%
Table 2. The test results of the propagation velocity of the excitation wave, the turbidity current velocity, and the excitation wave amplitude are compared with the simulation results.
From the above comparison, it can be observed that the simulated velocities of the excitation wave and the head of the turbidity current align well with the experimental results, indicating the rationality of using the numerical model established in this study for simulating the propagation velocity of the excitation wave induced by turbidity currents.
3.3. Analysis of Factors Affecting the Propagation Velocity of Excitation Waves
An analysis of the factors influencing the propagation velocity of excitation waves was conducted using numerical simulation. The reference model for wave velocity was based on the surge velocity model. The main factors affecting the propagation velocity of excitation waves were summarized, including the turbidity current density ρ, the thickness of the turbidity current source area d, the length of the turbidity current source area L, the depth at the initial flow of turbidity currents h, the canyon width l, and the initial velocity of the turbidity current v0 (as shown in Figure 2). The simulations were performed using a controlled variable approach for different parameters, and the velocity changes of the excitation wave were obtained, as shown in Table 3. The slope angle was fixed at 3°, and sensors were placed at intervals of 100 m starting from a distance of 500 m from the turbidity current source area (named Sensors 1, 2, 3). These sensors were used to extract surface elevation, density, and other relevant parameters at their respective locations. We can obtain the propagating velocity of excitation waves by measuring the time difference in surface elevation changes at the monitoring points. Similarly, we can determine the propagation velocity of turbidity currents by measuring the time difference in density changes.
Figure 2. Excitation wave velocity simulation model and parameters.
Group Order
Turbidity Current Density (kg/m3)
Length of Turbidity Source Area (m)
Canyon Width (m)
Thickness of Turbidity Source Area (m)
Depth (m)
Initial Velocity of Turbidity Current (m/s)
Propagation Velocity of Excitation Wave (m/s)
Excitation Wave Amplitude (m)
Velocity of Turbidity Current (m/s)
1
1600
1000
200
20
200
0
33.43
0.345
5.88
2
1500
1000
200
20
200
0
33.09
0.304
5.41
3
1400
1000
200
20
200
0
33.35
0.223
4.99
4
1300
1000
200
20
200
0
33.33
0.177
4.35
5
1200
1000
200
20
200
0
33.86
0.092
3.74
6
1600
1000
200
40
200
0
33.05
1.109
9.09
7
1600
1000
200
60
200
0
33.39
2.689
10.79
8
1600
1000
200
80
200
0
33.21
4.828
12.91
9
1600
1000
200
100
200
0
36.43
7.744
13.79
10
1600
200
200
20
200
0
32.93
0.181
5.58
11
1600
400
200
20
200
0
33.49
0.25
5.71
12
1600
600
200
20
200
0
33.06
0.278
5.79
13
1600
800
200
20
200
0
33.17
0.31
5.72
14
1600
1000
200
20
100
0
26.67
0.56
5.72
15
1600
1000
200
20
300
0
39.65
0.169
5.80
16
1600
1000
200
20
400
0
45.98
0.12
5.80
17
1600
1000
200
20
500
0
49.97
0.08
5.96
18
1600
1000
100
20
200
0
33.60
0.354
5.72
19
1600
1000
300
20
200
0
32.98
0.338
5.97
20
1600
1000
400
20
200
0
33.27
0.356
5.87
21
1600
1000
500
20
200
0
33.31
0.365
5.86
22
1600
1000
200
20
200
2
33.50
0.532
4.35
23
1600
1000
200
20
200
5
33.12
1.389
6.56
24
1600
1000
200
20
200
8
33.52
2.271
8.10
25
1600
1000
200
20
200
10
33.33
2.878
8.99
Table 3. Simulation results under different variables conditions.
The variations in surface elevation at three sensor locations in the simulated results of five different turbidity current density groups are presented in Figure 3.
Figure 3. Simulation of propagating velocity of excitation wave under the sole variable condition of turbulent current density. (Length of turbidity source area: 1000 m; canyon width: 200 m; thickness of turbidity source area: 20 m; depth: 200 m; initial velocity of turbidity current: 0 m/s).
Based on the simulation results described above, while keeping all other conditions constant, the impact of a single variable, namely, the turbidity current density, on the propagation velocity and amplitude of the excitation wave was analyzed. By fitting the data, the relationship between turbidity current density and the propagation velocity of turbidity currents as well as the amplitude of the excitation wave was obtained, as shown in Figure 4.
Figure 4. Relationship between turbidity current density and turbidity current velocity, as well as excitation wave amplitude.
The simulation results indicate that changes in turbidity current density, while keeping the other conditions constant, do not result in a change in the propagation velocity of the excitation waves. However, they do affect the amplitude of the excitation waves and the velocity of the turbidity current itself. The simulation reveals that within the selected density range, both the amplitude of the excitation waves and the velocity of the turbidity current increase with increasing turbidity current density. When the turbidity current density is equal to that of water (ρTurbidity current = ρWater), there is no turbidity current or excitation wave generation. Thus, the relationship between the turbidity current velocity (v) and density (ρ) is expressed as v = −34.80643 + 0.05082•ρ − 1.59286 × 10−5ρ2 (ρ > 1000, R2 = 0.994). Additionally, the relationship between the amplitude of the excitation waves (A) caused by turbidity currents and density (ρ) is expressed as A = −0.6021 + 5.9729 × 10−4ρ (ρ > 1000, R2 = 0.991).
3.3.2. The Influence of the Thickness of the Turbidity Source Area on the Propagation Velocity and Amplitude of Excitation Waves
The variations in surface elevation at three sensor locations in the simulated results of five different thickness of turbidity source area groups are presented in Figure 5.
Figure 5. Simulation of propagating velocity of excitation wave under the sole variable condition of thickness of turbidity source area. (Turbidity current density: 1600 kg/m3; length of turbidity source area: 1000 m; canyon width: 200 m; depth: 200 m; initial velocity of turbidity current: 0 m/s).
Based on the simulation results described above, while keeping all other conditions constant, the impact of a single variable, namely, the thickness of the turbidity source area, on the propagation velocity and amplitude of the excitation wave was analyzed. By fitting the data, the relationship between the thickness of the turbidity source area and the propagation velocity of the turbidity current as well as the amplitude of the excitation wave was obtained, as shown in Figure 6.
Figure 6. Relationship between thickness of turbidity source area and turbidity current velocity, as well as excitation wave amplitude.
Based on the simulated results mentioned above, it can be concluded that, while keeping the other conditions constant, changing only the thickness of the turbidity current source area does not affect the propagation velocity of the excitation waves. However, it does impact both the amplitude of the excitation waves and the velocity of the turbidity current itself. The simulation reveals that within the selected range of thickness values for the turbidity current source area, both the amplitude of the excitation waves and the velocity of the turbidity current increase with an increase in the thickness of the source area. Additionally, it is observed that when the length of the turbidity current source area is zero, neither the turbidity current nor the excitation waves are generated (i.e., no turbidity current is produced when hTurbidity current = 0). Therefore, the relationship between the velocity (v) of the turbidity current and its thickness (h) is expressed as v = 0.27983•h − 0.00146•h2 (h ≥ 0, R2 = 0.999). Similarly, the relationship between the amplitude (A) of the excitation waves caused by the turbidity current and its thickness (h) is A = −0.00375•h − 0.0008•h2 (h ≥ 0, R2 = 0.999).
3.3.3. The Influence of the Length of the Turbidity Source Area on the Propagation Velocity and Amplitude of Excitation Waves
The variations in surface elevation at three sensor locations in the simulated results of five different length of turbidity source area groups are presented in Figure 7.
Figure 7. Simulation of propagating velocity of excitation wave under the sole variable condition of length of turbidity source area. (Turbidity current density: 1600 kg/m3; canyon width: 200 m; thickness of turbidity source area: 20 m; depth: 200 m; initial velocity of turbidity current: 0 m/s).
Based on the simulation results described above, while keeping all other conditions constant, the impact of a single variable, namely, the length of the turbidity source area, on the propagation velocity and amplitude of the excitation wave was analyzed. By fitting the data, the relationship between the length of the turbidity source area and the amplitude of the excitation wave was obtained, as shown in Figure 8.
Figure 8. Relationship between length of turbidity source area and excitation wave amplitude.(Amplitude refers to the surface elevation change caused by the excitation wave).
Through simulations, it has been determined that within the chosen range of the length of the turbidity source area, the amplitude of the excitation waves increases with an increase in the length of the turbidity source area. When the length of the turbidity source area is zero, there is no turbidity current and no generation of excitation waves (i.e., when LTurbidity current = 0). Additionally, for large lengths of the turbidity source area, under the condition of sufficient sediment supply, the variations in surface elevation caused by the waves generated by turbidity currents are negligible. Therefore, the relationship between the amplitude of the excitation waves (A) generated by turbidity currents and the length of the turbidity source area (L) is expressed as follows: A = −0.3624 + 0.10305•ln(L − 6.15619) (L ≥ 0, R2 = 0.997).
3.3.4. The Influence of Depth on the Propagation Velocity and Amplitude of Excitation Waves
The variations in surface elevation at three sensor locations in the simulated results of five different depth groups are presented in Figure 9.
Figure 9. Simulation of propagation velocity of excitation wave under the sole variable condition of depth. (Turbidity current density: 1600 kg/m3; length of turbidity source area: 1000 m; canyon width: 200 m; thickness of turbidity source area: 20 m; initial velocity of turbidity current: 0 m/s).
Based on the simulation results described above, while keeping all other conditions constant, the impact of a single variable, namely, depth, on the propagation velocity and amplitude of the excitation wave was analyzed. By fitting the data, the relationship between depth and the propagation velocity of the excitation wave as well as the amplitude of the excitation wave was obtained, as shown in Figure 10.
Figure 10. Relationship between depth and propagating velocity of excitation wave, as well as excitation wave amplitude.
As the water depth approaches infinity, the excitation wave amplitude can only approach zero but cannot reach zero. Therefore, the characteristics of the excitation wave amplitude change with the water depth are similar to those of the velocity propagation of the excitation wave. The relationship between the velocity of the excitation wave induced by turbidity currents (vExcitation wave) and the water depth (H) can be described as vExcitation wave = −287.05446 + 48.59211•ln(H + 535.14863) (R2 = 0.998). The relationship between the excitation wave amplitude (A) and the water depth (H) can be expressed as A = 1.46573 − 0.22816•ln(H − 47.67563) (R2 = 0.985).
3.3.5. The Influence of the Canyon Width on the Propagation Velocity and Amplitude of Excitation Waves
The variations in surface elevation at three sensor locations in the simulated results of five different canyon width groups are presented in Figure 11.
Figure 11. Simulation of propagating velocity of excitation wave under the sole variable condition of canyon width. (Turbidity current density: 1600 kg/m3; length of turbidity source area: 1000 m; thickness of turbidity source area: 20 m; depth: 200 m; initial velocity of turbidity current: 0 m/s).
When the canyon width is taken as the single variable condition, changing the canyon width does not significantly affect the propagation velocity of excitation waves, the amplitude of excitation waves, and the velocity of turbidity currents. Therefore, it can be concluded that, without considering the impact of the differences in the terrain and sediment on the canyon width, the canyon width has no impact on the propagation of excitation waves and the movement of turbidity currents.
3.3.6. The Influence of the Initial Velocity of the Turbidity Current on the Propagation Velocity and Amplitude of Excitation Waves
The variations in surface elevation at three sensor locations in the simulated results of five different initial velocity of turbidity current groups are presented in Figure 12.
Figure 12. Simulation of propagating velocity of excitation wave under the sole variable condition of initial velocity of turbidity current. (Turbidity current density: 1600 kg/m3; length of turbidity source area: 1000 m; canyon width: 200 m; thickness of turbidity source area: 20 m; depth: 200 m).
Based on the simulation results described above, while keeping all other conditions constant, the impact of a single variable, namely, the initial velocity of the turbidity current, on the propagation velocity and amplitude of the excitation wave was analyzed. By fitting the data, the relationship between the initial velocity of the turbidity current and the amplitude of the excitation wave was obtained, as shown in Figure 13.
Figure 13. Relationship between initial velocity of turbidity current and excitation wave amplitude.
Based on the simulation, it is observed that within the selected range of the initial velocity of the turbidity current, the amplitude of the excitation wave increases linearly with the increase in the initial velocity of the turbidity current. Therefore, the relationship between the amplitude (A) of the excitation wave caused by the turbidity current and the initial velocity of the turbidity current (v0) can be expressed as A = 0.34 + 0.24084•v0 (A ≥ 0, R2 = 0.992).
Through controlling the simulation calculation of a single variable, it was found that there are several factors that can affect the amplitude of the excitation wave. These factors include the turbidity current density ρ, the thickness of the turbidity current source area d, the length of the turbidity current source area L, the water depth h, and the initial velocity of the turbidity current v0. In contrast, there are relatively few factors that influence the propagation velocity of the excitation wave. Within the selected parameter range, only the water depth can affect the propagation velocity of the excitation wave. The physical parameters of the turbidity current, including the turbidity current density ρ, the thickness of the turbidity current source area d, the length of the turbidity current source area L, the canyon width l, and the initial velocity of the turbidity current v0, have no direct influence on the propagation velocity of the excitation wave. Therefore, the turbidity current only serves as a triggering factor for the excitation wave and is not directly related to the propagation velocity of the excitation wave.
3.4. Analyze the Changes in Propagation Velocity of Excitation Waves along a Path
In order to further investigate the underlying truth behind the variation in the propagation velocity of the excitation wave, a discussion on whether there is velocity attenuation along the propagation path of the excitation wave is conducted. Since the seventh group of the excitation wave causes significant changes in surface elevation, the seventh group of the excitation wave is selected as the research object in order to study the variations in surface elevation along the propagation path of the excitation wave. The changes in surface elevation are extracted every 200 m along the sediment slope (with the first extraction point located 400 m away from the source area of the turbidity current). A total of six sets of surface elevation data are extracted (ranging from 400 m to 1400 m distance from the source area of the turbidity current), as shown in Figure 14.
Figure 14. Surface elevation changes during excitation wave propagation along sediment slopes.
The amplitudes and propagation velocities of the excitation wave at each point are shown in Table 4.
Distance from Turbidity Current Source Area (m)
Propagation Velocity of Excitation Wave (m/s)
Excitation Wave Amplitude (m)
400
33.34
2.524
600
36.79
2.596
800
37.13
2.589
1000
39.99
2.566
1200
40.04
2.542
1400
40.13
2.523
Table 4. Excitation wave velocity during the excitation wave propagation along the sediment slope.
From the table above, it can be observed that the amplitude of the excitation wave does not change while traveling along the slope. This indicates that the change in surface elevation caused by the propagation of the excitation wave does not attenuate. Furthermore, the propagation velocity of the excitation wave gradually increases, although the change is not very pronounced. This variation may be attributed to the change in the water depth caused by the sloping bed. To investigate this, a simulation was conducted in a straight channel with a length of 3000 m. Six sampling points were established from 400 m to 1400 m away from the turbidity current source area to extract the amplitude of the excitation wave. The results of the simulation are presented in Figure 15.
Figure 15. Surface elevation changes during wave propagation along a straight channel.
The amplitudes and propagation velocities of the excitation wave at each point are shown in Table 5.
Distance from Turbidity Current Source Area (m)
Propagation Velocity of Excitation Wave (m/s)
Excitation Wave Amplitude (m)
400
33.89
2.559
600
37.66
2.692
800
37.12
2.712
1000
36.92
2.717
1200
37.09
2.715
1400
37.48
2.718
Table 5. Excitation wave velocity during the propagation along the straight channel.
The data from the table above indicate that during the propagation of the excitation wave along a straight water channel, its velocity remains constant, except for a slight decrease at the initial point. This phenomenon may be attributed to the fact that in the starting phase, the excitation wave is not fully developed, and hence its velocity is relatively smaller. However, once it is fully developed, the propagation velocity of the excitation wave does not decrease in subsequent processes. Therefore, the propagation velocity of the excitation wave is only dependent on the real-time water depth of the wave. In future studies, we aim to explore the relationships between these influencing factors and other physical parameters, such as the speed of wave propagation, using the effective and accurate method of machine learning algorithms [45].
3.5. Expression of the Propagation Velocity of the Excitation Wave
The propagation of the excitation wave along a long distance does not experience an attenuation in velocity, as is the case with the propagation velocity of solitary waves. Referring to the estimated wave propagation velocity (the square of the propagation velocity is directly proportional to the water depth amplitude) [46], the wavelengths under different water depth conditions were extracted, as shown in Table 6.
Depth (m)
Propagation Velocity of Excitation Wave (m/s)
Excitation Wave Amplitude (m)
Excitation Wave Length (m)
100
26.67
0.56
2580
200
33.43
0.35
2850
300
39.65
0.17
3250
400
45.98
0.12
3600
500
49.97
0.08
4150
1000
66.67
0.04
6000
2000
90.91
0.02
9500
4000
165.84
0.3
17600
Table 6. Physical parameters of excitation wave under different water depth conditions.
From the simulation results of a single variable, the water depth, it could be seen that the wavelengths of the excitation waves were much larger than the water depth. Therefore, further simulations were conducted under water depth conditions ranging from 1000 m to 4000 m. Due to the minimal change in wave amplitude when the water depth reached 4000 m, it was not possible to observe a distinct waveform. However, through simulations with the thickness of the turbidity current source area as the single variable, it was found that an increase in the thickness of the source region led to a larger amplitude of the excitation waves, but it did not affect the wavelength of the excitation waves. Therefore, in order to better extract the wavelength of the excitation waves, the thickness of the source region in the simulation with a water depth of 4000 m was set to 200 m.
Through simulations at water depths of 1000 m and 4000 m, it is observed that the wavelengths of the excitation waves are much larger than the water depth, indicating that these waves belong to the category of shallow water waves. The amplitude of the excitation waves is relatively small compared to their wavelength, aligning with the small amplitude wave theory [47]. According to this theory, the wave velocity of shallow water waves is only dependent on the water depth (h) and gravity acceleration (g), regardless of the wave period. In the case of excitation waves induced by turbidity currents in deep water, the amplitudes of these waves are relatively small compared to the water depth. Referring to the expression for shallow water waves (when the relative water depth, which is the ratio of water depth to wavelength, is much smaller than 1/2), the wave velocity is denoted as 𝐶𝑠=√𝑔ℎ. This implies that the propagation velocity of the excitation waves is also solely related to the water depth. Therefore, a fitting of the square of the propagation velocity of the excitation waves (v2) and the water depth (h) was conducted (Figure 16).
Figure 16. The relationship between the propagation velocity of excitation wave and the depth.
Through fitting, the following can be obtained:
Through fitting, it can be discovered that the propagation model of the velocity of excitation waves is different from the shallow water wave theory. This is because turbidity currents, as granular materials, generate excitation waves by pushing the water in front of them with sediment particles underwater, which is different from the surges formed by solid blocks entering the ocean. Additionally, excitation waves formed by turbidity currents occur in an underwater environment, which may be the reason why the propagation velocity equation for the excitation waves behaves as if the velocity squared is equal to half the Earth’s gravity. This equation reveals the variation in the propagation velocity of the excitation wave with depth, explaining why the average velocity between the monitoring points in the field is greater than the instantaneous velocity measured at these points [41]. Further theoretical research on the propagation velocity of excitation waves requires subsequent field monitoring and the deployment of monitoring systems to more thoroughly investigate the fundamental causes.
4. Conclusions
This study aimed to investigate the velocity of turbidity current-induced excitation waves through numerical simulation. By fixing a single variable, different factors that could affect the propagation velocity and amplitude of the excitation waves were analyzed and discussed, leading to the following three conclusions:
Within the selected parameter range, there are several factors that can influence the amplitude of the excitation waves, including the turbidity current density ρ, the thickness of the turbidity current source area d, the length of the turbidity current source area L, the water depth h, and the initial velocity of the turbidity current v0.The amplitude of the excitation waves is positively correlated with the turbidity density, the thickness of the source area, the length of the source area, and the initial velocity, while it is negatively correlated with the water depth.
Within the selected parameter range, only the water depth can affect the propagation velocity of the excitation waves. As the water depth increases, the propagation velocity of the excitation waves also increases, and a relationship of v2 = 0.63gh (R2 = 0.967) is established between the square of the propagation velocity v2 and the water depth h.
During the propagation of the excitation waves, both the propagation velocity and the changes in surface elevation caused by the waves do not attenuate. Considering the relatively calm deep-sea environment, the high-speed propagation of the excitation waves and the resuspension of bottom sediments they cause not only complement the understanding of turbidity current motion patterns in canyons, but also provide new research directions for deep-sea sediment transport.
References
Azpiroz-Zabala, M.; Cartigny, M.J.B.; Talling, P.J.; Parsons, D.R.; Sumner, E.J.; Clare, M.A.; Simmons, S.M.; Cooper, C.; Pope, E.L. Newly recognized turbidity current structure can explain prolonged flushing of submarine canyons. Sci. Adv.2017, 3, e1700200.
Daly, R.A. Origin of submarine canyons. Am. J. Sci.1936, 5, 401–420.
Winterwerp, J. Stratification effects by fine suspended sediment at low, medium, and very high concentrations. J. Geophys. Res. Oceans.2006, 111, C5.
Nilsen, T.H.; Shew, R.D.; Steffens, G.S.; Studlick, J.R.J. Atlas of Deep-Water Outcrops; American Association of Petroleum Geologists: Tulsa, OK, USA, 2008.
Xu, J. Turbidity Current Research in the Past Century: An Overview. J. Ocean Univ. China2014, 44, 98–105.
Talling, P.J.; Allin, J.; Armitage, D.A.; Arnott, R.W.C.; Cartigny, M.J.B.; Clare, M.A.; Felletti, F.; Covault, J.A.; Girardclos, S.; Hansen, E.; et al. Key future directions for research on turbidity currents and their deposits. J. Sediment. Res.2015, 85, 153–169.
Maier, K.L.; Gales, J.A.; Paull, C.K.; Rosenberger, K.; Talling, P.J.; Simmons, S.M.; Gwiazda, R.; McGann, M.; Cartigny, M.J.; Lundsten, E. Linking direct measurements of turbidity currents to submarine canyon-floor deposits. Front. Earth Sci.2019, 7, 144.
Hughes Clarke, J.E. First wide-angle view of channelized turbidity currents links migrating cyclic steps to flow characteristics. Nat. Commun.2016, 7, 11896.
Normandeau, A.; Bourgault, D.; Neumeier, U.; Lajeunesse, P.; St-Onge, G.; Gostiaux, L.; Chavanne, C. Storm-induced turbidity currents on a sediment-starved shelf: Insight from direct monitoring and repeat seabed mapping of upslope migrating bedforms. Sedimentology2020, 67, 1045–1068.
Hill, P.R.; Lintern, D.G. Turbidity currents on the open slope of the Fraser Delta. Mar. Geol.2022, 445, 106738.
Summers, M. Review of deep-water submarine cable design. In Proceedings of the SubOptic 2001, Kyoto, Japan, 20–24 May 2001; p. 4.
Carter, L.; Gavey, R.; Talling, P.J.; Liu, J.T. Insights into Submarine Geohazards from Breaks in Subsea Telecommunication Cables. Oceanography2014, 27, 58–67.
Gavey, R.; Carter, L.; Liu, J.T.; Talling, P.J.; Hsu, R.; Pope, E.; Evans, G. Frequent sediment density flows during 2006 to 2015, triggered by competing seismic and weather events: Observations from subsea cable breaks off southern Taiwan. Mar. Geol.2017, 384, 147–158.
Carter, L.; Burnett, D.; Drew, S.; Hagadorn, L.; Marle, G.; Bartlett-Mcneil, D.; Irvine, N. Submarine Cables and the Oceans: Connecting the World; UNEP-WCMC Biodiversity Series 31; UNEP World Conservation Monitoring Centre: Cambridge, UK, 2010.
Qiu, W. Submarine cables cut after Taiwan earthquake in Dec 2006. Submar. Cable Netw.2011, 19.
Heezen, B.C.; Ewing, W.M. Turbidity currents and submarine slumps, and the 1929 Grand Banks [Newfoundland] earthquake. Am. J. Sci.1952, 250, 849–873.
Kuenen, P.H. Estimated size of the Grand Banks [Newfoundland] turbidity current. Am. J. Sci.1952, 250, 874–884.
Heezen, B.C.; Ericson, D.; Ewing, M. Further evidence for a turbidity current following the 1929 Grand Banks earthquake. Deep-Sea Res.1954, 1, 193–202.
Heezen, B.C. Whales entangled in deep sea cables. Deep-Sea Res.1957, 4, 105–115.
Piper, D.J.; Shor, A.N.; Farre, J.A.; O’Connell, S.; Jacobi, R. Sediment slides and turbidity currents on the Laurentian Fan: Sidescan sonar investigations near the epicenter of the 1929 Grand Banks earthquake. Geology1985, 13, 538–541.
Piper, D.J.; Cochonat, P.; Morrison, M.L. The sequence of events around the epicentre of the 1929 Grand Banks earthquake: Initiation of debris flows and turbidity current inferred from sidescan sonar. Sedimentology1999, 46, 79–97.
Houtz, R.; Wellman, H. Turbidity current at Kadavu Passage, Fiji. Geol. Mag.1962, 99, 57–62.
Heezen, B.C.; Ewing, M. Orleansville earthquake and turbidity currents. AAPG Bull.1955, 39, 2505–2514.
Krause, D.C.; White, W.C.; PIPER, D.J.W.; Heezen, B.C. Turbidity currents and cable breaks in the western New Britain Trench. Geol. Soc. Am. Bull.1970, 81, 2153–2160.
Piper, D.J.; Savoye, B. Processes of late Quaternary turbidity current flow and deposition on the Var deep-sea fan, north-west Mediterranean Sea. Sedimentology1993, 40, 557–582.
Soh, W.; Machiyama, H.; Shirasaki, Y.; Kasahara, J. Deep-sea floor instability as a cause of deepwater cable fault, off eastern part of Taiwan. AGU Fall Meet. Abstr.2004, 2, 1–8.
Cattaneo, A.; Babonneau, N.; Ratzov, G.; Dan-Unterseh, G.; Yelles, K.; Bracène, R.; De Lepinay, B.M.; Boudiaf, A.; Déverchère, J. Searching for the seafloor signature of the 21 May 2003 Boumerdès earthquake offshore central Algeria. Nat. Hazards Earth Syst. Sci.2012, 12, 2159–2172.
Hsu, S.-K.; Kuo, J.; Chung-Liang, L.; Ching-Hui, T.; Doo, W.-B.; Ku, C.-Y.; Sibuet, J.-C. Turbidity currents, submarine landslides and the 2006 Pingtung earthquake off SW Taiwan. Terr. Atmos. Ocean. Sci.2008, 19, 7.
Carter, L.; Milliman, J.D.; Talling, P.J.; Gavey, R.; Wynn, R.B. Near-synchronous and delayed initiation of long run-out submarine sediment flows from a record-breaking river flood, offshore Taiwan. Geophys. Res. Lett.2012, 39, L12603.
Clare, M.A.; Yeo, I.A.; Watson, S.; Wysoczanski, R.; Seabrook, S.; Mackay, K.; Hunt, J.E.; Lane, E.; Talling, P.J.; Pope, E.; et al. Fast and destructive density currents created by ocean-entering volcanic eruptions. Science2023, 381, 1085–1092.
Ren, Y.; Zhang, Y.; Xu, G.; Xu, X.; Wang, H.; Chen, Z. The failure propagation of weakly stable sediment: A reason for the formation of high-velocity turbidity currents in submarine canyons. J. Ocean. Limnol.2023, 41, 100–117.
Wang, Z.; Xu, J.; Talling, P.J.; Cartigny, M.J.B.; Simmons, S.M.; Gwiazda, R.; Paull, C.K.; Maier, K.L.; Parsons, D.R. Direct evidence of a high-concentration basal layer in a submarine turbidity current. Deep-Sea Res. Part I Oceanogr. Res. Pap.2020, 161, 103300.
Ren, Y.; Tian, H.; Chen, Z.; Xu, G.; Liu, L.; Li, Y. Two Kinds of Waves Causing the Resuspension of Deep-Sea Sediments: Excitation and Internal Solitary Waves. J. Ocean Univ. China2023, 22, 429–440.
Lambert, A.M.; Kelts, K.R.; Marshall, N.F. Measurements of density underflows from Walensee, Switzerland. Sedimentology1976, 23, 87–105.
Piper, D.J.W.; Shor, A.N.; Hughes Clarke, J.E. The 1929 “Grand Banks” earthquake, slump, and turbidity current. In Sedimentologic Consequences of Convulsive Geologic Events; Geological Society of America: Boulder, CO, USA, 1988; pp. 77–92.
Ren, Y.; Zhou, H.; Wang, H.; Wu, X.; Xu, G.; Meng, Q. Study on the critical sediment concentration determining the optimal transport capability of submarine sediment flows with different particle size composition. Mar. Geol.2023, 464, 107142.
Bagnold, R.A. Auto-suspension of transported sediment; turbidity currents. Proc. R. Soc. Lond. Ser. A Math. Phys. Sci.1962, 265, 315–319.
Parker, G. Conditions for the ignition of catastrophically erosive turbidity currents. Mar. Geol.2003, 46, 307–327.
Pantin, H.M. Interaction between velocity and effective density in turbidity flow: Phase-plane analysis, with criteria for autosuspension. Mar. Geol.1979, 31, 59–99.
Heerema, C.J.; Talling, P.J.; Cartigny, M.J.; Paull, C.K.; Bailey, L.; Simmons, S.M.; Parsons, D.R.; Clare, M.A.; Gwiazda, R.; Lundsten, E.; et al. What determines the downstream evolution of turbidity currents? Earth Planet. Sci. Lett.2020, 532, 116023.
Talling, P.J.; Cartigny, M.J.B.; Pope, E.; Baker, M.; Clare, M.A.; Heijnen, M.; Hage, S.; Parsons, D.R.; Simmons, S.M.; Paull, C.K.; et al. Detailed monitoring reveals the nature of submarine turbidity currents. Nat. Rev. Earth Environ.2023, 4, 642–658.
Heimsund, S. Numerical Simulation of Turbidity Currents: A New Perspective for Small-and Large-Scale Sedimentological Experiments; Sedimentology/Petroleum Geology; University of Bergen: Bergen, Norway, 2007.
Zhou, J.; Cenedese, C.; Williams, T.; Ball, M.; Venayagamoorthy, S.K.; Nokes, R.I. On the Propagation of Gravity Currents Over and through a Submerged Array of Circular Cylinders. J. Fluid Mech.2017, 831, 394–417.
Saha, S.; De, S.; Changdar, S. An Application of Machine Learning Algorithms on the Prediction of the Damage Level of Rubble-Mound Breakwaters. J. Offshore Mech. Arct. Eng.2024, 146, 011202.
Russell, J.S. Report on Waves: Made to the Meetings of the British Association; Richard and John E Taylor: London, UK, 1845.
Airy, G.B. Tides and Waves; B. Fellowes: London, UK, 1845.
Review on Blood Flow Dynamics in Lab-on-a-Chip Systems: An Engineering Perspective
Bin-Jie Lai
,
Li-Tao Zhu
,
Zhe Chen*
,
Bo Ouyang*
, and
Zheng-Hong Luo*
Abstract
다양한 수송 메커니즘 하에서, “LOC(lab-on-a-chip)” 시스템에서 유동 전단 속도 조건과 밀접한 관련이 있는 혈류 역학은 다양한 수송 현상을 초래하는 것으로 밝혀졌습니다.
본 연구는 적혈구의 동적 혈액 점도 및 탄성 거동과 같은 점탄성 특성의 역할을 통해 LOC 시스템의 혈류 패턴을 조사합니다. 모세관 및 전기삼투압의 주요 매개변수를 통해 LOC 시스템의 혈액 수송 현상에 대한 연구는 실험적, 이론적 및 수많은 수치적 접근 방식을 통해 제공됩니다.
전기 삼투압 점탄성 흐름에 의해 유발되는 교란은 특히 향후 연구 기회를 위해 혈액 및 기타 점탄성 유체를 취급하는 LOC 장치의 혼합 및 분리 기능 향상에 논의되고 적용됩니다. 또한, 본 연구는 보다 정확하고 단순화된 혈류 모델에 대한 요구와 전기역학 효과 하에서 점탄성 유체 흐름에 대한 수치 연구에 대한 강조와 같은 LOC 시스템 하에서 혈류 역학의 수치 모델링의 문제를 식별합니다.
전기역학 현상을 연구하는 동안 제타 전위 조건에 대한 보다 실용적인 가정도 강조됩니다. 본 연구는 모세관 및 전기삼투압에 의해 구동되는 미세유체 시스템의 혈류 역학에 대한 포괄적이고 학제적인 관점을 제공하는 것을 목표로 한다.
1.1. Microfluidic Flow in Lab-on-a-Chip (LOC) Systems
Over the past several decades, the ability to control and utilize fluid flow patterns at microscales has gained considerable interest across a myriad of scientific and engineering disciplines, leading to growing interest in scientific research of microfluidics.
(1) Microfluidics, an interdisciplinary field that straddles physics, engineering, and biotechnology, is dedicated to the behavior, precise control, and manipulation of fluids geometrically constrained to a small, typically submillimeter, scale.
(2) The engineering community has increasingly focused on microfluidics, exploring different driving forces to enhance working fluid transport, with the aim of accurately and efficiently describing, controlling, designing, and applying microfluidic flow principles and transport phenomena, particularly for miniaturized applications.
(3) This attention has chiefly been fueled by the potential to revolutionize diagnostic and therapeutic techniques in the biomedical and pharmaceutical sectorsUnder various driving forces in microfluidic flows, intriguing transport phenomena have bolstered confidence in sustainable and efficient applications in fields such as pharmaceutical, biochemical, and environmental science. The “lab-on-a-chip” (LOC) system harnesses microfluidic flow to enable fluid processing and the execution of laboratory tasks on a chip-sized scale. LOC systems have played a vital role in the miniaturization of laboratory operations such as mixing, chemical reaction, separation, flow control, and detection on small devices, where a wide variety of fluids is adapted. Biological fluid flow like blood and other viscoelastic fluids are notably studied among the many working fluids commonly utilized by LOC systems, owing to the optimization in small fluid sample volumed, rapid response times, precise control, and easy manipulation of flow patterns offered by the system under various driving forces.
(4)The driving forces in blood flow can be categorized as passive or active transport mechanisms and, in some cases, both. Under various transport mechanisms, the unique design of microchannels enables different functionalities in driving, mixing, separating, and diagnosing blood and drug delivery in the blood.
(5) Understanding and manipulating these driving forces are crucial for optimizing the performance of a LOC system. Such knowledge presents the opportunity to achieve higher efficiency and reliability in addressing cellular level challenges in medical diagnostics, forensic studies, cancer detection, and other fundamental research areas, for applications of point-of-care (POC) devices.
1.2. Engineering Approach of Microfluidic Transport Phenomena in LOC Systems
Different transport mechanisms exhibit unique properties at submillimeter length scales in microfluidic devices, leading to significant transport phenomena that differ from those of macroscale flows. An in-depth understanding of these unique transport phenomena under microfluidic systems is often required in fluidic mechanics to fully harness the potential functionality of a LOC system to obtain systematically designed and precisely controlled transport of microfluids under their respective driving force. Fluid mechanics is considered a vital component in chemical engineering, enabling the analysis of fluid behaviors in various unit designs, ranging from large-scale reactors to separation units. Transport phenomena in fluid mechanics provide a conceptual framework for analytically and descriptively explaining why and how experimental results and physiological phenomena occur. The Navier–Stokes (N–S) equation, along with other governing equations, is often adapted to accurately describe fluid dynamics by accounting for pressure, surface properties, velocity, and temperature variations over space and time. In addition, limiting factors and nonidealities for these governing equations should be considered to impose corrections for empirical consistency before physical models are assembled for more accurate controls and efficiency. Microfluidic flow systems often deviate from ideal conditions, requiring adjustments to the standard governing equations. These deviations could arise from factors such as viscous effects, surface interactions, and non-Newtonian fluid properties from different microfluid types and geometrical layouts of microchannels. Addressing these nonidealities supports the refining of theoretical models and prediction accuracy for microfluidic flow behaviors.
The analytical calculation of coupled nonlinear governing equations, which describes the material and energy balances of systems under ideal conditions, often requires considerable computational efforts. However, advancements in computation capabilities, cost reduction, and improved accuracy have made numerical simulations using different numerical and modeling methods a powerful tool for effectively solving these complex coupled equations and modeling various transport phenomena. Computational fluid dynamics (CFD) is a numerical technique used to investigate the spatial and temporal distribution of various flow parameters. It serves as a critical approach to provide insights and reasoning for decision-making regarding the optimal designs involving fluid dynamics, even prior to complex physical model prototyping and experimental procedures. The integration of experimental data, theoretical analysis, and reliable numerical simulations from CFD enables systematic variation of analytical parameters through quantitative analysis, where adjustment to delivery of blood flow and other working fluids in LOC systems can be achieved.
Numerical methods such as the Finite-Difference Method (FDM), Finite-Element-Method (FEM), and Finite-Volume Method (FVM) are heavily employed in CFD and offer diverse approaches to achieve discretization of Eulerian flow equations through filling a mesh of the flow domain. A more in-depth review of numerical methods in CFD and its application for blood flow simulation is provided in Section 2.2.2.
1.3. Scope of the Review
In this Review, we explore and characterize the blood flow phenomena within the LOC systems, utilizing both physiological and engineering modeling approaches. Similar approaches will be taken to discuss capillary-driven flow and electric-osmotic flow (EOF) under electrokinetic phenomena as a passive and active transport scheme, respectively, for blood transport in LOC systems. Such an analysis aims to bridge the gap between physical (experimental) and engineering (analytical) perspectives in studying and manipulating blood flow delivery by different driving forces in LOC systems. Moreover, the Review hopes to benefit the interests of not only blood flow control in LOC devices but also the transport of viscoelastic fluids, which are less studied in the literature compared to that of Newtonian fluids, in LOC systems.
Section 2 examines the complex interplay between viscoelastic properties of blood and blood flow patterns under shear flow in LOC systems, while engineering numerical modeling approaches for blood flow are presented for assistance. Sections 3 and 4 look into the theoretical principles, numerical governing equations, and modeling methodologies for capillary driven flow and EOF in LOC systems as well as their impact on blood flow dynamics through the quantification of key parameters of the two driving forces. Section 5 concludes the characterized blood flow transport processes in LOC systems under these two forces. Additionally, prospective areas of research in improving the functionality of LOC devices employing blood and other viscoelastic fluids and potentially justifying mechanisms underlying microfluidic flow patterns outside of LOC systems are presented. Finally, the challenges encountered in the numerical studies of blood flow under LOC systems are acknowledged, paving the way for further research.
Blood, an essential physiological fluid in the human body, serves the vital role of transporting oxygen and nutrients throughout the body. Additionally, blood is responsible for suspending various blood cells including erythrocytes (red blood cells or RBCs), leukocytes (white blood cells), and thrombocytes (blood platelets) in a plasma medium.Among the cells mentioned above, red blood cells (RBCs) comprise approximately 40–45% of the volume of healthy blood.
(7) An RBC possesses an inherent elastic property with a biconcave shape of an average diameter of 8 μm and a thickness of 2 μm. This biconcave shape maximizes the surface-to-volume ratio, allowing RBCs to endure significant distortion while maintaining their functionality.
(8,9) Additionally, the biconcave shape optimizes gas exchange, facilitating efficient uptake of oxygen due to the increased surface area. The inherent elasticity of RBCs allows them to undergo substantial distortion from their original biconcave shape and exhibits high flexibility, particularly in narrow channels.RBC deformability enables the cell to deform from a biconcave shape to a parachute-like configuration, despite minor differences in RBC shape dynamics under shear flow between initial cell locations. As shown in Figure 1(a), RBCs initiating with different resting shapes and orientations displaying display a similar deformation pattern
(10) in terms of its shape. Shear flow induces an inward bending of the cell at the rear position of the rim to the final bending position,
(11) resulting in an alignment toward the same position of the flow direction.
Figure 1. Images of varying deformation of RBCs and different dynamic blood flow behaviors. (a) The deforming shape behavior of RBCs at four different initiating positions under the same experimental conditions of a flow from left to right, (10) (b) RBC aggregation, (13) (c) CFL region. (18) Reproduced with permission from ref (10). Copyright 2011 Elsevier. Reproduced with permission from ref (13). Copyright 2022 The Authors, under the terms of the Creative Commons (CC BY 4.0) License https://creativecommons.org/licenses/by/4.0/. Reproduced with permission from ref (18). Copyright 2019 Elsevier.
The flexible property of RBCs enables them to navigate through narrow capillaries and traverse a complex network of blood vessels. The deformability of RBCs depends on various factors, including the channel geometry, RBC concentration, and the elastic properties of the RBC membrane.
(12) Both flexibility and deformability are vital in the process of oxygen exchange among blood and tissues throughout the body, allowing cells to flow in vessels even smaller than the original cell size prior to deforming.As RBCs serve as major components in blood, their collective dynamics also hugely affect blood rheology. RBCs exhibit an aggregation phenomenon due to cell to cell interactions, such as adhesion forces, among populated cells, inducing unique blood flow patterns and rheological behaviors in microfluidic systems. For blood flow in large vessels between a diameter of 1 and 3 cm, where shear rates are not high, a constant viscosity and Newtonian behavior for blood can be assumed. However, under low shear rate conditions (0.1 s
–1) in smaller vessels such as the arteries and venules, which are within a diameter of 0.2 mm to 1 cm, blood exhibits non-Newtonian properties, such as shear-thinning viscosity and viscoelasticity due to RBC aggregation and deformability. The nonlinear viscoelastic property of blood gives rise to a complex relationship between viscosity and shear rate, primarily influenced by the highly elastic behavior of RBCs. A wide range of research on the transient behavior of the RBC shape and aggregation characteristics under varied flow circumstances has been conducted, aiming to obtain a better understanding of the interaction between blood flow shear forces from confined flows.
For a better understanding of the unique blood flow structures and rheological behaviors in microfluidic systems, some blood flow patterns are introduced in the following section.
2.1.1. RBC Aggregation
RBC aggregation is a vital phenomenon to be considered when designing LOC devices due to its impact on the viscosity of the bulk flow. Under conditions of low shear rate, such as in stagnant or low flow rate regions, RBCs tend to aggregate, forming structures known as rouleaux, resembling stacks of coins as shown in Figure 1(b).
(13) The aggregation of RBCs increases the viscosity at the aggregated region,
(14) hence slowing down the overall blood flow. However, when exposed to high shear rates, RBC aggregates disaggregate. As shear rates continue to increase, RBCs tend to deform, elongating and aligning themselves with the direction of the flow.
(15) Such a dynamic shift in behavior from the cells in response to the shear rate forms the basis of the viscoelastic properties observed in whole blood. In essence, the viscosity of the blood varies according to the shear rate conditions, which are related to the velocity gradient of the system. It is significant to take the intricate relationship between shear rate conditions and the change of blood viscosity due to RBC aggregation into account since various flow driving conditions may induce varied effects on the degree of aggregation.
2.1.2. Fåhræus-Lindqvist Effect
The Fåhræus–Lindqvist (FL) effect describes the gradual decrease in the apparent viscosity of blood as the channel diameter decreases.
(16) This effect is attributed to the migration of RBCs toward the central region in the microchannel, where the flow rate is higher, due to the presence of higher pressure and asymmetric distribution of shear forces. This migration of RBCs, typically observed at blood vessels less than 0.3 mm, toward the higher flow rate region contributes to the change in blood viscosity, which becomes dependent on the channel size. Simultaneously, the increase of the RBC concentration in the central region of the microchannel results in the formation of a less viscous region close to the microchannel wall. This region called the Cell-Free Layer (CFL), is primarily composed of plasma.
(17) The combination of the FL effect and the following CFL formation provides a unique phenomenon that is often utilized in passive and active plasma separation mechanisms, involving branched and constriction channels for various applications in plasma separation using microfluidic systems.
2.1.3. Cell-Free Layer Formation
In microfluidic blood flow, RBCs form aggregates at the microchannel core and result in a region that is mostly devoid of RBCs near the microchannel walls, as shown in Figure 1(c).
(18) The region is known as the cell-free layer (CFL). The CFL region is often known to possess a lower viscosity compared to other regions within the blood flow due to the lower viscosity value of plasma when compared to that of the aggregated RBCs. Therefore, a thicker CFL region composed of plasma correlates to a reduced apparent whole blood viscosity.
(19) A thicker CFL region is often established following the RBC aggregation at the microchannel core under conditions of decreasing the tube diameter. Apart from the dependence on the RBC concentration in the microchannel core, the CFL thickness is also affected by the volume concentration of RBCs, or hematocrit, in whole blood, as well as the deformability of RBCs. Given the influence CFL thickness has on blood flow rheological parameters such as blood flow rate, which is strongly dependent on whole blood viscosity, investigating CFL thickness under shear flow is crucial for LOC systems accounting for blood flow.
2.1.4. Plasma Skimming in Bifurcation Networks
The uneven arrangement of RBCs in bifurcating microchannels, commonly termed skimming bifurcation, arises from the axial migration of RBCs within flowing streams. This uneven distribution contributes to variations in viscosity across differing sizes of bifurcating channels but offers a stabilizing effect. Notably, higher flow rates in microchannels are associated with increased hematocrit levels, resulting in higher viscosity compared with those with lower flow rates. Parametric investigations on bifurcation angle,
(21) and RBC dynamics, including aggregation and deformation,
(22) may alter the varying viscosity of blood and its flow behavior within microchannels.
2.2. Modeling on Blood Flow Dynamics
2.2.1. Blood Properties and Mathematical Models of Blood Rheology
Under different shear rate conditions in blood flow, the elastic characteristics and dynamic changes of the RBC induce a complex velocity and stress relationship, resulting in the incompatibility of blood flow characterization through standard presumptions of constant viscosity used for Newtonian fluid flow. Blood flow is categorized as a viscoelastic non-Newtonian fluid flow where constitutive equations governing this type of flow take into consideration the nonlinear viscometric properties of blood. To mathematically characterize the evolving blood viscosity and the relationship between the elasticity of RBC and the shear blood flow, respectively, across space and time of the system, a stress tensor (τ) defined by constitutive models is often coupled in the Navier–Stokes equation to account for the collective impact of the constant dynamic viscosity (η) and the elasticity from RBCs on blood flow.The dynamic viscosity of blood is heavily dependent on the shear stress applied to the cell and various parameters from the blood such as hematocrit value, plasma viscosity, mechanical properties of the RBC membrane, and red blood cell aggregation rate. The apparent blood viscosity is considered convenient for the characterization of the relationship between the evolving blood viscosity and shear rate, which can be defined by Casson’s law, as shown in eq 1.
𝜇=𝜏0𝛾˙+2𝜂𝜏0𝛾˙⎯⎯⎯⎯⎯⎯⎯√+𝜂�=�0�˙+2��0�˙+�
(1)where τ
0 is the yield stress–stress required to initiate blood flow motion, η is the Casson rheological constant, and γ̇ is the shear rate. The value of Casson’s law parameters under blood with normal hematocrit level can be defined as τ
0 = 0.0056 Pa and η = 0.0035 Pa·s.
(23) With the known property of blood and Casson’s law parameters, an approximation can be made to the dynamic viscosity under various flow condition domains. The Power Law model is often employed to characterize the dynamic viscosity in relation to the shear rate, since precise solutions exist for specific geometries and flow circumstances, acting as a fundamental standard for definition. The Carreau and Carreau–Yasuda models can be advantageous over the Power Law model due to their ability to evaluate the dynamic viscosity at low to zero shear rate conditions. However, none of the above-mentioned models consider the memory or other elastic behavior of blood and its RBCs. Some other commonly used mathematical models and their constants for the non-Newtonian viscosity property characterization of blood are listed in Table 1 below.
(24−26)Table 1. Comparison of Various Non-Newtonian Models for Blood Viscosity
The blood rheology is commonly known to be influenced by two key physiological factors, namely, the hematocrit value (H
t) and the fibrinogen concentration (c
f), with an average value of 42% and 0.252 gd·L
–1, respectively. Particularly in low shear conditions, the presence of varying fibrinogen concentrations affects the tendency for aggregation and rouleaux formation, while the occurrence of aggregation is contingent upon specific levels of hematocrit.
(28) modifies the Casson model through emphasizing its reliance on hematocrit and fibrinogen concentration parameter values, owing to the extensive knowledge of the two physiological blood parameters.The viscoelastic response of blood is heavily dependent on the elasticity of the RBC, which is defined by the relationship between the deformation and stress relaxation from RBCs under a specific location of shear flow as a function of the velocity field. The stress tensor is usually characterized by constitutive equations such as the Upper-Convected Maxwell Model
(30) to track the molecule effects under shear from different driving forces. The prominent non-Newtonian features, such as shear thinning and yield stress, have played a vital role in the characterization of blood rheology, particularly with respect to the evaluation of yield stress under low shear conditions. The nature of stress measurement in blood, typically on the order of 1 mPa, is challenging due to its low magnitude. The occurrence of the CFL complicates the measurement further due to the significant decrease in apparent viscosity near the wall over time and a consequential disparity in viscosity compared to the bulk region.In addition to shear thinning viscosity and yield stress, the formation of aggregation (rouleaux) from RBCs under low shear rates also contributes to the viscoelasticity under transient flow
(32) of whole blood. Given the difficulty in evaluating viscoelastic behavior of blood under low strain magnitudes and limitations in generalized Newtonian models, the utilization of viscoelastic models is advocated to encompass elasticity and delineate non-shear components within the stress tensor. Extending from the Oldroyd-B model, Anand et al.
(33) developed a viscoelastic model framework for adapting elasticity within blood samples and predicting non-shear stress components. However, to also address the thixotropic effects, the model developed by Horner et al.
(34) serves as a more comprehensive approach than the viscoelastic model from Anand et al. Thixotropy
(32) typically occurs from the structural change of the rouleaux, where low shear rate conditions induce rouleaux formation. Correspondingly, elasticity increases, while elasticity is more representative of the isolated RBCs, under high shear rate conditions. The model of Horner et al.
(34) considers the contribution of rouleaux to shear stress, taking into account factors such as the characteristic time for Brownian aggregation, shear-induced aggregation, and shear-induced breakage. Subsequent advancements in the model from Horner et al. often revolve around refining the three aforementioned key terms for a more substantial characterization of rouleaux dynamics. Notably, this has led to the recently developed mHAWB model
(35) and other model iterations to enhance the accuracy of elastic and viscoelastic contributions to blood rheology, including the recently improved model suggested by Armstrong et al.
Numerical simulation has become increasingly more significant in analyzing the geometry, boundary layers of flow, and nonlinearity of hyperbolic viscoelastic flow constitutive equations. CFD is a powerful and efficient tool utilizing numerical methods to solve the governing hydrodynamic equations, such as the Navier–Stokes (N–S) equation, continuity equation, and energy conservation equation, for qualitative evaluation of fluid motion dynamics under different parameters. CFD overcomes the challenge of analytically solving nonlinear forms of differential equations by employing numerical methods such as the Finite-Difference Method (FDM), Finite-Element Method (FEM), and Finite-Volume Method (FVM) to discretize and solve the partial differential equations (PDEs), allowing for qualitative reproduction of transport phenomena and experimental observations. Different numerical methods are chosen to cope with various transport systems for optimization of the accuracy of the result and control of error during the discretization process.FDM is a straightforward approach to discretizing PDEs, replacing the continuum representation of equations with a set of finite-difference equations, which is typically applied to structured grids for efficient implementation in CFD programs.
(37) However, FDM is often limited to simple geometries such as rectangular or block-shaped geometries and struggles with curved boundaries. In contrast, FEM divides the fluid domain into small finite grids or elements, approximating PDEs through a local description of physics.
(38) All elements contribute to a large, sparse matrix solver. However, FEM may not always provide accurate results for systems involving significant deformation and aggregation of particles like RBCs due to large distortion of grids.
(39) FVM evaluates PDEs following the conservation laws and discretizes the selected flow domain into small but finite size control volumes, with each grid at the center of a finite volume.
(40) The divergence theorem allows the conversion of volume integrals of PDEs with divergence terms into surface integrals of surface fluxes across cell boundaries. Due to its conservation property, FVM offers efficient outcomes when dealing with PDEs that embody mass, momentum, and energy conservation principles. Furthermore, widely accessible software packages like the OpenFOAM toolbox
(41) include a viscoelastic solver, making it an attractive option for viscoelastic fluid flow modeling.
The complexity in the blood flow simulation arises from deformability and aggregation that RBCs exhibit during their interaction with neighboring cells under different shear rate conditions induced by blood flow. Numerical models coupled with simulation programs have been applied as a groundbreaking method to predict such unique rheological behavior exhibited by RBCs and whole blood. The conventional approach of a single-phase flow simulation is often applied to blood flow simulations within large vessels possessing a moderate shear rate. However, such a method assumes the properties of plasma, RBCs and other cellular components to be evenly distributed as average density and viscosity in blood, resulting in the inability to simulate the mechanical dynamics, such as RBC aggregation under high-shear flow field, inherent in RBCs. To accurately describe the asymmetric distribution of RBC and blood flow, multiphase flow simulation, where numerical simulations of blood flows are often modeled as two immiscible phases, RBCs and blood plasma, is proposed. A common assumption is that RBCs exhibit non-Newtonian behavior while the plasma is treated as a continuous Newtonian phase.Numerous multiphase numerical models have been proposed to simulate the influence of RBCs on blood flow dynamics by different assumptions. In large-scale simulations (above the millimeter range), continuum-based methods are wildly used due to their lower computational demands.
(43) Eulerian multiphase flow simulations offer the solution of a set of conservation equations for each separate phase and couple the phases through common pressure and interphase exchange coefficients. Xu et al.
(44) utilized the combined finite-discrete element method (FDEM) to replicate the dynamic behavior and distortion of RBCs subjected to fluidic forces, utilizing the Johnson–Kendall–Roberts model
(45) to define the adhesive forces of cell-to-cell interactions. The iterative direct-forcing immersed boundary method (IBM) is commonly employed in simulations of the fluid–cell interface of blood. This method effectively captures the intricacies of the thin and flexible RBC membranes within various external flow fields.
(44) also adopts this approach to bridge the fluid dynamics and RBC deformation through IBM. Yoon and You utilized the Maxwell model to define the viscosity of the RBC membrane.
(47) It was discovered that the Maxwell model could represent the stress relaxation and unloading processes of the cell. Furthermore, the reduced flexibility of an RBC under particular situations such as infection is specified, which was unattainable by the Kelvin–Voigt model
(48) when compared to the Maxwell model in the literature. The Yeoh hyperplastic material model was also adapted to predict the nonlinear elasticity property of RBCs with FEM employed to discretize the RBC membrane using shell-type elements. Gracka et al.
(49) developed a numerical CFD model with a finite-volume parallel solver for multiphase blood flow simulation, where an updated Maxwell viscoelasticity model and a Discrete Phase Model are adopted. In the study, the adapted IBM, based on unstructured grids, simulates the flow behavior and shape change of the RBCs through fluid-structure coupling. It was found that the hybrid Euler–Lagrange (E–L) approach
(50) for the development of the multiphase model offered better results in the simulated CFL region in the microchannels.To study the dynamics of individual behaviors of RBCs and the consequent non-Newtonian blood flow, cell-shape-resolved computational models are often adapted. The use of the boundary integral method has become prevalent in minimizing computational expenses, particularly in the exclusive determination of fluid velocity on the surfaces of RBCs, incorporating the option of employing IBM or particle-based techniques. The cell-shaped-resolved method has enabled an examination of cell to cell interactions within complex ambient or pulsatile flow conditions
(51) surrounding RBC membranes. Recently, Rydquist et al.
(52) have looked to integrate statistical information from macroscale simulations to obtain a comprehensive overview of RBC behavior within the immediate proximity of the flow through introduction of respective models characterizing membrane shape definition, tension, bending stresses of RBC membranes.At a macroscopic scale, continuum models have conventionally been adapted for assessing blood flow dynamics through the application of elasticity theory and fluid dynamics. However, particle-based methods are known for their simplicity and adaptability in modeling complex multiscale fluid structures. Meshless methods, such as the boundary element method (BEM), smoothed particle hydrodynamics (SPH), and dissipative particle dynamics (DPD), are often used in particle-based characterization of RBCs and the surrounding fluid. By representing the fluid as discrete particles, meshless methods provide insights into the status and movement of the multiphase fluid. These methods allow for the investigation of cellular structures and microscopic interactions that affect blood rheology. Non-confronting mesh methods like IBM can also be used to couple a fluid solver such as FEM, FVM, or the Lattice Boltzmann Method (LBM) through membrane representation of RBCs. In comparison to conventional CFD methods, LBM has been viewed as a favorable numerical approach for solving the N–S equations and the simulation of multiphase flows. LBM exhibits the notable advantage of being amenable to high-performance parallel computing environments due to its inherently local dynamics. In contrast to DPD and SPH where RBC membranes are modeled as physically interconnected particles, LBM employs the IBM to account for the deformation dynamics of RBCs
(53,54) under shear flows in complex channel geometries.
(54,55) However, it is essential to acknowledge that the utilization of LBM in simulating RBC flows often entails a significant computational overhead, being a primary challenge in this context. Krüger et al.
(56) proposed utilizing LBM as a fluid solver, IBM to couple the fluid and FEM to compute the response of membranes to deformation under immersed fluids. This approach decouples the fluid and membranes but necessitates significant computational effort due to the requirements of both meshes and particles.Despite the accuracy of current blood flow models, simulating complex conditions remains challenging because of the high computational load and cost. Balachandran Nair et al.
(57) suggested a reduced order model of RBC under the framework of DEM, where the RBC is represented by overlapping constituent rigid spheres. The Morse potential force is adapted to account for the RBC aggregation exhibited by cell to cell interactions among RBCs at different distances. Based upon the IBM, the reduced-order RBC model is adapted to simulate blood flow transport for validation under both single and multiple RBCs with a resolved CFD-DEM solver.
(58) In the resolved CFD-DEM model, particle sizes are larger than the grid size for a more accurate computation of the surrounding flow field. A continuous forcing approach is taken to describe the momentum source of the governing equation prior to discretization, which is different from a Direct Forcing Method (DFM).
(59) As no body-conforming moving mesh is required, the continuous forcing approach offers lower complexity and reduced cost when compared to the DFM. Piquet et al.
(60) highlighted the high complexity of the DFM due to its reliance on calculating an additional immersed boundary flux for the velocity field to ensure its divergence-free condition.The fluid–structure interaction (FSI) method has been advocated to connect the dynamic interplay of RBC membranes and fluid plasma within blood flow such as the coupling of continuum–particle interactions. However, such methodology is generally adapted for anatomical configurations such as arteries
(63) where both the structural components and the fluid domain undergo substantial deformation due to the moving boundaries. Due to the scope of the Review being blood flow simulation within microchannels of LOC devices without deformable boundaries, the Review of the FSI method will not be further carried out.In general, three numerical methods are broadly used: mesh-based, particle-based, and hybrid mesh–particle techniques, based on the spatial scale and the fundamental numerical approach, mesh-based methods tend to neglect the effects of individual particles, assuming a continuum and being efficient in terms of time and cost. However, the particle-based approach highlights more of the microscopic and mesoscopic level, where the influence of individual RBCs is considered. A review from Freund et al.
(64) addressed the three numerical methodologies and their respective modeling approaches of RBC dynamics. Given the complex mechanics and the diverse levels of study concerning numerical simulations of blood and cellular flow, a broad spectrum of numerical methods for blood has been subjected to extensive review.
(65) offered an extensive review of the application of the DPD, SPH, and LBM for numerical simulations of RBC, while Rathnayaka et al.
(67) conducted a review of the particle-based numerical modeling for liquid marbles through drawing parallels to the transport of RBCs in microchannels. A comparative analysis between conventional CFD methods and particle-based approaches for cellular and blood flow dynamic simulation can be found under the review by Arabghahestani et al.
(69) offer an overview of both continuum-based models at micro/macroscales and multiscale particle-based models encompassing various length and temporal dimensions. Furthermore, these reviews deliberate upon the potential of coupling continuum-particle methods for blood plasma and RBC modeling. Arciero et al.
(70) investigated various modeling approaches encompassing cellular interactions, such as cell to cell or plasma interactions and the individual cellular phases. A concise overview of the reviews is provided in Table 2 for reference.
Table 2. List of Reviews for Numerical Approaches Employed in Blood Flow Simulation
Capillary driven (CD) flow is a pivotal mechanism in passive microfluidic flow systems
(9) such as the blood circulation system and LOC systems.
(71) CD flow is essentially the movement of a liquid to flow against drag forces, where the capillary effect exerts a force on the liquid at the borders, causing a liquid–air meniscus to flow despite gravity or other drag forces. A capillary pressure drops across the liquid–air interface with surface tension in the capillary radius and contact angle. The capillary effect depends heavily on the interaction between the different properties of surface materials. Different values of contact angles can be manipulated and obtained under varying levels of surface wettability treatments to manipulate the surface properties, resulting in different CD blood delivery rates for medical diagnostic device microchannels. CD flow techniques are appealing for many LOC devices, because they require no external energy. However, due to the passive property of liquid propulsion by capillary forces and the long-term instability of surface treatments on channel walls, the adaptability of CD flow in geometrically complex LOC devices may be limited.
3.2. Theoretical and Numerical Modeling of Capillary Driven Blood Flow
3.2.1. Theoretical Basis and Assumptions of Microfluidic Flow
The study of transport phenomena regarding either blood flow driven by capillary forces or externally applied forces under microfluid systems all demands a comprehensive recognition of the significant differences in flow dynamics between microscale and macroscale. The fundamental assumptions and principles behind fluid transport at the microscale are discussed in this section. Such a comprehension will lay the groundwork for the following analysis of the theoretical basis of capillary forces and their role in blood transport in LOC systems.
At the macroscale, fluid dynamics are often strongly influenced by gravity due to considerable fluid mass. However, the high surface to volume ratio at the microscale shifts the balance toward surface forces (e.g., surface tension and viscous forces), much larger than the inertial force. This difference gives rise to transport phenomena unique to microscale fluid transport, such as the prevalence of laminar flow due to a very low Reynolds number (generally lower than 1). Moreover, the fluid in a microfluidic system is often assumed to be incompressible due to the small flow velocity, indicating constant fluid density in both space and time.Microfluidic flow behaviors are governed by the fundamental principles of mass and momentum conservation, which are encapsulated in the continuity equation and the Navier–Stokes (N–S) equation. The continuity equation describes the conservation of mass, while the N–S equation captures the spatial and temporal variations in velocity, pressure, and other physical parameters. Under the assumption of the negligible influence of gravity in microfluidic systems, the continuity equation and the Eulerian representation of the incompressible N–S equation can be expressed as follows:
∇·𝐮⇀=0∇·�⇀=0
(7)
−∇𝑝+𝜇∇2𝐮⇀+∇·𝝉⇀−𝐅⇀=0−∇�+�∇2�⇀+∇·�⇀−�⇀=0
(8)Here, p is the pressure, u is the fluid viscosity,
𝝉⇀�⇀ represents the stress tensor, and F is the body force exerted by external forces if present.
3.2.2. Theoretical Basis and Modeling of Capillary Force in LOC Systems
The capillary force is often the major driving force to manipulate and transport blood without an externally applied force in LOC systems. Forces induced by the capillary effect impact the free surface of fluids and are represented not directly in the Navier–Stokes equations but through the pressure boundary conditions of the pressure term p. For hydrophilic surfaces, the liquid generally induces a contact angle between 0° and 30°, encouraging the spread and attraction of fluid under a positive cos θ condition. For this condition, the pressure drop becomes positive and generates a spontaneous flow forward. A hydrophobic solid surface repels the fluid, inducing minimal contact. Generally, hydrophobic solids exhibit a contact angle larger than 90°, inducing a negative value of cos θ. Such a value will result in a negative pressure drop and a flow in the opposite direction. The induced contact angle is often utilized to measure the wall exposure of various surface treatments on channel walls where different wettability gradients and surface tension effects for CD flows are established. Contact angles between different interfaces are obtainable through standard values or experimental methods for reference.
(72)For the characterization of the induced force by the capillary effect, the Young–Laplace (Y–L) equation
(73) is widely employed. In the equation, the capillary is considered a pressure boundary condition between the two interphases. Through the Y–L equation, the capillary pressure force can be determined, and subsequently, the continuity and momentum balance equations can be solved to obtain the blood filling rate. Kim et al.
(74) studied the effects of concentration and exposure time of a nonionic surfactant, Silwet L-77, on the performance of a polydimethylsiloxane (PDMS) microchannel in terms of plasma and blood self-separation. The study characterized the capillary pressure force by incorporating the Y–L equation and further evaluated the effects of the changing contact angle due to different levels of applied channel wall surface treatments. The expression of the Y–L equation utilized by Kim et al.
(9)where σ is the surface tension of the liquid and θ
b, θ
t, θ
l, and θ
r are the contact angle values between the liquid and the bottom, top, left, and right walls, respectively. A numerical simulation through Coventor software is performed to evaluate the dynamic changes in the filling rate within the microchannel. The simulation results for the blood filling rate in the microchannel are expressed at a specific time stamp, shown in Figure 2. The results portray an increasing instantaneous filling rate of blood in the microchannel following the decrease in contact angle induced by a higher concentration of the nonionic surfactant treated to the microchannel wall.
Figure 2. Numerical simulation of filling rate of capillary driven blood flow under various contact angle conditions at a specific timestamp. (74) Reproduced with permission from ref (74). Copyright 2010 Elsevier.
When in contact with hydrophilic or hydrophobic surfaces, blood forms a meniscus with a contact angle due to surface tension. The Lucas–Washburn (L–W) equation
(75) is one of the pioneering theoretical definitions for the position of the meniscus over time. In addition, the L–W equation provides the possibility for research to obtain the velocity of the blood formed meniscus through the derivation of the meniscus position. The L–W equation
(10)Here L(t) represents the distance of the liquid driven by the capillary forces. However, the generalized L–W equation solely assumes the constant physical properties from a Newtonian fluid rather than considering the non-Newtonian fluid behavior of blood. Cito et al.
(76) constructed an enhanced version of the L–W equation incorporating the power law to consider the RBC aggregation and the FL effect. The non-Newtonian fluid apparent viscosity under the Power Law model is defined as
𝜇=𝑘·(𝛾˙)𝑛−1�=�·(�˙)�−1
(11)where γ̇ is the strain rate tensor defined as
𝛾˙=12𝛾˙𝑖𝑗𝛾˙𝑗𝑖⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯√�˙=12�˙���˙��. The stress tensor term τ is computed as τ = μγ̇
(12)where k is the flow consistency index and n is the power law index, respectively. The power law index, from the Power Law model, characterizes the extent of the non-Newtonian behavior of blood. Both the consistency and power law index rely on blood properties such as hematocrit, the appearance of the FL effect, the formation of RBC aggregates, etc. The updated L–W equation computes the location and velocity of blood flow caused by capillary forces at specified time points within the LOC devices, taking into account the effects of blood flow characteristics such as RBC aggregation and the FL effect on dynamic blood viscosity.Apart from the blood flow behaviors triggered by inherent blood properties, unique flow conditions driven by capillary forces that are portrayed under different microchannel geometries also hold crucial implications for CD blood delivery. Berthier et al.
(77) studied the spontaneous Concus–Finn condition, the condition to initiate the spontaneous capillary flow within a V-groove microchannel, as shown in Figure 3(a) both experimentally and numerically. Through experimental studies, the spontaneous Concus–Finn filament development of capillary driven blood flow is observed, as shown in Figure 3(b), while the dynamic development of blood flow is numerically simulated through CFD simulation.
Figure 3. (a) Sketch of the cross-section of Berthier’s V-groove microchannel, (b) experimental view of blood in the V-groove microchannel, (78) (c) illustration of the dynamic change of the extension of filament from FLOW 3D under capillary flow at three increasing time intervals. (78) Reproduced with permission from ref (78). Copyright 2014 Elsevier.
Berthier et al.
(77) characterized the contact angle needed for the initiation of the capillary driving force at a zero-inlet pressure, through the half-angle (α) of the V-groove geometry layout, and its relation to the Concus–Finn filament as shown below:
(13)Three possible regimes were concluded based on the contact angle value for the initiation of flow and development of Concus–Finn filament:
𝜃>𝜃1𝜃1>𝜃>𝜃0𝜃0no SCFSCF without a Concus−Finn filamentSCF without a Concus−Finn filament{�>�1no SCF�1>�>�0SCF without a Concus−Finn filament�0SCF without a Concus−Finn filament
(14)Under Newton’s Law, the force balance with low Reynolds and Capillary numbers results in the neglect of inertial terms. The force balance between the capillary forces and the viscous force induced by the channel wall is proposed to derive the analytical fluid velocity. This relation between the two forces offers insights into the average flow velocity and the penetration distance function dependent on time. The apparent blood viscosity is defined by Berthier et al.
(23) given in eq 1. The research used the FLOW-3D program from Flow Science Inc. software, which solves transient, free-surface problems using the FDM in multiple dimensions. The Volume of Fluid (VOF) method
(79) is utilized to locate and track the dynamic extension of filament throughout the advancing interface within the channel ahead of the main flow at three progressing time stamps, as depicted in Figure 3(c).
The utilization of external forces, such as electric fields, has significantly broadened the possibility of manipulating microfluidic flow in LOC systems.
(80) Externally applied electric field forces induce a fluid flow from the movement of ions in fluid terms as the “electro-osmotic flow” (EOF).Unique transport phenomena, such as enhanced flow velocity and flow instability, induced by non-Newtonian fluids, particularly viscoelastic fluids, under EOF, have sparked considerable interest in microfluidic devices with simple or complicated geometries within channels.
(81) However, compared to the study of Newtonian fluids and even other electro-osmotic viscoelastic fluid flows, the literature focusing on the theoretical and numerical modeling of electro-osmotic blood flow is limited due to the complexity of blood properties. Consequently, to obtain a more comprehensive understanding of the complex blood flow behavior under EOF, theoretical and numerical studies of the transport phenomena in the EOF section will be based on the studies of different viscoelastic fluids under EOF rather than that of blood specifically. Despite this limitation, we believe these studies offer valuable insights that can help understand the complex behavior of blood flow under EOF.
4.1. EOF Phenomena
Electro-osmotic flow occurs at the interface between the microchannel wall and bulk phase solution. When in contact with the bulk phase, solution ions are absorbed or dissociated at the solid–liquid interface, resulting in the formation of a charge layer, as shown in Figure 4. This charged channel surface wall interacts with both negative and positive ions in the bulk sample, causing repulsion and attraction forces to create a thin layer of immobilized counterions, known as the Stern layer. The induced electric potential from the wall gradually decreases with an increase in the distance from the wall. The Stern layer potential, commonly termed the zeta potential, controls the intensity of the electrostatic interactions between mobile counterions and, consequently, the drag force from the applied electric field. Next to the Stern layer is the diffuse mobile layer, mainly composed of a mobile counterion. These two layers constitute the “electrical double layer” (EDL), the thickness of which is directly proportional to the ionic strength (concentration) of the bulk fluid. The relationship between the two parameters is characterized by a Debye length (λ
D), expressed as
𝜆𝐷=𝜖𝑘B𝑇2(𝑍𝑒)2𝑐0⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯√��=��B�2(��)2�0
(15)where ϵ is the permittivity of the electrolyte solution, k
B is the Boltzmann constant, T is the electron temperature, Z is the integer valence number, e is the elementary charge, and c
0 is the ionic density.
Figure 4. Schematic diagram of an electro-osmotic flow in a microchannel with negative surface charge. (82) Reproduced with permission from ref (82). Copyright 2012 Woodhead Publishing.
When an electric field is applied perpendicular to the EDL, viscous drag is generated due to the movement of excess ions in the EDL. Electro-osmotic forces can be attributed to the externally applied electric potential (ϕ) and the zeta potential, the system wall induced potential by charged walls (ψ). As illustrated in Figure 4, the majority of ions in the bulk phase have a uniform velocity profile, except for a shear rate condition confined within an extremely thin Stern layer. Therefore, EOF displays a unique characteristic of a “near flat” or plug flow velocity profile, different from the parabolic flow typically induced by pressure-driven microfluidic flow (Hagen–Poiseuille flow). The plug-shaped velocity profile of the EOF possesses a high shear rate above the Stern layer.Overall, the EOF velocity magnitude is typically proportional to the Debye Length (λ
D), zeta potential, and magnitude of the externally applied electric field, while a more viscous liquid reduces the EOF velocity.
4.2. Modeling on Electro-osmotic Viscoelastic Fluid Flow
4.2.1. Theoretical Basis of EOF Mechanisms
The EOF of an incompressible viscoelastic fluid is commonly governed by the continuity and incompressible N–S equations, as shown in eqs 7 and 8, where the stress tensor and the electrostatic force term are coupled. The electro-osmotic body force term F, representing the body force exerted by the externally applied electric force, is defined as
𝐹⇀=𝑝𝐸𝐸⇀�⇀=���⇀, where ρ
E and
𝐸⇀�⇀ are the net electric charge density and the applied external electric field, respectively.Numerous models are established to theoretically study the externally applied electric potential and the system wall induced potential by charged walls. The following Laplace equation, expressed as eq 16, is generally adapted and solved to calculate the externally applied potential (ϕ).
∇2𝜙=0∇2�=0
(16)Ion diffusion under applied electric fields, together with mass transport resulting from convection and diffusion, transports ionic solutions in bulk flow under electrokinetic processes. The Nernst–Planck equation can describe these transport methods, including convection, diffusion, and electro-diffusion. Therefore, the Nernst–Planck equation is used to determine the distribution of the ions within the electrolyte. The electric potential induced by the charged channel walls follows the Poisson–Nernst–Plank (PNP) equation, which can be written as eq 17.
i are the diffusion coefficient, ionic concentration, and ionic valence of the ionic species I, respectively. However, due to the high nonlinearity and numerical stiffness introduced by different lengths and time scales from the PNP equations, the Poisson–Boltzmann (PB) model is often considered the major simplified method of the PNP equation to characterize the potential distribution of the EDL region in microchannels. In the PB model, it is assumed that the ionic species in the fluid follow the Boltzmann distribution. This model is typically valid for steady-state problems where charge transport can be considered negligible, the EDLs do not overlap with each other, and the intrinsic potentials are low. It provides a simplified representation of the potential distribution in the EDL region. The PB equation governing the EDL electric potential distribution is described as
0 is the ion bulk concentration, z is the ionic valence, and ε
0 is the electric permittivity in the vacuum. Under low electric potential conditions, an even further simplified model to illustrate the EOF phenomena is the Debye–Hückel (DH) model. The DH model is derived by obtaining a charge density term by expanding the exponential term of the Boltzmann equation in a Taylor series.
4.2.2. EOF Modeling for Viscoelastic Fluids
Many studies through numerical modeling were performed to obtain a deeper understanding of the effect exhibited by externally applied electric fields on viscoelastic flow in microchannels under various geometrical designs. Bello et al.
(83) found that methylcellulose solution, a non-Newtonian polymer solution, resulted in stronger electro-osmotic mobility in experiments when compared to the predictions by the Helmholtz–Smoluchowski equation, which is commonly used to define the velocity of EOF of a Newtonian fluid. Being one of the pioneers to identify the discrepancies between the EOF of Newtonian and non-Newtonian fluids, Bello et al. attributed such discrepancies to the presence of a very high shear rate in the EDL, resulting in a change in the orientation of the polymer molecules. Park and Lee
(84) utilized the FVM to solve the PB equation for the characterization of the electric field induced force. In the study, the concept of fractional calculus for the Oldroyd-B model was adapted to illustrate the elastic and memory effects of viscoelastic fluids in a straight microchannel They observed that fluid elasticity and increased ratio of viscoelastic fluid contribution to overall fluid viscosity had a significant impact on the volumetric flow rate and sensitivity of velocity to electric field strength compared to Newtonian fluids. Afonso et al.
(85) derived an analytical expression for EOF of viscoelastic fluid between parallel plates using the DH model to account for a zeta potential condition below 25 mV. The study established the understanding of the electro-osmotic viscoelastic fluid flow under low zeta potential conditions. Apart from the electrokinetic forces, pressure forces can also be coupled with EOF to generate a unique fluid flow behavior within the microchannel. Sousa et al.
(86) analytically studied the flow of a standard viscoelastic solution by combining the pressure gradient force with an externally applied electric force. It was found that, at a near wall skimming layer and the outer layer away from the wall, macromolecules migrating away from surface walls in viscoelastic fluids are observed. In the study, the Phan-Thien Tanner (PTT) constitutive model is utilized to characterize the viscoelastic properties of the solution. The approach is found to be valid when the EDL is much thinner than the skimming layer under an enhanced flow rate. Zhao and Yang
(87) solved the PB equation and Carreau model for the characterization of the EOF mechanism and non-Newtonian fluid respectively through the FEM. The numerical results depict that, different from the EOF of Newtonian fluids, non-Newtonian fluids led to an increase of electro-osmotic mobility for shear thinning fluids but the opposite for shear thickening fluids.Like other fluid transport driving forces, EOF within unique geometrical layouts also portrays unique transport phenomena. Pimenta and Alves
(88) utilized the FVM to perform numerical simulations of the EOF of viscoelastic fluids considering the PB equation and the Oldroyd-B model, in a cross-slot and flow-focusing microdevices. It was found that electroelastic instabilities are formed due to the development of large stresses inside the EDL with streamlined curvature at geometry corners. Bezerra et al.
(89) used the FDM to numerically analyze the vortex formation and flow instability from an electro-osmotic non-Newtonian fluid flow in a microchannel with a nozzle geometry and parallel wall geometry setting. The PNP equation is utilized to characterize the charge motion in the EOF and the PTT model for non-Newtonian flow characterization. A constriction geometry is commonly utilized in blood flow adapted in LOC systems due to the change in blood flow behavior under narrow dimensions in a microchannel. Ji et al.
(90) recently studied the EOF of viscoelastic fluid in a constriction microchannel connected by two relatively big reservoirs on both ends (as seen in Figure 5) filled with the polyacrylamide polymer solution, a viscoelastic fluid, and an incompressible monovalent binary electrolyte solution KCl.
Figure 5. Schematic diagram of a negatively charged constriction microchannel connected to two reservoirs at both ends. An electro-osmotic flow is induced in the system by the induced potential difference between the anode and cathode. (90) Reproduced with permission from ref (90). Copyright 2021 The Authors, under the terms of the Creative Commons (CC BY 4.0) License https://creativecommons.org/licenses/by/4.0/.
In studying the EOF of viscoelastic fluids, the Oldroyd-B model is often utilized to characterize the polymeric stress tensor and the deformation rate of the fluid. The Oldroyd-B model is expressed as follows:
𝜏=𝜂p𝜆(𝐜−𝐈)�=�p�(�−�)
(19)where η
p, λ, c, and I represent the polymer dynamic viscosity, polymer relaxation time, symmetric conformation tensor of the polymer molecules, and the identity matrix, respectively.A log-conformation tensor approach is taken to prevent convergence difficulty induced by the viscoelastic properties. The conformation tensor (c) in the polymeric stress tensor term is redefined by a new tensor (Θ) based on the natural logarithm of the c. The new tensor is defined as
Θ=ln(𝐜)=𝐑ln(𝚲)𝐑Θ=ln(�)=�ln(�)�
(20)in which Λ is the diagonal matrix and R is the orthogonal matrix.Under the new conformation tensor, the induced EOF of a viscoelastic fluid is governed by the continuity and N–S equations adapting the Oldroyd-B model, which is expressed as
(21)where Ω and B represent the anti-symmetric matrix and the symmetric traceless matrix of the decomposition of the velocity gradient tensor ∇u, respectively. The conformation tensor can be recovered by c = exp(Θ). The PB model and Laplace equation are utilized to characterize the charged channel wall induced potential and the externally applied potential.The governing equations are numerically solved through the FVM by RheoTool,
(42) an open-source viscoelastic EOF solver on the OpenFOAM platform. A SIMPLEC (Semi-Implicit Method for Pressure Linked Equations-Consistent) algorithm was applied to solve the velocity-pressure coupling. The pressure field and velocity field were computed by the PCG (Preconditioned Conjugate Gradient) solver and the PBiCG (Preconditioned Biconjugate Gradient) solver, respectively.Ranging magnitudes of an applied electric field or fluid concentration induce both different streamlines and velocity magnitudes at various locations and times of the microchannel. In the study performed by Ji et al.,
(90) notable fluctuation of streamlines and vortex formation is formed at the upper stream entrance of the constriction as shown in Figure 6(a) and (b), respectively, due to the increase of electrokinetic effect, which is seen as a result of the increase in polymeric stress (τ
xx).
(90) The contraction geometry enhances the EOF velocity within the constriction channel under high E
app condition (600 V/cm). Such phenomena can be attributed to the dependence of electro-osmotic viscoelastic fluid flow on the system wall surface and bulk fluid properties.
Figure 6. Schematic diagram of vortex formation and streamlines of EOF depicting flow instability at (a) 1.71 s and (b) 1.75 s. Spatial distribution of the elastic normal stress at (c) high Eapp condition. Streamline of an electro-osmotic flow under Eapp of 600 V/cm (90) for (d) non-Newtonian and (e) Newtonian fluid through a constriction geometry. Reproduced with permission from ref (90). Copyright 2021 The Authors, under the terms of the Creative Commons (CC BY 4.0) License https://creativecommons.org/licenses/by/4.0/.
As elastic normal stress exceeds the local shear stress, flow instability and vortex formation occur. The induced elastic stress under EOF not only enhances the instability of the flow but often generates an irregular secondary flow leading to strong disturbance.
(92) It is also vital to consider the effect of the constriction layout of microchannels on the alteration of the field strength within the system. The contraction geometry enhances a larger electric field strength compared with other locations of the channel outside the constriction region, resulting in a higher velocity gradient and stronger extension on the polymer within the viscoelastic solution. Following the high shear flow condition, a higher magnitude of stretch for polymer molecules in viscoelastic fluids exhibits larger elastic stresses and enhancement of vortex formation at the region.
(93)As shown in Figure 6(c), significant elastic normal stress occurs at the inlet of the constriction microchannel. Such occurrence of a polymeric flow can be attributed to the dominating elongational flow, giving rise to high deformation of the polymers within the viscoelastic fluid flow, resulting in higher elastic stress from the polymers. Such phenomena at the entrance result in the difference in velocity streamline as circled in Figure 6(d) compared to that of the Newtonian fluid at the constriction entrance in Figure 6(e).
(90) The difference between the Newtonian and polymer solution at the exit, as circled in Figure 6(d) and (e), can be attributed to the extrudate swell effect of polymers
(94) within the viscoelastic fluid flow. The extrudate swell effect illustrates that, as polymers emerge from the constriction exit, they tend to contract in the flow direction and grow in the normal direction, resulting in an extrudate diameter greater than the channel size. The deformation of polymers within the polymeric flow at both the entrance and exit of the contraction channel facilitates the change in shear stress conditions of the flow, leading to the alteration in streamlines of flows for each region.
4.3. EOF Applications in LOC Systems
4.3.1. Mixing in LOC Systems
Rather than relying on the micromixing controlled by molecular diffusion under low Reynolds number conditions, active mixers actively leverage convective instability and vortex formation induced by electro-osmotic flows from alternating current (AC) or direct current (DC) electric fields. Such adaptation is recognized as significant breakthroughs for promotion of fluid mixing in chemical and biological applications such as drug delivery, medical diagnostics, chemical synthesis, and so on.
(95)Many researchers proposed novel designs of electro-osmosis micromixers coupled with numerical simulations in conjunction with experimental findings to increase their understanding of the role of flow instability and vortex formation in the mixing process under electrokinetic phenomena. Matsubara and Narumi
(96) numerically modeled the mixing process in a microchannel with four electrodes on each side of the microchannel wall, which generated a disruption through unstable electro-osmotic vortices. It was found that particle mixing was sensitive to both the convection effect induced by the main and secondary vortex within the micromixer and the change in oscillation frequency caused by the supplied AC voltage when the Reynolds number was varied. Qaderi et al.
(97) adapted the PNP equation to numerically study the effect of the geometry and zeta potential configuration of the microchannel on the mixing process with a combined electro-osmotic pressure driven flow. It was reported that the application of heterogeneous zeta potential configuration enhances the mixing efficiency by around 23% while the height of the hurdles increases the mixing efficiency at most 48.1%. Cho et al.
(98) utilized the PB model and Laplace equation to numerically simulate the electro-osmotic non-Newtonian fluid mixing process within a wavy and block layout of microchannel walls. The Power Law model is adapted to describe the fluid rheological characteristic. It was found that shear-thinning fluids possess a higher volumetric flow rate, which could result in poorer mixing efficiency compared to that of Newtonian fluids. Numerous studies have revealed that flow instability and vortex generation, in particular secondary vortices produced by barriers or greater magnitudes of heterogeneous zeta potential distribution, enhance mixing by increasing bulk flow velocity and reducing flow distance.To better understand the mechanism of disturbance formed in the system due to externally applied forces, known as electrokinetic instability, literature often utilize the Rayleigh (Ra) number,
(22)where γ is the conductivity ratio of the two streams and can be written as
𝛾=𝜎el,H𝜎el,L�=�el,H�el,L. The Ra number characterizes the ratio between electroviscous and electro-osmotic flow. A high Ra
v value often results in good mixing. It is evident that fluid properties such as the conductivity (σ) of the two streams play a key role in the formation of disturbances to enhance mixing in microsystems. At the same time, electrokinetic parameters like the zeta potential (ζ) in the Ra number is critical in the characterization of electro-osmotic velocity and a slip boundary condition at the microchannel wall.To understand the mixing result along the channel, the concentration field can be defined and simulated under the assumption of steady state conditions and constant diffusion coefficient for each of the working fluid within the system through the convection–diffusion equation as below:
∂𝑐𝒊∂𝑡+∇⇀(𝑐𝑖𝑢⇀−𝐷𝑖∇⇀𝑐𝒊)=0∂��∂�+∇⇀(���⇀−��∇⇀��)=0
(23)where c
i is the species concentration of species i and D
i is the diffusion coefficient of the corresponding species.The standard deviation of concentration (σ
sd) can be adapted to evaluate the mixing quality of the system.
(97) The standard deviation for concentration at a specific portion of the channel may be calculated using the equation below:
m are the non-dimensional concentration profile and the mean concentration at the portion, respectively. C* is the non-dimensional concentration and can be calculated as
𝐶∗=𝐶𝐶ref�*=��ref, where C
ref is the reference concentration defined as the bulk solution concentration. The mean concentration profile can be calculated as
𝐶m=∫10(𝐶∗(𝑦∗)d𝑦∗∫10d𝑦∗�m=∫01(�*(�*)d�*∫01d�*. With the standard deviation of concentration, the mixing efficiency
sd,0 is the standard derivation of the case of no mixing. The value of the mixing efficiency is typically utilized in conjunction with the simulated flow field and concentration field to explore the effect of geometrical and electrokinetic parameters on the optimization of the mixing results.
Viscoelastic fluids such as blood flow in LOC systems are an essential topic to proceed with diagnostic analysis and research through microdevices in the biomedical and pharmaceutical industries. The complex blood flow behavior is tightly controlled by the viscoelastic characteristics of blood such as the dynamic viscosity and the elastic property of RBCs under various shear rate conditions. Furthermore, the flow behaviors under varied driving forces promote an array of microfluidic transport phenomena that are critical to the management of blood flow and other adapted viscoelastic fluids in LOC systems. This review addressed the blood flow phenomena, the complicated interplay between shear rate and blood flow behaviors, and their numerical modeling under LOC systems through the lens of the viscoelasticity characteristic. Furthermore, a theoretical understanding of capillary forces and externally applied electric forces leads to an in-depth investigation of the relationship between blood flow patterns and the key parameters of the two driving forces, the latter of which is introduced through the lens of viscoelastic fluids, coupling numerical modeling to improve the knowledge of blood flow manipulation in LOC systems. The flow disturbances triggered by the EOF of viscoelastic fluids and their impact on blood flow patterns have been deeply investigated due to their important role and applications in LOC devices. Continuous advancements of various numerical modeling methods with experimental findings through more efficient and less computationally heavy methods have served as an encouraging sign of establishing more accurate illustrations of the mechanisms for multiphase blood and other viscoelastic fluid flow transport phenomena driven by various forces. Such progress is fundamental for the manipulation of unique transport phenomena, such as the generated disturbances, to optimize functionalities offered by microdevices in LOC systems.
The following section will provide further insights into the employment of studied blood transport phenomena to improve the functionality of micro devices adapting LOC technology. A discussion of the novel roles that external driving forces play in microfluidic flow behaviors is also provided. Limitations in the computational modeling of blood flow and electrokinetic phenomena in LOC systems will also be emphasized, which may provide valuable insights for future research endeavors. These discussions aim to provide guidance and opportunities for new paths in the ongoing development of LOC devices that adapt blood flow.
5.2. Future Directions
5.2.1. Electro-osmosis Mixing in LOC Systems
Despite substantial research, mixing results through flow instability and vortex formation phenomena induced by electro-osmotic mixing still deviate from the effective mixing results offered by chaotic mixing results such as those seen in turbulent flows. However, recent discoveries of a mixing phenomenon that is generally observed under turbulent flows are found within electro-osmosis micromixers under low Reynolds number conditions. Zhao
(99) experimentally discovered a rapid mixing process in an AC applied micromixer, where the power spectrum of concentration under an applied voltage of 20 V
p-p induces a −5/3 slope within a frequency range. This value of the slope is considered as the O–C spectrum in macroflows, which is often visible under relatively high Re conditions, such as the Taylor microscale Reynolds number Re > 500 in turbulent flows.
(100) However, the Re value in the studied system is less than 1 at the specific location and applied voltage. A secondary flow is also suggested to occur close to microchannel walls, being attributed to the increase of convective instability within the system.Despite the experimental phenomenon proposed by Zhao et al.,
(99) the range of effects induced by vital parameters of an EOF mixing system on the enhanced mixing results and mechanisms of disturbance generated by the turbulent-like flow instability is not further characterized. Such a gap in knowledge may hinder the adaptability and commercialization of the discovery of micromixers. One of the parameters for further evaluation is the conductivity gradient of the fluid flow. A relatively strong conductivity gradient (5000:1) was adopted in the system due to the conductive properties of the two fluids. The high conductivity gradients may contribute to the relatively large Rayleigh number and differences in EDL layer thickness, resulting in an unusual disturbance in laminar flow conditions and enhanced mixing results. However, high conductivity gradients are not always achievable by the working fluids due to diverse fluid properties. The reliance on turbulent-like phenomena and rapid mixing results in a large conductivity gradient should be established to prevent the limited application of fluids for the mixing system. In addition, the proposed system utilizes distinct zeta potential distributions at the top and bottom walls due to their difference in material choices, which may be attributed to the flow instability phenomena. Further studies should be made on varying zeta potential magnitude and distribution to evaluate their effect on the slip boundary conditions of the flow and the large shear rate condition close to the channel wall of EOF. Such a study can potentially offer an optimized condition in zeta potential magnitude through material choices and geometrical layout of the zeta potential for better mixing results and manipulation of mixing fluid dynamics. The two vital parameters mentioned above can be varied with the aid of numerical simulation to understand the effect of parameters on the interaction between electro-osmotic forces and electroviscous forces. At the same time, the relationship of developed streamlines of the simulated velocity and concentration field, following their relationship with the mixing results, under the impact of these key parameters can foster more insight into the range of impact that the two parameters have on the proposed phenomena and the microfluidic dynamic principles of disturbances.
In addition, many of the current investigations of electrokinetic mixers commonly emphasize the fluid dynamics of mixing for Newtonian fluids, while the utilization of biofluids, primarily viscoelastic fluids such as blood, and their distinctive response under shear forces in these novel mixing processes of LOC systems are significantly less studied. To develop more compatible microdevice designs and efficient mixing outcomes for the biomedical industry, it is necessary to fill the knowledge gaps in the literature on electro-osmotic mixing for biofluids, where properties of elasticity, dynamic viscosity, and intricate relationship with shear flow from the fluid are further considered.
5.2.2. Electro-osmosis Separation in LOC Systems
Particle separation in LOC devices, particularly in biological research and diagnostics, is another area where disturbances may play a significant role in optimization.
(101) Plasma analysis in LOC systems under precise control of blood flow phenomena and blood/plasma separation procedures can detect vital information about infectious diseases from particular antibodies and foreign nucleic acids for medical treatments, diagnostics, and research,
(102) offering more efficient results and simple operating procedures compared to that of the traditional centrifugation method for blood and plasma separation. However, the adaptability of LOC devices for blood and plasma separation is often hindered by microchannel clogging, where flow velocity and plasma yield from LOC devices is reduced due to occasional RBC migration and aggregation at the filtration entrance of microdevices.
(103)It is important to note that the EOF induces flow instability close to microchannel walls, which may provide further solutions to clogging for the separation process of the LOC systems. Mohammadi et al.
(104) offered an anti-clogging effect of RBCs at the blood and plasma separating device filtration entry, adjacent to the surface wall, through RBC disaggregation under high shear rate conditions generated by a forward and reverse EOF direction.
Further theoretical and numerical research can be conducted to characterize the effect of high shear rate conditions near microchannel walls toward the detachment of binding blood cells on surfaces and the reversibility of aggregation. Through numerical modeling with varying electrokinetic parameters to induce different degrees of disturbances or shear conditions at channel walls, it may be possible to optimize and better understand the process of disrupting the forces that bind cells to surface walls and aggregated cells at filtration pores. RBCs that migrate close to microchannel walls are often attracted by the adhesion force between the RBC and the solid surface originating from the van der Waals forces. Following RBC migration and attachment by adhesive forces adjacent to the microchannel walls as shown in Figure 7, the increase in viscosity at the region causes a lower shear condition and encourages RBC aggregation (cell–cell interaction), which clogs filtering pores or microchannels and reduces flow velocity at filtration region. Both the impact that shear forces and disturbances may induce on cell binding forces with surface walls and other cells leading to aggregation may suggest further characterization. Kinetic parameters such as activation energy and the rate-determining step for cell binding composition attachment and detachment should be considered for modeling the dynamics of RBCs and blood flows under external forces in LOC separation devices.
Figure 7. Schematic representations of clogging at a microchannel pore following the sequence of RBC migration, cell attachment to channel walls, and aggregation. (105) Reproduced with permission from ref (105). Copyright 2018 The Authors under the terms of the Creative Commons (CC BY 4.0) License https://creativecommons.org/licenses/by/4.0/.
5.2.3. Relationship between External Forces and Microfluidic Systems
In blood flow, a thicker CFL suggests a lower blood viscosity, suggesting a complex relationship between shear stress and shear rate, affecting the blood viscosity and blood flow. Despite some experimental and numerical studies on electro-osmotic non-Newtonian fluid flow, limited literature has performed an in-depth investigation of the role that applied electric forces and other external forces could play in the process of CFL formation. Additional studies on how shear rates from external forces affect CFL formation and microfluidic flow dynamics can shed light on the mechanism of the contribution induced by external driving forces to the development of a separate phase of layer, similar to CFL, close to the microchannel walls and distinct from the surrounding fluid within the system, then influencing microfluidic flow dynamics.One of the mechanisms of phenomena to be explored is the formation of the Exclusion Zone (EZ) region following a “Self-Induced Flow” (SIF) phenomenon discovered by Li and Pollack,
(106) as shown in Figure 8(a) and (b), respectively. A spontaneous sustained axial flow is observed when hydrophilic materials are immersed in water, resulting in the buildup of a negative layer of charges, defined as the EZ, after water molecules absorb infrared radiation (IR) energy and break down into H and OH
+–.
Figure 8. Schematic representations of (a) the Exclusion Zone region and (b) the Self Induced Flow through visualization of microsphere movement within a microchannel. (106) Reproduced with permission from ref (106). Copyright 2020 The Authors under the terms of the Creative Commons (CC BY 4.0) License https://creativecommons.org/licenses/by/4.0/.
Despite the finding of such a phenomenon, the specific mechanism and role of IR energy have yet to be defined for the process of EZ development. To further develop an understanding of the role of IR energy in such phenomena, a feasible study may be seen through the lens of the relationships between external forces and microfluidic flow. In the phenomena, the increase of SIF velocity under a rise of IR radiation resonant characteristics is shown in the participation of the external electric field near the microchannel walls under electro-osmotic viscoelastic fluid flow systems. The buildup of negative charges at the hydrophilic surfaces in EZ is analogous to the mechanism of electrical double layer formation. Indeed, research has initiated the exploration of the core mechanisms for EZ formation through the lens of the electrokinetic phenomena.
(107) Such a similarity of the role of IR energy and the transport phenomena of SIF with electrokinetic phenomena paves the way for the definition of the unknown SIF phenomena and EZ formation. Furthermore, Li and Pollack
(106) suggest whether CFL formation might contribute to a SIF of blood using solely IR radiation, a commonly available source of energy in nature, as an external driving force. The proposition may be proven feasible with the presence of the CFL region next to the negatively charged hydrophilic endothelial glycocalyx layer, coating the luminal side of blood vessels.
(108) Further research can dive into the resonating characteristics between the formation of the CFL region next to the hydrophilic endothelial glycocalyx layer and that of the EZ formation close to hydrophilic microchannel walls. Indeed, an increase in IR energy is known to rapidly accelerate EZ formation and SIF velocity, depicting similarity to the increase in the magnitude of electric field forces and greater shear rates at microchannel walls affecting CFL formation and EOF velocity. Such correlation depicts a future direction in whether SIF blood flow can be observed and characterized theoretically further through the lens of the relationship between blood flow and shear forces exhibited by external energy.
The intricate link between the CFL and external forces, more specifically the externally applied electric field, can receive further attention to provide a more complete framework for the mechanisms between IR radiation and EZ formation. Such characterization may also contribute to a greater comprehension of the role IR can play in CFL formation next to the endothelial glycocalyx layer as well as its role as a driving force to propel blood flow, similar to the SIF, but without the commonly assumed pressure force from heart contraction as a source of driving force.
5.3. Challenges
Although there have been significant improvements in blood flow modeling under LOC systems over the past decade, there are still notable constraints that may require special attention for numerical simulation applications to benefit the adaptability of the designs and functionalities of LOC devices. Several points that require special attention are mentioned below:
1.
The majority of CFD models operate under the relationship between the viscoelasticity of blood and the shear rate conditions of flow. The relative effect exhibited by the presence of highly populated RBCs in whole blood and their forces amongst the cells themselves under complex flows often remains unclearly defined. Furthermore, the full range of cell populations in whole blood requires a much more computational load for numerical modeling. Therefore, a vital goal for future research is to evaluate a reduced modeling method where the impact of cell–cell interaction on the viscoelastic property of blood is considered.
2.
Current computational methods on hemodynamics rely on continuum models based upon non-Newtonian rheology at the macroscale rather than at molecular and cellular levels. Careful considerations should be made for the development of a constructive framework for the physical and temporal scales of micro/nanoscale systems to evaluate the intricate relationship between fluid driving forces, dynamic viscosity, and elasticity.
3.
Viscoelastic fluids under the impact of externally applied electric forces often deviate from the assumptions of no-slip boundary conditions due to the unique flow conditions induced by externally applied forces. Furthermore, the mechanism of vortex formation and viscoelastic flow instability at laminar flow conditions should be better defined through the lens of the microfluidic flow phenomenon to optimize the prediction of viscoelastic flow across different geometrical layouts. Mathematical models and numerical methods are needed to better predict such disturbance caused by external forces and the viscoelasticity of fluids at such a small scale.
4.
Under practical situations, zeta potential distribution at channel walls frequently deviates from the common assumption of a constant distribution because of manufacturing faults or inherent surface charges prior to the introduction of electrokinetic influence. These discrepancies frequently lead to inconsistent surface potential distribution, such as excess positive ions at relatively more negatively charged walls. Accordingly, unpredicted vortex formation and flow instability may occur. Therefore, careful consideration should be given to these discrepancies and how they could trigger the transport process and unexpected results of a microdevice.
Zhe Chen – Department of Chemical Engineering, School of Chemistry and Chemical Engineering, State Key Laboratory of Metal Matrix Composites, Shanghai Jiao Tong University, Shanghai 200240, P. R. China; Email: zaccooky@sjtu.edu.cn
Bo Ouyang – Department of Chemical Engineering, School of Chemistry and Chemical Engineering, State Key Laboratory of Metal Matrix Composites, Shanghai Jiao Tong University, Shanghai 200240, P. R. China; Email: bouy93@sjtu.edu.cn
Zheng-Hong Luo – Department of Chemical Engineering, School of Chemistry and Chemical Engineering, State Key Laboratory of Metal Matrix Composites, Shanghai Jiao Tong University, Shanghai 200240, P. R. China; https://orcid.org/0000-0001-9011-6020; Email: luozh@sjtu.edu.cn
Authors
Bin-Jie Lai – Department of Chemical Engineering, School of Chemistry and Chemical Engineering, State Key Laboratory of Metal Matrix Composites, Shanghai Jiao Tong University, Shanghai 200240, P. R. China; https://orcid.org/0009-0002-8133-5381
Li-Tao Zhu – Department of Chemical Engineering, School of Chemistry and Chemical Engineering, State Key Laboratory of Metal Matrix Composites, Shanghai Jiao Tong University, Shanghai 200240, P. R. China; https://orcid.org/0000-0001-6514-8864
NotesThe authors declare no competing financial interest.
This work was supported by the National Natural Science Foundation of China (No. 22238005) and the Postdoctoral Research Foundation of China (No. GZC20231576).
the field of technological and scientific study that investigates fluid flow in channels with dimensions between 1 and 1000 μm
Lab-on-a-Chip Technology
the field of research and technological development aimed at integrating the micro/nanofluidic characteristics to conduct laboratory processes on handheld devices
Computational Fluid Dynamics (CFD)
the method utilizing computational abilities to predict physical fluid flow behaviors mathematically through solving the governing equations of corresponding fluid flows
Shear Rate
the rate of change in velocity where one layer of fluid moves past the adjacent layer
Viscoelasticity
the property holding both elasticity and viscosity characteristics relying on the magnitude of applied shear stress and time-dependent strain
Electro-osmosis
the flow of fluid under an applied electric field when charged solid surface is in contact with the bulk fluid
Vortex
the rotating motion of a fluid revolving an axis line
1Neethirajan, S.; Kobayashi, I.; Nakajima, M.; Wu, D.; Nandagopal, S.; Lin, F. Microfluidics for food, agriculture and biosystems industries. Lab Chip2011, 11 (9), 1574– 1586, DOI: 10.1039/c0lc00230eViewGoogle Scholar
2Whitesides, G. M. The origins and the future of microfluidics. Nature2006, 442 (7101), 368– 373, DOI: 10.1038/nature05058ViewGoogle Scholar
3Burklund, A.; Tadimety, A.; Nie, Y.; Hao, N.; Zhang, J. X. J. Chapter One – Advances in diagnostic microfluidics; Elsevier, 2020; DOI: DOI: 10.1016/bs.acc.2019.08.001 .ViewGoogle Scholar
4Abdulbari, H. A. Chapter 12 – Lab-on-a-chip for analysis of blood. In Nanotechnology for Hematology, Blood Transfusion, and Artificial Blood; Denizli, A., Nguyen, T. A., Rajan, M., Alam, M. F., Rahman, K., Eds.; Elsevier, 2022; pp 265– 283.ViewGoogle Scholar
5Vladisavljević, G. T.; Khalid, N.; Neves, M. A.; Kuroiwa, T.; Nakajima, M.; Uemura, K.; Ichikawa, S.; Kobayashi, I. Industrial lab-on-a-chip: Design, applications and scale-up for drug discovery and delivery. Advanced Drug Delivery Reviews2013, 65 (11), 1626– 1663, DOI: 10.1016/j.addr.2013.07.017ViewGoogle Scholar
6Kersaudy-Kerhoas, M.; Dhariwal, R.; Desmulliez, M. P. Y.; Jouvet, L. Hydrodynamic blood plasma separation in microfluidic channels. Microfluid. Nanofluid.2010, 8 (1), 105– 114, DOI: 10.1007/s10404-009-0450-5ViewGoogle Scholar
7Popel, A. S.; Johnson, P. C. Microcirculation and Hemorheology. Annu. Rev. Fluid Mech.2005, 37 (1), 43– 69, DOI: 10.1146/annurev.fluid.37.042604.133933ViewGoogle Scholar
8Fedosov, D. A.; Peltomäki, M.; Gompper, G. Deformation and dynamics of red blood cells in flow through cylindrical microchannels. Soft Matter2014, 10 (24), 4258– 4267, DOI: 10.1039/C4SM00248BViewGoogle Scholar
9Chakraborty, S. Dynamics of capillary flow of blood into a microfluidic channel. Lab Chip2005, 5 (4), 421– 430, DOI: 10.1039/b414566fViewGoogle Scholar
10Tomaiuolo, G.; Guido, S. Start-up shape dynamics of red blood cells in microcapillary flow. Microvascular Research2011, 82 (1), 35– 41, DOI: 10.1016/j.mvr.2011.03.004ViewGoogle Scholar
11Sherwood, J. M.; Dusting, J.; Kaliviotis, E.; Balabani, S. The effect of red blood cell aggregation on velocity and cell-depleted layer characteristics of blood in a bifurcating microchannel. Biomicrofluidics2012, 6 (2), 24119, DOI: 10.1063/1.4717755ViewGoogle Scholar
12Nader, E.; Skinner, S.; Romana, M.; Fort, R.; Lemonne, N.; Guillot, N.; Gauthier, A.; Antoine-Jonville, S.; Renoux, C.; Hardy-Dessources, M.-D. Blood Rheology: Key Parameters, Impact on Blood Flow, Role in Sickle Cell Disease and Effects of Exercise. Frontiers in Physiology2019, 10, 01329, DOI: 10.3389/fphys.2019.01329ViewGoogle Scholar
13Trejo-Soto, C.; Lázaro, G. R.; Pagonabarraga, I.; Hernández-Machado, A. Microfluidics Approach to the Mechanical Properties of Red Blood Cell Membrane and Their Effect on Blood Rheology. Membranes2022, 12 (2), 217, DOI: 10.3390/membranes12020217ViewGoogle Scholar
14Wagner, C.; Steffen, P.; Svetina, S. Aggregation of red blood cells: From rouleaux to clot formation. Comptes Rendus Physique2013, 14 (6), 459– 469, DOI: 10.1016/j.crhy.2013.04.004ViewGoogle Scholar
15Kim, H.; Zhbanov, A.; Yang, S. Microfluidic Systems for Blood and Blood Cell Characterization. Biosensors2023, 13 (1), 13, DOI: 10.3390/bios13010013ViewGoogle Scholar
16Fåhræus, R.; Lindqvist, T. THE VISCOSITY OF THE BLOOD IN NARROW CAPILLARY TUBES. American Journal of Physiology-Legacy Content1931, 96 (3), 562– 568, DOI: 10.1152/ajplegacy.1931.96.3.562ViewGoogle Scholar
17Ascolese, M.; Farina, A.; Fasano, A. The Fåhræus-Lindqvist effect in small blood vessels: how does it help the heart?. J. Biol. Phys.2019, 45 (4), 379– 394, DOI: 10.1007/s10867-019-09534-4ViewGoogle Scholar
18Bento, D.; Fernandes, C. S.; Miranda, J. M.; Lima, R. In vitro blood flow visualizations and cell-free layer (CFL) measurements in a microchannel network. Experimental Thermal and Fluid Science2019, 109, 109847, DOI: 10.1016/j.expthermflusci.2019.109847ViewGoogle Scholar
19Namgung, B.; Ong, P. K.; Wong, Y. H.; Lim, D.; Chun, K. J.; Kim, S. A comparative study of histogram-based thresholding methods for the determination of cell-free layer width in small blood vessels. Physiological Measurement2010, 31 (9), N61, DOI: 10.1088/0967-3334/31/9/N01ViewGoogle Scholar
20Hymel, S. J.; Lan, H.; Fujioka, H.; Khismatullin, D. B. Cell trapping in Y-junction microchannels: A numerical study of the bifurcation angle effect in inertial microfluidics. Phys. Fluids (1994)2019, 31 (8), 082003, DOI: 10.1063/1.5113516ViewGoogle Scholar
21Li, X.; Popel, A. S.; Karniadakis, G. E. Blood-plasma separation in Y-shaped bifurcating microfluidic channels: a dissipative particle dynamics simulation study. Phys. Biol.2012, 9 (2), 026010, DOI: 10.1088/1478-3975/9/2/026010ViewGoogle Scholar
22Yin, X.; Thomas, T.; Zhang, J. Multiple red blood cell flows through microvascular bifurcations: Cell free layer, cell trajectory, and hematocrit separation. Microvascular Research2013, 89, 47– 56, DOI: 10.1016/j.mvr.2013.05.002ViewGoogle Scholar
23Shibeshi, S. S.; Collins, W. E. The Rheology of Blood Flow in a Branched Arterial System. Appl. Rheol2005, 15 (6), 398– 405, DOI: 10.1515/arh-2005-0020ViewGoogle Scholar
24Sequeira, A.; Janela, J. An Overview of Some Mathematical Models of Blood Rheology. In A Portrait of State-of-the-Art Research at the Technical University of Lisbon; Pereira, M. S., Ed.; Springer Netherlands: Dordrecht, 2007; pp 65– 87.ViewGoogle Scholar
25Walburn, F. J.; Schneck, D. J. A constitutive equation for whole human blood. Biorheology1976, 13, 201– 210, DOI: 10.3233/BIR-1976-13307ViewGoogle Scholar
26Quemada, D. A rheological model for studying the hematocrit dependence of red cell-red cell and red cell-protein interactions in blood. Biorheology1981, 18, 501– 516, DOI: 10.3233/BIR-1981-183-615ViewGoogle Scholar
27Varchanis, S.; Dimakopoulos, Y.; Wagner, C.; Tsamopoulos, J. How viscoelastic is human blood plasma?. Soft Matter2018, 14 (21), 4238– 4251, DOI: 10.1039/C8SM00061AViewGoogle Scholar
28Apostolidis, A. J.; Moyer, A. P.; Beris, A. N. Non-Newtonian effects in simulations of coronary arterial blood flow. J. Non-Newtonian Fluid Mech.2016, 233, 155– 165, DOI: 10.1016/j.jnnfm.2016.03.008ViewGoogle Scholar
29Luo, X. Y.; Kuang, Z. B. A study on the constitutive equation of blood. J. Biomech.1992, 25 (8), 929– 934, DOI: 10.1016/0021-9290(92)90233-QViewGoogle Scholar
30Oldroyd, J. G.; Wilson, A. H. On the formulation of rheological equations of state. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences1950, 200 (1063), 523– 541, DOI: 10.1098/rspa.1950.0035ViewGoogle Scholar
31Prado, G.; Farutin, A.; Misbah, C.; Bureau, L. Viscoelastic transient of confined red blood cells. Biophys J.2015, 108 (9), 2126– 2136, DOI: 10.1016/j.bpj.2015.03.046ViewGoogle Scholar
32Huang, C. R.; Pan, W. D.; Chen, H. Q.; Copley, A. L. Thixotropic properties of whole blood from healthy human subjects. Biorheology1987, 24 (6), 795– 801, DOI: 10.3233/BIR-1987-24630ViewGoogle Scholar
33Anand, M.; Kwack, J.; Masud, A. A new generalized Oldroyd-B model for blood flow in complex geometries. International Journal of Engineering Science2013, 72, 78– 88, DOI: 10.1016/j.ijengsci.2013.06.009ViewGoogle Scholar
34Horner, J. S.; Armstrong, M. J.; Wagner, N. J.; Beris, A. N. Investigation of blood rheology under steady and unidirectional large amplitude oscillatory shear. J. Rheol.2018, 62 (2), 577– 591, DOI: 10.1122/1.5017623ViewGoogle Scholar
35Horner, J. S.; Armstrong, M. J.; Wagner, N. J.; Beris, A. N. Measurements of human blood viscoelasticity and thixotropy under steady and transient shear and constitutive modeling thereof. J. Rheol.2019, 63 (5), 799– 813, DOI: 10.1122/1.5108737ViewGoogle Scholar
36Armstrong, M.; Tussing, J. A methodology for adding thixotropy to Oldroyd-8 family of viscoelastic models for characterization of human blood. Phys. Fluids2020, 32 (9), 094111, DOI: 10.1063/5.0022501ViewGoogle Scholar
37Crank, J.; Nicolson, P. A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Mathematical Proceedings of the Cambridge Philosophical Society1947, 43 (1), 50– 67, DOI: 10.1017/S0305004100023197ViewGoogle Scholar
38Clough, R. W. Original formulation of the finite element method. Finite Elements in Analysis and Design1990, 7 (2), 89– 101, DOI: 10.1016/0168-874X(90)90001-UViewGoogle Scholar
39Liu, W. K.; Liu, Y.; Farrell, D.; Zhang, L.; Wang, X. S.; Fukui, Y.; Patankar, N.; Zhang, Y.; Bajaj, C.; Lee, J.Immersed finite element method and its applications to biological systems. Computer Methods in Applied Mechanics and Engineering2006, 195 (13), 1722– 1749, DOI: 10.1016/j.cma.2005.05.049ViewGoogle Scholar
40Lopes, D.; Agujetas, R.; Puga, H.; Teixeira, J.; Lima, R.; Alejo, J. P.; Ferrera, C. Analysis of finite element and finite volume methods for fluid-structure interaction simulation of blood flow in a real stenosed artery. International Journal of Mechanical Sciences2021, 207, 106650, DOI: 10.1016/j.ijmecsci.2021.106650ViewGoogle Scholar
41Favero, J. L.; Secchi, A. R.; Cardozo, N. S. M.; Jasak, H. Viscoelastic flow analysis using the software OpenFOAM and differential constitutive equations. J. Non-Newtonian Fluid Mech.2010, 165 (23), 1625– 1636, DOI: 10.1016/j.jnnfm.2010.08.010ViewGoogle Scholar
42Pimenta, F.; Alves, M. A. Stabilization of an open-source finite-volume solver for viscoelastic fluid flows. J. Non-Newtonian Fluid Mech.2017, 239, 85– 104, DOI: 10.1016/j.jnnfm.2016.12.002ViewGoogle Scholar
43Chee, C. Y.; Lee, H. P.; Lu, C. Using 3D fluid-structure interaction model to analyse the biomechanical properties of erythrocyte. Phys. Lett. A2008, 372 (9), 1357– 1362, DOI: 10.1016/j.physleta.2007.09.067ViewGoogle Scholar
44Xu, D.; Kaliviotis, E.; Munjiza, A.; Avital, E.; Ji, C.; Williams, J. Large scale simulation of red blood cell aggregation in shear flows. J. Biomech.2013, 46 (11), 1810– 1817, DOI: 10.1016/j.jbiomech.2013.05.010ViewGoogle Scholar
45Johnson, K. L.; Kendall, K.; Roberts, A. Surface energy and the contact of elastic solids. Proceedings of the royal society of London. A. mathematical and physical sciences1971, 324 (1558), 301– 313, DOI: 10.1098/rspa.1971.0141ViewGoogle Scholar
46Shi, L.; Pan, T.-W.; Glowinski, R. Deformation of a single red blood cell in bounded Poiseuille flows. Phys. Rev. E2012, 85 (1), 016307, DOI: 10.1103/PhysRevE.85.016307ViewGoogle Scholar
47Yoon, D.; You, D. Continuum modeling of deformation and aggregation of red blood cells. J. Biomech.2016, 49 (11), 2267– 2279, DOI: 10.1016/j.jbiomech.2015.11.027ViewGoogle Scholar
48Mainardi, F.; Spada, G. Creep, relaxation and viscosity properties for basic fractional models in rheology. European Physical Journal Special Topics2011, 193 (1), 133– 160, DOI: 10.1140/epjst/e2011-01387-1ViewGoogle Scholar
49Gracka, M.; Lima, R.; Miranda, J. M.; Student, S.; Melka, B.; Ostrowski, Z. Red blood cells tracking and cell-free layer formation in a microchannel with hyperbolic contraction: A CFD model validation. Computer Methods and Programs in Biomedicine2022, 226, 107117, DOI: 10.1016/j.cmpb.2022.107117ViewGoogle Scholar
50Aryan, H.; Beigzadeh, B.; Siavashi, M. Euler-Lagrange numerical simulation of improved magnetic drug delivery in a three-dimensional CT-based carotid artery bifurcation. Computer Methods and Programs in Biomedicine2022, 219, 106778, DOI: 10.1016/j.cmpb.2022.106778ViewGoogle Scholar
51Czaja, B.; Závodszky, G.; Azizi Tarksalooyeh, V.; Hoekstra, A. G. Cell-resolved blood flow simulations of saccular aneurysms: effects of pulsatility and aspect ratio. J. R Soc. Interface2018, 15 (146), 20180485, DOI: 10.1098/rsif.2018.0485ViewGoogle Scholar
52Rydquist, G.; Esmaily, M. A cell-resolved, Lagrangian solver for modeling red blood cell dynamics in macroscale flows. J. Comput. Phys.2022, 461, 111204, DOI: 10.1016/j.jcp.2022.111204ViewGoogle Scholar
53Dadvand, A.; Baghalnezhad, M.; Mirzaee, I.; Khoo, B. C.; Ghoreishi, S. An immersed boundary-lattice Boltzmann approach to study the dynamics of elastic membranes in viscous shear flows. Journal of Computational Science2014, 5 (5), 709– 718, DOI: 10.1016/j.jocs.2014.06.006ViewGoogle Scholar
54Krüger, T.; Holmes, D.; Coveney, P. V. Deformability-based red blood cell separation in deterministic lateral displacement devices─A simulation study. Biomicrofluidics2014, 8 (5), 054114, DOI: 10.1063/1.4897913ViewGoogle Scholar
55Takeishi, N.; Ito, H.; Kaneko, M.; Wada, S. Deformation of a Red Blood Cell in a Narrow Rectangular Microchannel. Micromachines2019, 10 (3), 199, DOI: 10.3390/mi10030199ViewGoogle Scholar
56Krüger, T.; Varnik, F.; Raabe, D. Efficient and accurate simulations of deformable particles immersed in a fluid using a combined immersed boundary lattice Boltzmann finite element method. Computers & Mathematics with Applications2011, 61 (12), 3485– 3505, DOI: 10.1016/j.camwa.2010.03.057ViewGoogle Scholar
57Balachandran Nair, A. N.; Pirker, S.; Umundum, T.; Saeedipour, M. A reduced-order model for deformable particles with application in bio-microfluidics. Computational Particle Mechanics2020, 7 (3), 593– 601, DOI: 10.1007/s40571-019-00283-8ViewGoogle Scholar
58Balachandran Nair, A. N.; Pirker, S.; Saeedipour, M. Resolved CFD-DEM simulation of blood flow with a reduced-order RBC model. Computational Particle Mechanics2022, 9 (4), 759– 774, DOI: 10.1007/s40571-021-00441-xViewGoogle Scholar
60Piquet, A.; Roussel, O.; Hadjadj, A. A comparative study of Brinkman penalization and direct-forcing immersed boundary methods for compressible viscous flows. Computers & Fluids2016, 136, 272– 284, DOI: 10.1016/j.compfluid.2016.06.001ViewGoogle Scholar
61Akerkouch, L.; Le, T. B. A Hybrid Continuum-Particle Approach for Fluid-Structure Interaction Simulation of Red Blood Cells in Fluid Flows. Fluids2021, 6 (4), 139, DOI: 10.3390/fluids6040139ViewGoogle Scholar
62Barker, A. T.; Cai, X.-C. Scalable parallel methods for monolithic coupling in fluid-structure interaction with application to blood flow modeling. J. Comput. Phys.2010, 229 (3), 642– 659, DOI: 10.1016/j.jcp.2009.10.001ViewGoogle Scholar
63Cetin, A.; Sahin, M. A monolithic fluid-structure interaction framework applied to red blood cells. International Journal for Numerical Methods in Biomedical Engineering2019, 35 (2), e3171 DOI: 10.1002/cnm.3171ViewGoogle Scholar
64Freund, J. B. Numerical Simulation of Flowing Blood Cells. Annu. Rev. Fluid Mech.2014, 46 (1), 67– 95, DOI: 10.1146/annurev-fluid-010313-141349ViewGoogle Scholar
65Ye, T.; Phan-Thien, N.; Lim, C. T. Particle-based simulations of red blood cells─A review. J. Biomech.2016, 49 (11), 2255– 2266, DOI: 10.1016/j.jbiomech.2015.11.050ViewGoogle Scholar
66Arabghahestani, M.; Poozesh, S.; Akafuah, N. K. Advances in Computational Fluid Mechanics in Cellular Flow Manipulation: A Review. Applied Sciences2019, 9 (19), 4041, DOI: 10.3390/app9194041ViewGoogle Scholar
67Rathnayaka, C. M.; From, C. S.; Geekiyanage, N. M.; Gu, Y. T.; Nguyen, N. T.; Sauret, E. Particle-Based Numerical Modelling of Liquid Marbles: Recent Advances and Future Perspectives. Archives of Computational Methods in Engineering2022, 29 (5), 3021– 3039, DOI: 10.1007/s11831-021-09683-7ViewGoogle Scholar
68Li, X.; Vlahovska, P. M.; Karniadakis, G. E. Continuum- and particle-based modeling of shapes and dynamics of red blood cells in health and disease. Soft Matter2013, 9 (1), 28– 37, DOI: 10.1039/C2SM26891DViewGoogle Scholar
69Beris, A. N.; Horner, J. S.; Jariwala, S.; Armstrong, M. J.; Wagner, N. J. Recent advances in blood rheology: a review. Soft Matter2021, 17 (47), 10591– 10613, DOI: 10.1039/D1SM01212FViewGoogle Scholar
70Arciero, J.; Causin, P.; Malgaroli, F. Mathematical methods for modeling the microcirculation. AIMS Biophysics2017, 4 (3), 362– 399, DOI: 10.3934/biophy.2017.3.362ViewGoogle Scholar
71Maria, M. S.; Chandra, T. S.; Sen, A. K. Capillary flow-driven blood plasma separation and on-chip analyte detection in microfluidic devices. Microfluid. Nanofluid.2017, 21 (4), 72, DOI: 10.1007/s10404-017-1907-6ViewGoogle Scholar
72Huhtamäki, T.; Tian, X.; Korhonen, J. T.; Ras, R. H. A. Surface-wetting characterization using contact-angle measurements. Nat. Protoc.2018, 13 (7), 1521– 1538, DOI: 10.1038/s41596-018-0003-zViewGoogle Scholar
73Young, T., III. An essay on the cohesion of fluids. Philosophical Transactions of the Royal Society of London1805, 95, 65– 87, DOI: 10.1098/rstl.1805.0005ViewGoogle Scholar
74Kim, Y. C.; Kim, S.-H.; Kim, D.; Park, S.-J.; Park, J.-K. Plasma extraction in a capillary-driven microfluidic device using surfactant-added poly(dimethylsiloxane). Sens. Actuators, B2010, 145 (2), 861– 868, DOI: 10.1016/j.snb.2010.01.017ViewGoogle Scholar
75Washburn, E. W. The Dynamics of Capillary Flow. Physical Review1921, 17 (3), 273– 283, DOI: 10.1103/PhysRev.17.273ViewGoogle Scholar
76Cito, S.; Ahn, Y. C.; Pallares, J.; Duarte, R. M.; Chen, Z.; Madou, M.; Katakis, I. Visualization and measurement of capillary-driven blood flow using spectral domain optical coherence tomography. Microfluid Nanofluidics2012, 13 (2), 227– 237, DOI: 10.1007/s10404-012-0950-6ViewGoogle Scholar
77Berthier, E.; Dostie, A. M.; Lee, U. N.; Berthier, J.; Theberge, A. B. Open Microfluidic Capillary Systems. Anal Chem.2019, 91 (14), 8739– 8750, DOI: 10.1021/acs.analchem.9b01429ViewGoogle Scholar
78Berthier, J.; Brakke, K. A.; Furlani, E. P.; Karampelas, I. H.; Poher, V.; Gosselin, D.; Cubizolles, M.; Pouteau, P. Whole blood spontaneous capillary flow in narrow V-groove microchannels. Sens. Actuators, B2015, 206, 258– 267, DOI: 10.1016/j.snb.2014.09.040ViewGoogle Scholar
79Hirt, C. W.; Nichols, B. D. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys.1981, 39 (1), 201– 225, DOI: 10.1016/0021-9991(81)90145-5ViewGoogle Scholar
80Chen, J.-L.; Shih, W.-H.; Hsieh, W.-H. AC electro-osmotic micromixer using a face-to-face, asymmetric pair of planar electrodes. Sens. Actuators, B2013, 188, 11– 21, DOI: 10.1016/j.snb.2013.07.012ViewGoogle Scholar
81Zhao, C.; Yang, C. Electrokinetics of non-Newtonian fluids: A review. Advances in Colloid and Interface Science2013, 201-202, 94– 108, DOI: 10.1016/j.cis.2013.09.001ViewGoogle Scholar
82Oh, K. W. 6 – Lab-on-chip (LOC) devices and microfluidics for biomedical applications. In MEMS for Biomedical Applications; Bhansali, S., Vasudev, A., Eds.; Woodhead Publishing, 2012; pp 150– 171.ViewGoogle Scholar
83Bello, M. S.; De Besi, P.; Rezzonico, R.; Righetti, P. G.; Casiraghi, E. Electroosmosis of polymer solutions in fused silica capillaries. ELECTROPHORESIS1994, 15 (1), 623– 626, DOI: 10.1002/elps.1150150186ViewGoogle Scholar
84Park, H. M.; Lee, W. M. Effect of viscoelasticity on the flow pattern and the volumetric flow rate in electroosmotic flows through a microchannel. Lab Chip2008, 8 (7), 1163– 1170, DOI: 10.1039/b800185eViewGoogle Scholar
85Afonso, A. M.; Alves, M. A.; Pinho, F. T. Analytical solution of mixed electro-osmotic/pressure driven flows of viscoelastic fluids in microchannels. J. Non-Newtonian Fluid Mech.2009, 159 (1), 50– 63, DOI: 10.1016/j.jnnfm.2009.01.006ViewGoogle Scholar
86Sousa, J. J.; Afonso, A. M.; Pinho, F. T.; Alves, M. A. Effect of the skimming layer on electro-osmotic─Poiseuille flows of viscoelastic fluids. Microfluid. Nanofluid.2011, 10 (1), 107– 122, DOI: 10.1007/s10404-010-0651-yViewGoogle Scholar
87Zhao, C.; Yang, C. Electro-osmotic mobility of non-Newtonian fluids. Biomicrofluidics2011, 5 (1), 014110, DOI: 10.1063/1.3571278ViewGoogle Scholar
88Pimenta, F.; Alves, M. A. Electro-elastic instabilities in cross-shaped microchannels. J. Non-Newtonian Fluid Mech.2018, 259, 61– 77, DOI: 10.1016/j.jnnfm.2018.04.004ViewGoogle Scholar
89Bezerra, W. S.; Castelo, A.; Afonso, A. M. Numerical Study of Electro-Osmotic Fluid Flow and Vortex Formation. Micromachines (Basel)2019, 10 (12), 796, DOI: 10.3390/mi10120796ViewGoogle Scholar
90Ji, J.; Qian, S.; Liu, Z. Electroosmotic Flow of Viscoelastic Fluid through a Constriction Microchannel. Micromachines (Basel)2021, 12 (4), 417, DOI: 10.3390/mi12040417ViewGoogle Scholar
91Zhao, C.; Yang, C. Exact solutions for electro-osmotic flow of viscoelastic fluids in rectangular micro-channels. Applied Mathematics and Computation2009, 211 (2), 502– 509, DOI: 10.1016/j.amc.2009.01.068ViewGoogle Scholar
92Gerum, R.; Mirzahossein, E.; Eroles, M.; Elsterer, J.; Mainka, A.; Bauer, A.; Sonntag, S.; Winterl, A.; Bartl, J.; Fischer, L. Viscoelastic properties of suspended cells measured with shear flow deformation cytometry. Elife2022, 11, e78823, DOI: 10.7554/eLife.78823ViewGoogle Scholar
93Sadek, S. H.; Pinho, F. T.; Alves, M. A. Electro-elastic flow instabilities of viscoelastic fluids in contraction/expansion micro-geometries. J. Non-Newtonian Fluid Mech.2020, 283, 104293, DOI: 10.1016/j.jnnfm.2020.104293ViewGoogle Scholar
94Spanjaards, M.; Peters, G.; Hulsen, M.; Anderson, P. Numerical Study of the Effect of Thixotropy on Extrudate Swell. Polymers2021, 13 (24), 4383, DOI: 10.3390/polym13244383ViewGoogle Scholar
95Rashidi, S.; Bafekr, H.; Valipour, M. S.; Esfahani, J. A. A review on the application, simulation, and experiment of the electrokinetic mixers. Chemical Engineering and Processing – Process Intensification2018, 126, 108– 122, DOI: 10.1016/j.cep.2018.02.021ViewGoogle Scholar
96Matsubara, K.; Narumi, T. Microfluidic mixing using unsteady electroosmotic vortices produced by a staggered array of electrodes. Chemical Engineering Journal2016, 288, 638– 647, DOI: 10.1016/j.cej.2015.12.013ViewGoogle Scholar
97Qaderi, A.; Jamaati, J.; Bahiraei, M. CFD simulation of combined electroosmotic-pressure driven micro-mixing in a microchannel equipped with triangular hurdle and zeta-potential heterogeneity. Chemical Engineering Science2019, 199, 463– 477, DOI: 10.1016/j.ces.2019.01.034ViewGoogle Scholar
98Cho, C.-C.; Chen, C.-L.; Chen, C. o.-K. Mixing enhancement in crisscross micromixer using aperiodic electrokinetic perturbing flows. International Journal of Heat and Mass Transfer2012, 55 (11), 2926– 2933, DOI: 10.1016/j.ijheatmasstransfer.2012.02.006ViewGoogle Scholar
99Zhao, W.; Yang, F.; Wang, K.; Bai, J.; Wang, G. Rapid mixing by turbulent-like electrokinetic microflow. Chemical Engineering Science2017, 165, 113– 121, DOI: 10.1016/j.ces.2017.02.027ViewGoogle Scholar
100Tran, T.; Chakraborty, P.; Guttenberg, N.; Prescott, A.; Kellay, H.; Goldburg, W.; Goldenfeld, N.; Gioia, G. Macroscopic effects of the spectral structure in turbulent flows. Nat. Phys.2010, 6 (6), 438– 441, DOI: 10.1038/nphys1674ViewGoogle Scholar
101Toner, M.; Irimia, D. Blood-on-a-chip. Annu. Rev. Biomed Eng.2005, 7, 77– 103, DOI: 10.1146/annurev.bioeng.7.011205.135108ViewGoogle Scholar
102Maria, M. S.; Rakesh, P. E.; Chandra, T. S.; Sen, A. K. Capillary flow of blood in a microchannel with differential wetting for blood plasma separation and on-chip glucose detection. Biomicrofluidics2016, 10 (5), 054108, DOI: 10.1063/1.4962874ViewGoogle Scholar
103Tripathi, S.; Varun Kumar, Y. V. B.; Prabhakar, A.; Joshi, S. S.; Agrawal, A. Passive blood plasma separation at the microscale: a review of design principles and microdevices. Journal of Micromechanics and Microengineering2015, 25 (8), 083001, DOI: 10.1088/0960-1317/25/8/083001ViewGoogle Scholar
104Mohammadi, M.; Madadi, H.; Casals-Terré, J. Microfluidic point-of-care blood panel based on a novel technique: Reversible electroosmotic flow. Biomicrofluidics2015, 9 (5), 054106, DOI: 10.1063/1.4930865ViewGoogle Scholar
105Kang, D. H.; Kim, K.; Kim, Y. J. An anti-clogging method for improving the performance and lifespan of blood plasma separation devices in real-time and continuous microfluidic systems. Sci. Rep2018, 8 (1), 17015, DOI: 10.1038/s41598-018-35235-4ViewGoogle Scholar
106Li, Z.; Pollack, G. H. Surface-induced flow: A natural microscopic engine using infrared energy as fuel. Science Advances2020, 6 (19), eaba0941 DOI: 10.1126/sciadv.aba0941ViewGoogle Scholar
107Mercado-Uribe, H.; Guevara-Pantoja, F. J.; García-Muñoz, W.; García-Maldonado, J. S.; Méndez-Alcaraz, J. M.; Ruiz-Suárez, J. C. On the evolution of the exclusion zone produced by hydrophilic surfaces: A contracted description. J. Chem. Phys.2021, 154 (19), 194902, DOI: 10.1063/5.0043084ViewGoogle Scholar
108Yalcin, O.; Jani, V. P.; Johnson, P. C.; Cabrales, P. Implications Enzymatic Degradation of the Endothelial Glycocalyx on the Microvascular Hemodynamics and the Arteriolar Red Cell Free Layer of the Rat Cremaster Muscle. Front Physiol2018, 9, 168, DOI: 10.3389/fphys.2018.00168ViewGoogle Scholar
화성 미션 애플리케이션을 위한 NERVA 파생 원자로 냉각수 채널 모델은 1.3m NERVA 파생 원자로(NDR) 냉각수 채널의 전산유체역학(CFD) 연구 결과를 제시합니다. CFD 코드 FLOW-3D는 NDR 코어를 통과하는 기체 수소의 흐름을 모델링하는 데 사용되었습니다. 수소는 냉각제 채널을 통해 노심을 통과하여 원자로의 냉각제 및 로켓의 추진제 역할을 합니다. 수소는 고밀도/저온 상태로 채널에 들어가고 저밀도/고온 상태로 빠져나오므로 압축성 모델을 사용해야 합니다. 기술 문서의 설계 사양이 모델에 사용되었습니다. 채널 길이에 걸친 압력 강하가 이전에 추정한 것(0.9MPa)보다 높은 것으로 확인되었으며, 이는 더 강력한 냉각수 펌프가 필요하고 설계 사양을 재평가해야 함을 나타냅니다.
NERVA-Derived Reactor Coolant Channel Model for Mars Mission Applications presents the results of a computational fluid dynamics (CFD) study of a 1.3m NERVA-Derived Reactor (NDR) coolant channel; The CFD code FLOW-3D was used to model the flow of gaseous hydrogen through the core of a NDR. Hydrogen passes through the core by way of coolant channels, acting as the coolant for the reactor as well as the propellant for the rocket. Hydrogen enters the channel in a high density/low temperature state and exits in a low density/high temperature state necessitating the use of a compressible model. Design specifications from a technical paper were used for the model; It was determined that the pressure drop across the length of the channel was higher than previously estimated (0.9 MPa), indicating the possible need for more powerful coolant pumps and a re-evaluation of the design specifications.
Figure 1 Nuclear Rocket Schematic DiagramFigure 2 Fuel Element – Tip ViewFigure 3 Fuel Element – Tie-Tube Structure (Tie-tubes are black)Figure 5 Three-Dimensional Coolant Channel ModelFigure 6 Two-Dimensional Coolant Channel Model
REFERENCES
Anderson, J. D., Jr., (1990) Modern Compressible Flow, 2d ed., McGraw-Hill, New York. Avallone E. A. and T. Baumeister III, eds., (1987) Mark’s Standard Handbookfor Mechanical Engineers, 9th ed., McGraw-Hill, New York. Bennett, G. L. and T. J. Miller (1992) “Nuclear Propulsion: A Key Transportation Technology for the Exploration of Mars,” Proceedings o f the 9th Symposium on Space Nuclear Power Systems, CONF-920104, M. S. El-Genk and M. D. Hoover, eds., American Institute of Physics, New York, AIP Conference Proceedings No. 246, 2: 383-388. Black, D. L., and S. V. Gunn (1991) “A Technical Summary of Engine and Reactor Subsystem Design Performance during the NERVA Program,” AIAA-91-3450, American Institute of Aeronautics and Astronautics, Washington, D. C. Borowski, S. K., et al. (1992) “Nuclear Thermal Rockets: Key to Moon-Mars Exploration,” Aerospace America, July 1992, pp. 34(5). Borowski, S. K., et al. (1993) “ Nuclear Thermal Rocket/Vehicle Design Options for Future NASA Missions to the Moon and Mars,” AIAA-93-4170, American Institute of Aeronautics and Astronautics, Washington, D. C. Borowski, S. K., et al. (1994) “Nuclear Thermal Rocket/Stage Technology Options for NASA’s Future Human Exploration Missions to the Moon and Mars,” Proceedings o f the 11th Symposium on Space Nuclear Power and Propulsion, CONF-940101, M. S. El-Genk and M. D. Hoover, eds., American Institute of Physics, New York, NY, AIP Conference Proceedings No. 301, 2: 745 – 758. Burmeister, L. C. (1993) Convective Heat Transfer, 2d ed., John Wiley & Sons, New York. Chi, J., R. Holman, and B. Pierce (1989) “Nerva Derivative Reactors for Thermal and Electrical Propulsion,” AIAA-89-2770, American Institute of Aeronautics and Astronautics, Washington, D. C. FIDAP (1993) FIDAP 7.0 User’s Manual, Fluid Dynamics International, Inc. FL0W-3D (1994) FL0W-3D Version 6.0 Quick Reference Guide, Flow Science, Inc., Los Alamos, NM. Hill, P. G. and C. R. Peterson (1970) Mechanics and Thermodynamics o f Propulsion, Addison-Wesley, Reading, MA. Lamarsh, J. R. (1983) Introduction to Nuclear Engineering, 2d ed., Addison-Wesley, Reading, MA. Nassersharif, B. (1991) Notes from a Nuclear Propulsion Short Course, 3-5 January 1991, American Institute of Physics. Nassersharif, B., E. Porta, and D. Hailes (1994) “A Proposal Entitled: Scenario Based Design of Nuclear Propulsion for Manned Mars Mission,” NSCEE, Las Vegas, NV. Shepard, K., et al. (1992) “A Split Sprint Mission to Mars,” Proceedings o f the 9th Symposium on Space Nuclear Power Systems, CONF-920104, M. S. El-Genk and M. D. Hoover, eds., American Institute of Physics, New York, AIP Conference Proceedings No. 246, 1: 58 – 63. Sutton, G. P. (1986) Rocket Propulsion Elements: An Introduction to the Engineering o f Rockets, 5th ed., John Wiley & Sons, New York. U.S. President (1989) “Remarks on the 20th Anniversary of the Apollo 11 Moon Landing July 20, 1989,” Administration o f George Bush, Office of the Federal Register. National Archives and Records Service, 1989, Washington D. C., George Bush, 1989, p. 992. VSAERO (1994) VSAERO User’s Manual E.5, Analytical Methods, Inc., Redmond, WA. White, F. M. (1991) Viscous Fluid Flow, 2d ed., McGraw-Hill, Inc., New York. Zweig, H. R. and M. H. Cooper (1993) “NERVA-Derived Rocket Module for Solar System Exploration,” AIAA-93-2110, American Institute of Aeronautics and Astronautics, Washington, D. C.
aNational Cheng Kung University, Department of Mechanical Engineering, Tainan, Taiwan
bNational Cheng Kung University, Academy of Innovative Semiconductor and Sustainable Manufacturing, Tainan, Taiwan
cJum-bo Co., Ltd, Xinshi District, Tainan, Taiwan
Abstract
워블 전략이 포함된 펄스 레이저 용접(PLW) 방법을 사용하여 알루미늄 및 구리 이종 랩 조인트의 제조를 위한 최적의 가공 매개변수에 대해 실험 및 수치 조사가 수행됩니다. 피크 레이저 출력과 접선 용접 속도의 대표적인 조합 43개를 선택하기 위해 원형 패킹 설계 알고리즘이 먼저 사용됩니다.
선택한 매개변수는 PLW 프로세스의 전산유체역학(CFD) 모델에 제공되어 용융 풀 형상(즉, 인터페이스 폭 및 침투 깊이) 및 구리 농도를 예측합니다. 시뮬레이션 결과는 설계 공간 내에서 PLW 매개변수의 모든 조합에 대한 용융 풀 형상 및 구리 농도를 예측하기 위해 3개의 대리 모델을 교육하는 데 사용됩니다.
마지막으로, 대체 모델을 사용하여 구성된 처리 맵은 용융 영역에 균열이나 기공이 없고 향상된 기계적 및 전기적 특성이 있는 이종 조인트를 생성하는 PLW 매개변수를 결정하기 위해 세 가지 품질 기준에 따라 필터링됩니다.
제안된 최적화 접근법의 타당성은 최적의 용접 매개변수를 사용하여 생성된 실험 샘플의 전단 강도, 금속간 화합물(IMC) 형성 및 전기 접촉 저항을 평가하여 입증됩니다.
결과는 최적의 매개변수가 1209N의 높은 전단 강도와 86µΩ의 낮은 전기 접촉 저항을 생성함을 확인합니다. 또한 용융 영역에는 균열 및 기공과 같은 결함이 없습니다.
An experimental and numerical investigation is performed into the optimal processing parameters for the fabrication of aluminum and copper dissimilar lap joints using a pulsed laser welding (PLW) method with a wobble strategy. A circle packing design algorithm is first employed to select 43 representative combinations of the peak laser power and tangential welding speed. The selected parameters are then supplied to a computational fluidic dynamics (CFD) model of the PLW process to predict the melt pool geometry (i.e., interface width and penetration depth) and copper concentration. The simulation results are used to train three surrogate models to predict the melt pool geometry and copper concentration for any combination of the PLW parameters within the design space. Finally, the processing maps constructed using the surrogate models are filtered in accordance with three quality criteria to determine the PLW parameters that produce dissimilar joints with no cracks or pores in the fusion zone and enhanced mechanical and electrical properties. The validity of the proposed optimization approach is demonstrated by evaluating the shear strength, intermetallic compound (IMC) formation, and electrical contact resistance of experimental samples produced using the optimal welding parameters. The results confirm that the optimal parameters yield a high shear strength of 1209 N and a low electrical contact resistance of 86 µΩ. Moreover, the fusion zone is free of defects, such as cracks and pores.
Fig. 1. Schematic illustration of Al-Cu lap-joint arrangementFig. 2. Machine setup (MFQS-150W_1500WFig. 5. Lap-shear mechanical tests: (a) experimental setup and specimen dimensions, and (b) two different failures of lap-joint welding.
N. Thi Tien et al.
Fig. 9. Simulation and experimental results for melt pool profile. (a) Simulation results for melt pool cross-section, and (b) OM image of melt pool cross-section.
(Note that laser processing parameter of 830 W and 565 mm/s is chosen.).
CrossRefView Record in ScopusGoogle Scholar[11]S. Smith, J. Blackburn, M. Gittos, P. de Bono, and P. Hilton, “Welding of dissimilar metallic materials using a scanned laser beam,” in International Congress on Applications of Lasers & Electro-Optics, 2013, vol. 2013, no. 1: Laser Institute of America, pp. 493-502.
Welding strategies for joining copper and aluminum by fast oscillating, high quality laser beam
High-Power Laser Materials Processing: Applications, Diagnostics, and Systems IX, vol. 11273, International Society for Optics and Photonics (2020), p. 112730C
Experimental investigation on the effect of spot diameter on continuous-wave laser welding of copper and aluminum thin sheets for battery manufacturing
Systematic approach for determining optimal processing parameters to produce parts with high density in selective laser melting process
Int. J. Adv. Manuf. Technol., 105 (10) (2019), pp. 4443-4460 View PDF
CrossRefView Record in ScopusGoogle Scholar[32]A. Ascari, A. Fortunato, E. Liverani, and A. Lutey, “Application of different pulsed laser sources to dissimilar welding of Cu and Al alloys,” in Proceedings of Lasers in Manufacturing Conference (LIM), 2019.
FLOW-3D는 작은 하수 처리 시스템부터 대형 수력 발전 프로젝트까지 수처리 및 환경 산업에 직면한 광범위한 문제를 해결할 수 있는 뛰어난 CFD 소프트웨어 입니다. FLOW-3D는 시뮬레이션의 복잡성을 감소시키고 최적의 솔루션에 대해 노력을 집중할 수 있도록 해줍니다. 이를 통해 통해 파악된 가치 있는 통찰력은 귀하의 상당한 시간과 비용을 절약 할 수 있습니다.
실제 지형을 적용하여 3차원 shallow water hybrid model을 이용한 댐 붕괴 시뮬레이션
FLOW-3D는 자유표면 흐름이 있는 수치해석 알고리듬에 의해 유동의 표면이 시공간적으로 변하는 모사를 위한 이상적인 도구라고 할 수 있습니다. 자유 표면은 물과 공기 같은 높은 비율의 밀도 변화를 가지는 유체들 사이의 특정한 경계를 일컫습니다. 자유 표면 흐름을 모델링하는 것은 일반적인 유동방정식과 난류 모델이 결합된 고급 알고리즘을 필요로 합니다. 이 기능은 FLOW-3D로 하여금 침수 구조에 의해 형성된 방수, 수력 점프 및 수면 변화의 흐름의 궤적을 포착 할 수 있습니다.
September 2023 DOI:10.30955/gnc2023.00436 Conference: 18th International Conference on Environmental Science and Technology CEST2023, 30 August to 2 September 2023, Athens, ...
이 강력한 신제품은 FLOW-3D POST의 기능을 FLOW-3D 제품군 전반으로 확장합니다.
경로 추적 개선 사항 사용자는 이제 즉석에서 경로 추적 재료 속성을 추가, 편집 및 조정할 수 있으므로 복잡한 기술적 결과를 더 많은 청중에게 사실적인 렌더링으로 더욱 쉽게 전달할 수 있습니다.
History data 계산기
이제 FLOW-3D POST 내에서 History data에 대한 수학적 연산이 직접 가능합니다. 사용자는 업데이트된 파이썬 계산기를 이력 데이터에 적용하여 프로브 및 플럭스 표면과 같은 측정 장치의 시계열 데이터에 대한 연산을 간소화할 수 있습니다.
EXODUS 파일 포맷 성능 향상
이 릴리스는 EXODUS 객체의 렌더링을 더 부드럽고 현실적으로 개선합니다.
JSON / EXODUS 파일 다시 로드: 이제 사용자는 시뮬레이션이 실행되는 동안 JSON / EXODUS 파일을 다시 로드하여 시뮬레이션 워크플로우를 중단하지 않고도 진화하는 데이터를 시각화하고 분석할 수 있습니다.
수정된 다공성 시각화: EXODUS 파일 형식의 다공성 출력은 수축 다공성을 해소하고 주조물 내부의 누출 경로를 더 잘 시각화할 수 있도록 수정되었습니다.
EXODUS 출력의 측면 표면은 2025R1 이후 FLOW-3D POST에서 평활화할 수 있습니다.
FLOW-3D WELD 및 FLOW-3D AM 지원
유체, 용융 영역, 열원, 반사 및 입자를 위한 새로운 사전 구성 객체는 FLOW-3D WELD 및 FLOW-3D AM 시뮬레이션의 시각화를 용이하게 합니다. 일반적으로 사용되는 출력의 주석은 FLOW-3D POST에서 결과 파일을 열면 자동으로 제공되므로 후처리 워크플로우가 가속화됩니다.
일반적으로 사용되는 출력을 쉽게 볼 수 있어 사용자가 데이터 해석과 분석에 집중할 수 있습니다.
FLOW-3D POST 2024R1 의 새로운 기능
FLOW-3D POST 2024R1은 EXODUS II 기반의 결과를 확장하여 유체-구조 상호작용 및 열 응력을 시각화할 수 있는 기능을 제공합니다.
또한, 사용자는 이제 삼각형 격자 래스터 및 LandXML 파일을 시각화할 수 있어 모델링 영역을 둘러싼 지형을 더 쉽게 확인할 수 있습니다. 이를 통해 시뮬레이션에 대한 더 나은 컨텍스트를 제공하고, 결과에 집중할 수 있도록 돕습니다.
모델링 영역 내 지형(왼쪽)과 삼각형 지형(오른쪽)의 비교. 모델링 영역에는 산과 하류 계곡이 포함되지 않은 반면, 삼각형 지형은 이를 포함하여 더 우수한 컨텍스트와 명확성을 제공합니다.
마지막으로, 주조 사용자들은 기포 발생에 영향을 받는 지역과 냉각이 필요한 영역을 식별하는 데 도움이 되는 새로운 출력을 보게 될 것입니다.
FLOW-3D POST 2023R2 의 새로운 기능
새로운 결과 파일 형식
FLOW-3D POST 2023R2는 EXODUS II 형식을 기반으로 하는 완전히 새로운 결과 파일 형식을 도입하여 더 빠른 후처리를 가능하게 합니다. 이 새로운 파일 형식은 크고 복잡한 시뮬레이션의 후처리 작업에 소요되는 시간을 크게 줄이는 동시에(평균 최대 5배!) 다른 시각화 도구와의 연결성을 향상시킵니다.
FLOW-3D POST 2023R2 에서 사용자는 이제 선택한 데이터를 flsgrf , EXODUS II 또는 flsgrf 및 EXODUS II 파일 형식 으로 쓸 수 있습니다 . 새로운 EXODUS II 파일 형식은 각 객체에 대해 유한 요소 메쉬를 활용하므로 사용자는 다른 호환 가능한 포스트 프로세서 및 FEA 코드를 사용 하여 FLOW-3D 결과를 열 수도 있습니다. 새로운 워크플로우를 통해 사용자는 크고 복잡한 사례를 신속하게 시각화하고 임의 슬라이싱, 볼륨 렌더링 및 통계를 사용하여 보조 정보를 추출할 수 있습니다.
새로운 결과 파일 형식은 hydr3d 솔버의 성능을 저하시키지 않으면서 flsgrf 에 비해 시각화 작업 흐름에서 놀라운 속도 향상을 자랑합니다.
이 흥미로운 새로운 개발은 결과 분석의 속도와 유연성이 향상되어 원활한 시뮬레이션 경험을 제공합니다.
또한 FLOW-3D POST 2023R2 는 최신 버전의 ParaView로 업그레이드되었으며 ParaView 5.11.1 과 관련된 개선 사항을 제공합니다 .
새로운 시각화 기능
임의의 클립 및 슬라이스를 매끄럽게 만듭니다.
EXODUS II 파일 형식을 사용하면 사용자는 모든 방향에서 부드러운 슬라이스를 생성할 수 있으므로 보고 싶은 대로 정확히 흐름을 시각화하는 것이 더 쉬워집니다.
호형 위어 위의 흐름 방향에 맞춰 정렬된 슬라이스입니다. Surface LIC 표현에서 매끄러운 표면과 유선형을 확인하세요.
모델 출력의 더 나은 정량화
EXODUS II 파일은 체적 개체이므로 흐름의 특성을 더 쉽게 정량화할 수 있습니다. 예를 들어, 아래 표시된 주조 응고 시뮬레이션에서 오른쪽 패널은 히스토그램을 사용하여 주조의 다공성 분포를 설명할 수 있는 방법을 보여줍니다. 마찬가지로 접촉 탱크의 예는 시간이 지남에 따라 소독제 및 병원체 농도 분포가 어떻게 변화하는지 보여주므로 설계 요구 사항이 충족되었는지 여부를 보여주는 데 도움이 됩니다.
향상된 광선 추적
광선 추적은 기술적인 청중과 비기술적인 청중 모두에게 결과를 전달하는 데 유용한 도구이며 EXODUS II 파일 형식에서 사용할 수 있는 체적 데이터는 이 시각화 방법과 잘 작동합니다.
FLOW-3D POST 의 뛰어난 광선 추적 기능을 보여주는 병 채우기 시뮬레이션
Surface LIC로 유동장 표현
새로운 Surface LIC 시각화 도구는 흐름 선단이 함께 모이는 재순환 및 불감대뿐만 아니라 온도, 오염 물질 등의 일반적인 이동을 강조하여 흐름장을 시각화하는 데 도움이 됩니다.
FLOW-3D POST 의 새로운 EXODUS II 파일 형식 및 Surface LIC 표현의 예
애니메이션 유선형
애니메이션 유선형은 표준 보기에서 보기 어려울 수 있는 흐름의 내부 구조에 대한 세부 정보를 시각화하는 데 도움이 됩니다.
FLOW-3D POST 2023R1 의 새로운 기능
FLOW-3D POST 2023R1은 기본 MP4 지원을 갖춘 업데이트된 ParaView 엔진, 쉬운 설치를 위한 자동 종속성 테스트 기능을 갖춘 간소화된 Linux 설치 프로그램, Windows 11 및 RHEL 8 지원을 특징으로 합니다.
단위 표시
단위는 엔지니어링 분석 결과를 해석하고 전달하는 핵심 부분입니다. FLOW-3D POST 2023R1 에서는 단위가 결과 파일에서 자동으로 판독되고 공간 및 히스토리 플롯의 범례에 설정되므로 시뮬레이션 결과를 쉽게 해석하고 전달할 수 있습니다.
자동 PQ 2 플롯
FLOW-3D CAST는 수년 동안 PQ 2 분석을 통해 HPDC 기계 성능에 대한 정보를 제공해 왔으며 이제 FLOW-3D POST 에서 시각화를 지원하도록 이 기능을 확장했습니다. PQ 2 정보는 사전 정의된 플롯에 자동으로 요약되므로 플롯의 가시성을 전환하여 기계가 주조 작업을 수행하는 방식을 확인하기만 하면 됩니다 . 추가적인 이점은 데이터와 시간을 비교하여 압력이 기계 성능을 초과하는 시기를 확인할 수도 있다는 것입니다.
FLOW-3D POST 2022R1은 FLOW-3D 의 포스트 프로세서 에 세 가지 중요한 개발을 제공합니다. 즉, 간소화된 2D 슬라이싱, ParaView의 Python 도구를 사용한 고급 자동화, 향상된 포스트 프로세싱 렌더링 속도입니다.
2D 슬라이싱 기능
2D 슬라이싱 기능이 확장되고 간소화되어 작업이 더욱 간단해지고 강력해졌습니다. FLOW-3D POST 사용자는 이제 슬라이스 표면의 벡터 표현과 여러 색상 변수를 사용하여 2D 슬라이스를 빠르게 생성할 수 있습니다. 이 2분짜리 비디오는 새로운 2D 슬라이스 기능의 예를 제공합니다.
파이썬 도구
2022R1에 ParaView의 Python 도구가 추가되면 FLOW-3D POST 의 자동화 기능이 확장 되어 반복 작업을 자동화하는 매크로는 물론 클릭 한 번으로 전체 결과 세트를 생성하는 일괄 후처리도 포함됩니다. 특정하거나 정교한 유형의 후처리, 시뮬레이션 후 시뮬레이션을 표시하려는 경우 출력을 표준화하고 후처리 작업을 자동화할 수 있는 이러한 새로운 기능을 통해 엄청난 이점을 얻을 수 있습니다.
일괄 후처리를 사용하면 후처리 작업을 사전 정의하는 스크립트 또는 상태 파일을 사용하여 명령줄에서 후처리할 수 있으므로 DOE, 매개변수 스윕 또는 자동화된 워크플로우로 인한 여러 결과 파일에 대한 이미지 및 애니메이션 생성이 용이해집니다. 배치 스크립트 또는 상태 파일을 다양한 결과 파일이나 시뮬레이션 결과 파일의 전체 작업 공간에 적용하여 각 사례에 대해 원하는 출력을 빠르고 일관되게 생성할 수 있습니다. 또한 단일 결과 파일에 대한 일련의 다양한 시각화 출력을 생성하는 데 활용할 수도 있습니다.
우리는 또한 후처리 속도에 대해 연구해 왔으며 FLOW-3D POST 2022R1은 일반적으로 FLOW-3D POST v1.1 보다 10%-30% 더 빠르지 만 정확한 속도 향상은 시뮬레이션 및 출력 세부 사항에 따라 다릅니다. 오른쪽의 몇 가지 예는 성능 향상을 보여줍니다.