part a the experimental technique density gradient centrifugation

• Alerion TM Microbubble Cell Separation System
• Cell Therapy Product Family
• Dead Cell Removal Products
• Red Blood Cell (RBC) Depletion Products
• Streptavidin Microbubbles
• Partnerships
• Applications Overview
• Cell Enrichment for Immunology
• Optimizing Sample Prep for Flow Cytometry
• Explore all Applications
• Leukopak Product
• GMP Quality
• Case Studies
• Improved Efficiency
• Unmatched Scalability
• Leukopak Human T Cell Isolation Kit
• Cloning Mouse B Cells
• RBC Depletion in PBMCs
• Isolation of Rare Cell Types
• RNA Extraction
• Ultimate Guide to Microbubbles
• Articles and Blog
• RBC Depletion for CTC Enrichment
• Cell Separation Techniques
• Cell Separation Overview
• B-Cell Overview
• T-Cell Overview
• RBC Overview
• Video Library
• About Akadeum
• Microbubbles Origin Story
• Latest News
• Our Technology
• Our Advisors
• Partner with Akadeum

Home / Resources / Blog

Density Gradient Centrifugation for Cell Separation

Updated on Nov 28, 2023 Share

Optimizing a cell separation method is critical for many researchers. There are various ways to sort different types of cells, and the best methodology often depends on the cell subset. One of the quicker and more cost-effective ways to sort a sample based only on physical characteristics is using centrifugation.

What Is Density Gradient Centrifugation?

Centrifugation allows researchers to separate cells based on their shape and size. Samples are placed into a centrifuge—a machine designed to spin liquid solutions at a set speed. The rotation causes the mixture to experience a centrifugal force that pushes larger particles from the center toward the bottom and smaller ones to the top. The larger components experience a stronger centrifugal pull force than the smaller components.

The process is similar to density gradient centrifugation, where samples are also placed into a centrifuge, but the end goal is not to sort them by size. Spinning from the centrifuge causes more dense particles to move to the bottom of the tube because these particles have more mass and are carried further by their inertia. Less dense particles then settle higher within the sample. This process creates a sorted solution layered by particle density from least to most, top to bottom.

Principles of Density Gradient Centrifugation

Each particle has a specific set of physical characteristics; the properties of its biological makeup that can be used for separation and isolation. Density gradient centrifugation focuses on two characteristics—size and density.

The length of time required for this process depends upon the particles’ size. Larger particles will reach their position of stability earlier, whereas smaller particles take longer to pass through the larger particle zone and assume a position deeper in the gradient.

Density Gradient Reagents

In density gradient centrifugation, distinct reagents with specific densities are used to isolate or separate cells. Not only can these products speed up the process, but they can also increase purity and throughput. By keeping particles from clustering—creating a distinct interface—the quality of reagents can greatly increase the efficiency of density gradient centrifugation.

Centrifuge Components

The centrifuge’s spinning mechanism is composed of two distinct parts: the rotor assembly and the electric motor. The electric motor provides the rotation power and the rotor assembly transmits the rotation energy into the rotor tubes. The rotor assembly will hold the density gradient tubes, which are prepped with a media of known density to aid in layer separation.

Density Gradient Media

Density gradient centrifugation uses the natural differences in leukocyte density to separate blood cells into specific cell populations. Using a gradient media of known density aids in developing distinct layers and provides natural media barriers between cell populations.

As the motor spins at very high speeds, the chosen density gradient media are pushed upward by centrifugal force through the whole blood solution. Each media has a set density of its own, such as that used to remove red blood cells from whole blood. In this instance, the density is between that of red blood cells and the PBMCs (human peripheral blood mononuclear cells). This centrifugation-facilitated separation depends on selecting an appropriate medium to divide the sample into the desired populations.

Applications of Density Gradient Centrifugation

Gradient centrifugation is used to purify samples into their component populations. While cell separation is a popular application of density gradient centrifugation, it can also be used to determine the unknown densities of different particles.

Centrifugation of all types benefits researchers because it harvests substances for additional experimentation or medical uses. It helps to remove dead cells and debris in samples so that specific groups of cells or particles may be effectively isolated and studied.

Cell separation methods—such as buoyancy-activated cell sorting (BACS™) or FACS (fluorescence-activated cell sorting)—can be supplemented with centrifugation to increase overall efficiency.

Learn more about how BACS microbubbles can be used with centrifugation to increase productivity.

Differential Centrifugation

Differential centrifugation is another centrifugation separation method that is based on a particle’s mass. Since different-sized cells already behave differently, the process is done with no reagent or medium.

Differential centrifugation is sometimes considered a simpler form of centrifugation. It separates cells and organelles while density gradient centrifugation is used for molecules and particles.

The main difference between the two centrifugation methods is the type of physical properties on which the process is based. Differential centrifugation is more straightforward but density gradient centrifugation can sort much smaller particles for increased specificity. An example of differential centrifugation is the preparation of a buffy coat from whole blood.

Difference Between Density Gradient and Differential Centrifugation

Density gradient centrifugation focuses on separating cells based on their density, rather than mass and size alone, allowing scientists to sort cells that share basic sizes, such as white and red blood cells. While differential centrifugation does not require special separation reagents, density gradients utilize density-known medium reagents to aid in forming distinct cell population layers during the centrifugation process.

Density Gradient Centrifugation and Other Cell Isolation Methods

There are many blood separation techniques , such as flow cytometry and single-cell type isolation available. Flow cytometry utilizes fluorescence, size, and granularity of cell populations to distinguish cells. Flow cytometry is ideal for separating populations of cells based on multiple parameters. Single-cell isolation allows researchers to isolate specific cell subsets. Selecting the right separation technique for different desired final cell populations is determined by evaluating each method’s effectiveness .

Microbubbles: A Novel Approach to Immune Cell Separation & Enrichment

Are you running into barriers and limitations caused by inadequate sample preparation using outdated technology and workflows? Stop losing your precious samples!

Akadeum’s microbubble approach to immune cell enrichment aims to maximize yield, reduce processing time, and maintain the health and physiology of immune cells. And we do all of this while eliminating the need for extra equipment like magnets and columns. Akadeum uses a fast, easy workflow that takes place directly in the sample container—no pour-offs, no exposing delicate cells to external forces or rare earth magnets, and no harsh chemicals. Shop all products to see how Akadeum can improve overall yield and maximize efficiency.

Can microbubbles be pelleted by centrifugation? In short, no. At least, not a pellet that…

Fluorescence-activated cell sorting (FACS) is an advanced variant of flow cytometry that leverages fluorescent labels…

Cell isolation—also referred to as cell separation or cell sorting—is the process of isolating one…

View Resources

• Testimonials

Equilibrium Density Gradient Centrifugation in Cesium Chloride Solutions Developed by Matthew Meselson and Franklin Stahl

Matthew Meselson, Franklin Stahl, and Jerome Vinograd, developed cesium chloride, or CsCl, density gradient centrifugation in the 1950s at the California Institute of Technology, or Caltech, in Pasadena, California. Density gradient centrifugation enables scientists to separate substances based on size, shape, and density. Meselson and Stahl invented a specific type of density gradient centrifugation, called isopycnic centrifugation that used a solution of cesium chloride to separate DNA molecules based on density alone. When Meselson and Stahl developed the technique in the mid-1950s, scientists had no other way to separate macromolecules that were of similar size but varied in density. Meselson and Stahl employed their method to determine how DNA replicates, became known as the Meselson-Stahl experiment. Density gradient centrifugation using cesium salts allowed scientists to isolate DNA and other macromolecules by density alone.

Density gradient centrifugation requires the use of a centrifuge, an instrument that spins mixtures in a rotor to concentrate or separate materials. The spinning causes sample solutions in tube or bottle shaped containers to experience a centrifugal force that pushes samples away from the center of the rotor toward the bottom of the tube. Centrifugal force causes components of a mixture to separate by size because larger components experience greater centrifugal force than smaller components. Thus, centrifugal force pushes the larger components of a mixture farther from the rotor and closer to the bottom of the tube.

While many laboratory procedures use conventional centrifuges, density gradient centrifugation requires a special type of centrifuge called an analytical ultracentrifuge or ultracentrifuge. Theodor Svedberg at the Uppsala University in Uppsala, Sweden invented analytical centrifuges in the mid-1920s, which contributed to him winning the Nobel Prize in Chemistry in 1926. Analytical ultracentrifuges can spin samples around five times faster than standard centrifuges. The faster a centrifuge spins, the more force exerted on the sample. The increased force not only provides clearer separation by condensing components of a mixture, and causes solutions to form density gradients during centrifugation. The definition of density is the amount of matter per unit volume of a substance. Dense substances have more matter in a given space than less dance substances. During centrifugation, a density gradient forms in solution when the density of that solution gradually increases moving farther from the centrifuge rotor. In addition to forming density gradients and separating substances, many analytical ultracentrifuges enable scientists to analyze their samples during the process of centrifugation. When performing density gradient centrifugation, researchers use the ability to observe their samples during centrifugation and to watch the order that the components of their samples separate over time.

Before scientists invented density gradient centrifugation in the 1950s, scientists used differential centrifugation to separate substances in a mixture to identify, characterize, and analyze the mixture´s components. Differential centrifugation separates particles of different sizes or shape and isolates them from the liquid in a mixture. Differential centrifugation works well if the particles are greatly different in size. However, scientists could not use differential centrifugation to separate macromolecules of similar sizes or shape, such as DNA molecules.

Several researchers began developing methods to overcome the size and shape problem. Myron Brakke, a protein scientist at the Brooklyn Botanic Garden in Brooklyn, New York, invented the first kind of density gradient centrifugation in 1951. Brakke’s method, called rate zonal centrifugation, separated particles based on size and shape in a relatively dense solution along a density gradient. The solution formed a density gradient for two reasons. First, the centrifugal force pushing down on the solution increased the pressure moving farther from the rotor, it thereby condensed more of the solution in a smaller amount of space. Second, the solution closer to the rotor compressed the solution farther from the rotor. In a dense solution, particles move slow enough for their varying speeds to be measurable. Rate zonal centrifugation used the different speeds at which particles moved away from the rotor, determined by size and shape, to separate components of a solution. However, it could not separate macromolecules based on density alone, because density does not necessarily depend on the size and shape of a particle.

Scientists started to understand DNA structure and knew that genes were made of DNA when Meselson and Stahl conducted their DNA research. However, scientists did not understand how DNA copied itself to pass on the inherited traits stored in the DNA. Meselson and Stahl created their density gradient centrifugation method to better study DNA. In 1953, James Watson and Francis Crick, two scientists at the University of Cambridge in Cambridge, England, proposed a mechanism by which DNA replicated itself. That self-replication of DNA explained how the molecule could pass on its genetic information. Throughout the mid-1950s, scientists proposed different mechanisms to explain the same self-replication phenomenon. Meselson and Stahl used density gradient centrifugation to analyze DNA during replication and test the different replication mechanisms.

The density gradient centrifugation method that Meselson and Stahl developed separated DNA based on density alone. Meselson and Stahl, two postdoctoral fellows at Caltech, hypothesized that they could tag DNA as it replicated and then trace the tag over many replication cycles to see what role the original DNA molecules played in replication. The researchers suggested tagging the original DNA with heavy elements that would not significantly alter the DNA´s chemical properties. While heavy elements would increase the weight of the DNA, they would not increase the size of the molecule, thereby only increasing the density of the original DNA. Meselson and Stahl planned to stop tagging DNA in subsequent replication cycles, which meant that once DNA replicated, the samples would contain some of the original heavy DNA and some light DNA. Therefore, Meselson and Stahl needed to develop a new separation technique capable of separating DNA based on density rather than size or shape.

In developing their new density gradient centrifugation technique, Meselson and Stahl first chose a dense solution where the DNA molecules would separate. The solution needed to be dense enough for the DNA molecules to float before centrifugation. Because Meselson and Stahl studied DNA from viruses and bacteria, the solution also needed to be mild enough so the microbes could survive. In 1956, Meselson and Stahl chose a cesium salt solution. When salt dissolves in a liquid the volume changes slightly but the mass increases resulting in a greater density solution. Because cesium is a heavy element, a cesium salt solution is much denser than the density of most salt solutions and the cesium salt solution did not affect viruses or DNA. Cesium salt solutions had ideal properties for studying DNA.

After picking their dense solution, Meselson and Stahl needed to determine which centrifuge to use for their technique. The researchers used centrifuges manufactured by the Specialized Instrument Corporation, or Spinco, a company founded by Edward Greydon Pickles and Maurice Hanafin in 1946, that became part of Beckman Instruments in 1954. Meselson and Stahl chose the Spinco Model E analytical ultracentrifuge, a model of ultracentrifuge built in 1956. At six feet high and over seven feet long, the Spinco Model E analytical ultracentrifuge spun samples up to 60,000 rotations per minute, which was on the higher end for ultracentrifuges at the time. In addition, the Spinco Model E ultracentrifuge obtained those speeds with very little error and contained a temperature and pressure controlled chamber for the samples. Those qualities ensured that the DNA separated at the right rate and that it retained its structure due to a consistent temperature. For learning how to operate the ultracentrifuge, Meselson and Stahl consulted Jerome Vinograd, the ultracentrifuge expert at Caltech.

When spun rapidly in an ultracentrifuge, cesium chloride spontaneously formed a density gradient essential for Meselson and Stahl´s study of DNA. Because Meselson and Stahl subjected the cesium chloride solution to a large centrifugal force over many hours, the force gradually pushed the cesium chloride salt particles toward the bottom of the tube, which resulted in a higher concentration of cesium chloride farther from the rotor. Cesium Chloride never formed a complete pellet at the bottom of the test tube because of its natural tendency to diffuse, or re-dissolve, back into solution. The solution of the bottom of the tube became slightly denser than the solution at the top. The density of the cesium chloride solution increased along a gradient down the tube.

The cesium chloride density had a density range greater than the difference in densities between the heavy and light DNA that Meselson and Stahl aimed to separate. Therefore, instead of centrifuging until the heavy DNA sank to the bottom of the container and the light DNA rose to the top of the container, Meselson and Stahl centrifuged the samples until the DNA suspended in separate parts of the solution and stopped moving at a point called equilibrium.

Density gradient centrifugation as developed by Meselson and Stahl employed a process called equilibrium sedimentation or isopycnic centrifugation. Equilibrium sedimentation is the process by which particles in a solution reach a point where they reach their isopycnic position and stop moving. Equilibrium sedimentation separates Meselson and Stahl’s centrifugation method from other types of density gradient centrifugation that do not allow particles to reach equilibrium where the density of the particles equals the density of the solution around them. The centrifugal force pushing the particles down equals the force of the solution pushing up, causing the particles to stop moving in the solution. The particles do not reach the bottom of the tube, as in other forms of centrifugation. Meselson and Stahl could determine the density of different DNA molecules based on where they reached equilibrium and stopped moving, because the molecules would stop when their density matched the density of the surrounding solution.

Meselson and Stahl conducted the Meselson-Stahl experiment in 1957 and 1958 with the successful density gradient centrifugation separation of DNA molecules of differing densities. To perform the technique, Meselson and Stahl first added one part lysed Escherichia coli bacteria cells containing heavy and light DNA to seven parts seven molar cesium chloride solution. To prevent the DNA molecules from breaking down in an acidic solution, the researchers maintained a constant level of pH six. Meselson and Stahl then centrifuged the samples for twenty-four hours at twenty-five degrees Celsius and 44,770 rotations per minute. Meselson and Stahl chose the speed of 44,770 rotations per minute for two reasons. First, if they spun the samples too fast, the samples would leak into the rotor chamber. Second, if they spun the samples too slow, equilibrium would take much longer to reach. The speed chosen, 44,770 rotations per minute, balanced the two. After centrifugation, Meselson and Stahl observed the formation of concentrated DNA in the form of dark bands across the sample container. The placement of the bands depended on the density of the DNA. Heavy and light DNA formed separate, distinct bands, with the heavier DNA further down the density gradient.

The particle bands formed during density gradient centrifugation have Gaussian or normal distributions that are symmetric about the mean average value. For the DNA band, the number of DNA molecules greater than the mean density equals the number of DNA molecules less than the mean density. The width of the band is inversely proportional to the centrifugal force. In other words, the faster the centrifuge spins, the sharper, narrower, and more defined the bands become. The bands are sharper because there is more centrifugal force to counteract diffusion. The width of the band relates inversely to the rate the density increases along the density gradient. The steeper the density gradient is, the sharper the band, which allow for more accurate measurements because the densities of the DNA molecules in the band are closer to the mean density. The width of the band relates to the molecular weight and mean density of the molecules in that band. Since Meselson and Stahl knew the approximate molecular weight of their DNA molecules, they used the bandwidth to calculate the densities of the DNA molecules they separated.

Density gradient centrifugation that utilizes equilibrium sedimentation in a cesium chloride solution assisted many scientists in many different experiments throughout the twentieth and twenty-first centuries. Because the method enabled Meselson and Stahl to separate DNA molecules based on density alone, and not particle size or shape, the researchers were able to label parental DNA with denser elements and analyze the distribution of the heavy DNA over many replication cycles. Meselson and Stahl’s findings allowed them to determine how DNA replicates. In 1960, at Caltech, Nobel laureates François Jacob and Sydney Brenner used the same density gradient centrifugation method to show the presence of messenger RNA in cells, which scientists later determined to be an essential molecule for converting genetic information in DNA to expressed proteins. Scientists also later used density gradient centrifugation in a cesium solution to purify plasmids, circular pieces of bacterial DNA used for genetic engineering.

The type of density gradient centrifugation developed by Meselson and Stahl enabled scientists to analyze DNA and other nucleotides like RNA in a way they could not do previously. Unlike other forms of centrifugation, density gradient centrifugation in cesium salts separates molecules based on density only, rather than the size of the molecule or how fast that molecule travels in solution. Molecules like DNA remain undamaged in cesium chloride solutions, thereby enabling scientists to isolate, purify, and analyze those types of molecules. As of 2017, scientists still use density gradient centrifugation with cesium chloride to analyze DNA and RNA.

• The American Chemical Society. "Specialized Instruments Corporation: Rapid Centrifugal Preparation." Analytical Chemistry 26 (1954): 53A.
• Brakke, Myron K. "Density Gradient Centrifugation: A New Separation Technique." Journal of the American Chemical Society 73 (1951): 1847–8.
• Brenner, Sydney, François Jacob, and Mathew Meselson. "An Unstable Intermediate Carrying Information from Genes to Ribosomes for Protein Synthesis." Nature 190 (1961): 576–80.
• Elzen, Boelie. "Scientists and Rotors: The Development of Biochemical Ultracentrifuges." PhD diss., University of Twente, at Enschede, Netherlands, 1988.
• Hinton, Richard, and Miloslav Dobrota. "Density Gradient Centrifugation." Laboratory Techniques in Biochemistry and Molecular Biology 6 (1978): 8–290.
• Holmes, Frederic L. Meselson, Stahl, and the Replication of DNA: A History of "The Most Beautiful Experiment in Biology." New Haven: Yale University Press, 2001.
• Koehler, Christopher S. W. "Developing the Ultracentrifuge" American Chemical Society 2003: 63–6. https://pubs.acs.org/subscribe/archive/tcaw/12/i02/pdf/203chronicles.pdf (Accessed October 17, 2017).
• Price, Carl A. Centrifugation in Density Gradients. N.Y., New York: Academic Press, 1982.
• Meselson, Matthew, and Franklin W. Stahl. "The Replication of DNA in Escherichia Coli." Proceedings of the National Academy of Sciences 44 (1958): 671–82.
• Meselson, Matthew, Stahl, Franklin W., and Jerome Vinograd. "Equilibrium Sedimentation of Macromolecules in Density Gradients." Proceedings of the National Academy of Sciences 43, (1957): 581–8.
• Shackelford, Rodney E. "Molecular Pathology DNA Purification Cesium Chloride (CsCl) Density Gradient Centrifugation." Pathology Outlines 2012. http://www.pathologyoutlines.com/topic/moleculardnapurcesium.html (Accessed October 19, 2016).
• Svedberg, Theodor, and K. O. Pedersen. The Ultracentrifuge, London and New York: Oxford University Press, 1940.
• Watson, James D., and Francis H C Crick. "Genetical Implications of the Structure of Deoxyribonucleic Acid." Nature 171 (1953): 964–7.

How to cite

Articles rights and graphics.

Copyright Arizona Board of Regents Licensed as Creative Commons Attribution-NonCommercial-Share Alike 3.0 Unported (CC BY-NC-SA 3.0)  

Last modified

Share this page.

An official website of the United States government

The .gov means it’s official. Federal government websites often end in .gov or .mil. Before sharing sensitive information, make sure you’re on a federal government site.

The site is secure. The https:// ensures that you are connecting to the official website and that any information you provide is encrypted and transmitted securely.

• Publications
• Account settings

Preview improvements coming to the PMC website in October 2024. Learn More or Try it out now .

• Advanced Search
• Journal List
• Biochem Biophys Rep
• v.36; 2023 Dec
• PMC10709023

Guidelines for an optimized differential centrifugation of cells

Literature reviews reveal a significant deficiency in conceptual comprehension concerning centrifugation, a crucial step in both medical and research protocols. The arbitrary fluctuations in centrifugal forces present a potential threat to the reproducibility of results. To address this, we propose concise guidelines that integrate key factors such as temperature, osmolarity, fluid volume, and viscosity. These guidelines aim to enhance comprehension of optimal sedimentation conditions for cell suspensions. Additionally, we introduce a standardized protocol for determining the optimal RCF and centrifugation time. The goal is to maximize sedimentation efficiency while minimizing cell damage, contributing to a universally applicable and reproducible method in centrifugation practices.

Graphical abstract

• • Lack of conceptual understanding of centrifugation endangers result reproducibility in medical and research procedures.
• • Temperature, osmolarity, fluid volume, and viscosity play crucial roles in achieving optimal sedimentation of cells.
• • An effective protocol for utilizing centrifugal force to attain optimal cell sedimentation is required.

1. Introduction

Centrifugation is a frequently utilized technique in medical and research laboratories for isolating distinct cell populations and cell components from their medium, as well as for precipitating cells in bodily fluids or separating fluids with differing densities by means of a centrifugal force greater than that of gravity, i.e., g. Centrifuges can be classified based on their intended usage (differential versus density gradient), rotor design (fixed-angle, swinging-bucket, and right-angle), or rotational speed, including very low (≤4,000× g ), low (up to 10,000× g ), high (up to 50,000× g ), and ultra (≥50,000× g ) speeds. In ‘differential centrifugation’, dispersions consisting of components with varying densities or sizes are separated based on their sedimentation rates in a density-homogenized suspension. On the other hand, ‘density gradient centrifugation’ separates components into narrow zones that are layered atop each other, depending on their size and mass (rate-zonal) or density (isopycnic). Herein, we focus on differential centrifugation, which is widely employed in separating cells or subcellular fractions in a suspension.

Despite the significant role of differential centrifugation in medical and research procedures, a wide range of centrifugal forces are frequently employed without explicit justification (as outlined in Table 1 ), thereby posing a risk of potential damage or loss of cells that can compromise the reproducibility of results. For instance, shear forces arising from compaction can affect cell surface-sensitive properties, such as cell deposition on surfaces and cell surface charges [ 1 ], and can also lead to a significant reduction in the viability [ 2 ] and virulence of bacterial cells [ 3 ] at higher forces. The outcome is also influenced by the type of resuspension medium and salt concentration [ 2 ]. In addition, high centrifugal force-induced platelet activation has been shown to cause erroneous coagulation test outcomes [ 4 ]. Conversely, the use of low centrifugal force (or time) may result in the loss of cells in the supernatant. These effects are compounded during sequential centrifugation steps, particularly when changes are made to the volume of cell suspensions, temperature, or salt concentration. Optimizing the key parameters in cell centrifugation will lead to improved separation efficiency and higher quality of isolated cellular constituents, with significant implications for advancing cell biology and pharmaceutical research.

Potential discrepancies between recommended and commonly used values for centrifugal force and sedimentation time in cell suspension: Various combinations of values of volume, cell numbers, RCF, and sedimentation time, including scenarios where parameters are minimized or maximized, result in a wide range of experimental conditions. As a consequence, the disparities between these conditions can become noticeably accentuated. The data presented herein is indicative of numerous scientific publications and diverse cellular sources, including ATCC, Sigma, and Invitrogen, among others [ [5] , [6] , [7] ].

Here, we present important factors such as temperature, osmolarity, fluid volume, and viscosity to provide a comprehensive understanding of optimal sedimentation conditions of cell suspensions. Subsequently, we outline a succinct protocol for achieving optimized sedimentation of cell samples in a specific experimental setting.

2. Differential centrifugation

In the process of differential centrifugation ( Fig. 1 A), the duration of sedimentation of a dispersed sample can be determined through the utilization of Eq. (1) : [ 9 ]

• η : viscosity of a suspension ( kg . m − 1 . s − 1 )
• l : pathlength of suspension in a centrifuge tube ( m )
• d : average diameter of dispersion (e.g., cell) ( m )
• ρ 0 and ρ : densities of solvent and dispersion, respectively ( kg . m − 3 )
• G : centrifugal force ( m . s − 2 ) which is the ‘relative centrifugal force’, RCF (unitless) times g (9.8 m . s − 2 )

Experimental principles of cell sedimentation utilizing a centrifugal force. ( A) Fixed-angle (left) and swinging-bucket centrifuges (right) are compared, where differential centrifugation separates dispersions of varying densities or sizes in a homogenized suspension based on their sedimentation rates. The centrifugal force, denoted by G , is displayed using arrows. ( B) Water viscosity at different temperatures is displayed.

The unit of measure ‘revolutions per minute’ (RPM; min − 1 ) is commonly used in scientific practice, albeit erroneously. A more appropriate metric is RCF, which can be determined using nomograms supplied by centrifuge manufacturers, or by utilizing the following equation: RCF = ( 1.118 × 10 − 3 ) × r × RPM 2 , where ‘ r ’ denotes the distance from the central axis of the centrifuge in meters ( Fig. 1 A).

It is important to note that Eq. (1) accounts for the maximum value of r (i.e., l  =  r max ), as illustrated in Fig. 2 . However, particles located closer to the central axis (i.e., those with smaller values of l ) experience reduced centrifugal force, and therefore require a longer time to sediment. And, particles in proximity to the bottom of the centrifuge tube (i.e., those with larger values of l ) experience greater centrifugal force during sedimentation, causing them to be compressed by the tube wall before complete sedimentation is achieved. This phenomenon becomes particularly relevant in experiments involving multiple distinct populations [ 10 ]. Despite this issue, density gradient centrifugation enables the separation of components into narrow zones based on their size, mass, or density, which stack on top of each other. In the present work, we focus on a simpler scenario in which the aforementioned complication can be largely disregarded, and Eq. (1) can be utilized in a simplified format during differential centrifugation.

Regional RCF ( A ) and sedimentation time ( B ) at varying distances from the rotation axis.

3. Impactful parameters in cell centrifugation

The terms relating to solvent and dispersion within Eq. (1) may be categorized as either independent (i.e., ρ 0 , ρ and d ) or dependent (i.e., η ) upon experimental variables. Factors such as temperature, salt concentration, and the volume fraction of the dispersed phase (e.g., cells) can exert an influence upon η , as well as two additional dependent variables - namely, the correction coefficient ( β ) for η , and the zeta potential (ξ) of the cellular surface. The precise effects of such variables are detailed as follows:

3.1. Temperature

Although it is advisable to conduct experiments at lower temperatures (e.g., 4°C) in order to minimize temperature-induced biological effects, such as protein degradation, it is worth noting that alterations in temperature can also have an impact on the centrifugation process, owing to the temperature-dependence of η [ 11 ].

• h : Planck's constant (6.626 × 10 −34  J.s)
• N A : Avogadro's number (6.022 × 10 23  mol −1 )
• V ‾ : average molar volume (18 × 10 −6  m 3 . mol −1 for water)
• Δ G : Gibbs free energy of activation for viscous flow
• R : gas constant (8.314 J K −1 . mol −1 )

As indicated in Fig. 1 B, the information illustrates that with the viscosity of water at 4°C and 25°C being 1.49 g.m −1 . s −1 and 1.11 g.m −1 . s −1 , respectively, and holding all other factors in Eq. (1) steady, a shift in temperature from 4°C to 25°C can result in a 25% alteration in the optimal sedimentation time or RCF required for the sedimentation of a dispersed substance.

3.2. Osmolarity

Biological media and solvents are composed of various concentrations of ions, such as Na + , K + , Ca 2+ , Mg 2+ , NH 4 + , Cl − , SO 4 2− , PO 4 3− etc., that serve to mitigate the effects of osmotic shock. In the context of serial centrifugation, it is conceivable that fluctuations in the concentrations of cations and anions may occur as a result of the use of differing media or buffers. These changes have the potential to influence centrifugation yield via two distinct pathways:

3.2.1. Viscosity

The interplay between ions and water molecules results in alterations to the structure of water, thereby affecting its viscosity. These ions can be classified as either kosmotropes, which are structure makers, or chaotropes, which are structure breakers [ 12 ]. In the realm of low electrolyte concentrations, where the concentration does not exceed 100 mM, the Jones-Dole equation is the prevailing equation utilized to establish a correlation between the concentration dependence of an electrolyte solution's viscosity: [ 13 ]

In Eq. (3) , the symbols η and η 0 represent the dynamic viscosities of the solution and the pure solvent, respectively, while C salt indicates the molar concentration of the electrolyte. The coefficient ‘A’ is linked to the conductance of ions under conditions of infinite dilution, in accordance with the Debye-Huckel theory [ 14 ]. Conversely, the coefficient ‘B’ is determined through experimental data fitting and characterizes the extent of water structuring. For kosmotropic ions, this coefficient possesses a positive value, while chaotropic ions are associated with a negative value [ 15 , 16 ].

The viscosity of electrolyte solutions is known to exhibit anomalous concentration dependence, particularly at higher ionic concentrations [ 12 ]. Notably, for certain electrolytes such as NaCl(aq) or LiCl(aq), the viscosity ( η ) increases monotonically with the concentration of salt ( C salt ). In contrast, for other electrolytes, such as KCl(aq) or KBr(aq), the viscosity initially increases up to a maximum value (at approximately 50–100 mM), then decreases to a minimum value (at approximately 500–750 mM), and finally increases monotonically at higher concentrations [ 17 , 18 ]. Esteves et al. have proposed a more comprehensive model for computing the viscosity of binary strong electrolyte solutions [ 19 ].

In summary, altering the type and/or molar concentration of ions can affect the η value of a suspension, which, in turn, impacts the optimal sedimentation time. Typically, variations in the viscosity of water by 1–5% are observed when the electrolyte concentration of the suspension changes within the range of 100–500 mM [ 17 , 18 ].

3.2.2. Cell-cell repulsion

The surface charge density of prokaryotic [ 20 ] and eukaryotic cells [ 21 ] is contingent upon the physiological state and molecular composition of the cell membrane. Notably, the lipopolysaccharide-coated surface of Gram-negative bacteria exhibits a negative charge density (6.6 ± 1.3 nm −2 for E. coli ) that is up to seven times greater than that of the protein surface layer of Gram-positive cells at physiological pH (1.0 ± 0.2 nm −2 for L. rhamnosus ) [ 22 ]. To determine the cell surface charge, zeta-potential (ξ) is typically utilized, which measures the electrical potential of the interfacial region between the cell surface and the local environment. The ξ values usually range from −50 to −20 mV for different cell types [ 22 , 23 ] or extracellular vesicles [ 24 ] under varying salt and detergent concentrations.

The ξ term serves as an indicator of the repulsive interaction between cells, whereby at zero potential, the conditions for cellular aggregation are maximized. The repulsive force between surface-charged dispersions is roughly proportional to the square of the ξ value (according to Coulomb's law) [ 25 ]. The ionic strength of a suspension has a significant impact on the ξ value, wherein the presence of ions results in a less negative ξ value, particularly for ions with higher valency, such as Ca 2+ [ 1 , [22] , [23] , [24] , 26 ]

To summarize, a higher concentration of salts leads to a more rapid sedimentation of cells. Altering the salinity during serial centrifugation steps may result in the loss of cells in the supernatant or possibly damage the cells due to the application of excessive centrifugal force.

3.3. Volume fraction of cells

In order to account for the influence of the volume fraction of the dispersed phase ( φ , unitless) on the viscosity ( η ) as described in Box 1 , it is necessary to apply a correction coefficient for viscosity, β (unitless), to Eq. (1) [ 27 ]. As φ increases, particularly for lower φ * values, the β value increases. For instance, when the volume fraction of a suspension of eukaryotic cells changes from 10 3  cells.ml −1 ( φ  ≈ 6.5 × 10 −5 ) to 10 4  cells.ml −1 ( φ  ≈ 6.5 × 10 −4 ), the optimal sedimentation time or RCF will increase by 22.5 % (if φ *  = 0.01) or by 2.3 % (if φ *  = 0.1) (see Box 1 ). Therefore, employing constant values for time and RCF during consecutive centrifugation steps will result in either loss of cells (at higher φ values) or excessive centrifugal force on the cells (at lower φ values), which is particularly pronounced for cells with lower φ * values.

The βdependence onφandφ*.

The β coefficient (unitless) incorporates the influence of the volume fraction ( φ ) on the intrinsic viscosity ([η]), where [η] (unitless) is defined as the limit of ( η - η 0 )/( φ η 0 ) as φ approaches zero. Here, η 0 and η represent the viscosities of the pure solvent and the suspension, respectively: [ 27 , 28 ]

Theoretical computations of the Huggins coefficient, denoted as k H (unitless), anticipate a range of values spanning from 5.0 to 5.2 for repulsive spherical particles. Nonetheless, the effects of Brownian motion on diminutive dispersions escalate the predicted values to 6.0 to 6.2 [ 28 ]. The coefficient β is reliant not only on the concentration of the dispersed phase but also on the ratio of the dispersed phase viscosity to that of the pure solvent. As a result, β diverges as the volume fraction of dispersed phase (denoted as φ ) approaches the critical value φ * (described by Eq. (9) ). φ * is interpreted as the volume fraction at which the dispersion undergoes a transition from a fluid to a non-fluid state. The precise value of φ * depends on the particle shape, the propensity for aggregation, and the rate of shear applied to the suspension [ 27 ].

A proposition is put forth for a crossover equation that smoothly transitions between the semi-dilute and concentrated regimes of a suspension, incorporating both Eq. (8) and Eq. (9) . This proposal is put forward with the aim of providing a more systematic and rigorous framework for analyzing the behavior of suspensions in these regimes: [ 27 ]

The term β can be expressed as a function dependent solely on the variables φ and φ * , given that the [ η ] is proportional to φ * . As an instance, [ η ] can be approximated as 1.7/ φ * [ 29 ].

Alt-text: Box 1

4. A standardized protocol for optimizing cellular sedimentation via centrifugation

The absence of a standardized protocol that comprehensively considers all crucial parameters involved in cell centrifugation, coupled with the prevalent practice of disseminating suboptimal protocols that are not tailored to specific cell samples, underscores the imperative for a universal protocol. Accordingly, we propose a succinct yet universally applicable protocol to achieve optimal sedimentation of cell samples within a specified experimental framework. Our proposed protocol to determine optimal centrifugation time and RCF values, involves three steps, as follow:

• Step (1): Assessing the potential impact of the centrifugation on the overall experimental procedure

As depicted in Fig. 3 , the suggested guideline proves beneficial in scenarios involving (i) time-sensitive experiments, requiring a swift completion of the sedimentation step, or (ii) force-sensitive cells where there's a concern that higher centrifugal forces might negatively impact cellular physiology. For situation (i), one should experimentally establish the RCF value for a brief, fixed centrifugation time (e.g., 1 min) that yields the maximum cell pellet (see Fig. 3 A). Conversely, in case (ii), it becomes crucial to ascertain the time value for a fixed, low RCF (e.g., 100), resulting in the maximum cell pellet ( Fig. 3 B). Note that several other experimental variables, including the solvent's viscosity ( η ) and density ( ρ 0 ), cell size ( d ), cell density ( ρ ), and volume fraction ( φ and φ * ), as well as the suspension's length ( l ), temperature, and salt concentration ( C salt ), are all maintained constant.

• Step (2): Experimentally identifying the optimal values for RCF and centrifugation time

A broadly applicable guideline is provided for establishing the ideal RCF in time-sensitive experiments ( A ) or the optimal sedimentation time for cells sensitive to force ( B ). This guideline maintains constant values for solvent and dispersion parameters, including temperature, tube length, salt concentration, and volume fraction of cells. ( C ) The magnitude of the l value will vary based on the volume of the suspension being subjected to centrifugation, which, in turn, depends on the shape of the centrifugation tube (i.e., Eppendorf, conical, or bottle tubes).

After determining the optimal approach between utilizing constant centrifugation time to establish the optimal RCF or determining the optimal centrifugation time by maintaining a constant RCF (Step 1), it is essential to conduct a sample experiment. Fig. 4 presents a straightforward protocol that is applicable to various cell types under different conditions.

A straightforward protocol for achieving optimal sedimentation of cells via centrifugal force. Experimental illustration demonstrating the centrifugation process for a type of ( A ) eukaryotic cells (e.g., Human Embryonic Kidney, HEK cell) and ( B ) bacterial cells (e.g., Pseudomonas aeruginosa ). Schematic representations highlight sedimented cells at the bottom of Eppendorf tubes. The presence of un-sedimented cells grown on LB agar plates (part ( B )), even after 32 min of centrifugation, underscores the challenge of achieving 100 % sedimentation for small-sized bacterial cells.

Fig. 4 A presents an example of employing differential centrifugation to separate 1 ml of 10 6  cells/ml of a force-sensitive eukaryotic cell in an Eppendorf tube. The process involved applying a low RCF of 200 at 4°C, maintaining physiological salt concentration. At various time points, 5 μl of the supernatant was scrutinized using a bright-field transmission microscope. Schematic representations were utilized to enhance the visibility of the pellets, which represent the sedimented cells. These findings indicate that a centrifugation time of 2 min is sufficient to precipitate the majority of cells. Note that this is significantly lower than the times (and RCF) typically employed by users, as documented in Table 1 .

In Fig. 4 B, a similar protocol was followed using a bacterial cell type with a cell density of 10 9  cells/ml. The centrifugation involved an RCF of 2,000 and was conducted under two different osmolarities. The findings indicated that pellets were distinctly visible following 1 min of centrifugation at a salt concentration of 500 mM, whereas at a lower salinity of 50 mM, pellets became apparent only after 2 min (i.e., as expected due to the impact of osmolarity on cell-cell repulsion). 2 μl of the supernatants at each time point of centrifugation were then transferred to an LB agar plate for overnight growth. As shown, even after 32 min of centrifugation, there were still un-sedimented cells in the supernatant that grew on the plate. This suggests that achieving 100% sedimentation is challenging for small-sized cells like bacteria, even with prolonged centrifugation. Consequently, it is highly recommended to avoid extended periods (or high RCF) in such cases.

• Step (3): To obtain the optimal RCF or centrifugation time values under diverse conditions

Upon conducting experimental analysis of a particular cell sample, the optimal values of RCF and sedimentation time (RCF 1 and t 1 ) can be determined. Subsequently, utilizing Eq. (4) or (5) enables the acquisition of optimum RCF or time values across varying parameters where the cell type and solvent remain constant, but the temperature, salt concentration, volume fraction, or suspension's length may vary. Both equations are derived from Eq. (1) under two distinct conditions represented by subscripts 1 and 2. Note that terms such as 6 π , d 2 , and ( ρ − ρ 0 ) in Eq. (1) are eliminated as their values remain unchanged under both conditions.

It is worth noting that the l value varies differently with the volume of a suspension, depending on the centrifuge tube type (see Fig. 3 C).

As already mentioned, altering influential parameters such as temperature, salinity, volume fraction of cells, or even the centrifugation tube can have a substantial impact on the optimized values of RCF and sedimentation time. It is important to note that determining the exact influence of those variables on viscosity of cell suspensions, zeta-potential of cells, and β parameter can be challenging and time-consuming. Consequently, it is strongly advised to maintain a consistent set of values for these parameters throughout a serial centrifugation process. By doing so, the equations (4) , (5) can be simplified to equations (6) , (7) as follows:

As an illustration, in a time-sensitive experiment, if an optimized certification time of 1 min and an RCF of 2,500 are determined for a specific cell sample, extending the length of the cell suspension by two times and applying a sedimentation time of 30 s would necessitate an RCF of 10,000, assuming that all other parameters remain unaltered (refer to Eq. (7) ). Alternatively, in the context of optimizing certification time, if a specific force-sensitive cell sample exhibits an ideal configuration when subjected to a 5-min centrifugation with an RCF of 300, by increasing the duration of the cell suspension twofold and employing an RCF of 200, the optimal sedimentation time extends to 15 min, assuming all other parameters remain unchanged (refer to Eq. (6) ).

5. Conclusions

A diverse range of centrifugal forces are often employed without explicit justification ( Table 1 ), putting cells at risk of potential damage or loss. Shear forces, which are triggered by compaction, can impact cell surface-sensitive properties. Higher forces can also lead to significant reductions in viability of cells. The outcome is also affected by the type of resuspension medium and salt concentration. Conversely, utilizing low RCF (or time) may result in cell loss in the supernatant. These effects are amplified during serial centrifugation steps, particularly when changes are made to the volume of cell suspensions, temperature, or salt concentration. For example, if one were to switch from utilizing 5 ml, 8°C, and higher salt concentrations (e.g., ξ ≈ −50 mV) to 10 ml, 4°C, and lower salt concentrations (e.g., ξ ≈ −35 mV), the optimal sedimentation time or RCF could differ by up to 400% (according to Eq. (4) or (5)). Therefore, utilizing an unoptimized centrifugation method has the potential to generate non-reproducible experimental results, a crucial aspect that is frequently overlooked by users. Optimizing the influential parameters in cell centrifugation will enhance separation efficacy and improve the quality of isolated cellular constituents, with significant implications for advancing cell biology and pharmaceutical research.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by NIH R01 EY026171 and EY032956 (GWL).

The Isolation of Satellite DNA by Density Gradient Centrifugation

Series: Methods In Molecular Biology > Book: Nucleic Acids

Protocol | DOI: 10.1385/0-89603-064-4:21

• Department of Biological Chemistry, School of Medicine, University of California, Davis, CA, USA

The term satellite DNA is used for a DNA component that gives a sharp band in a density gradient and can be resolved from the broader main band of DNA in the gradient. The usual gradient material is CsCl in aqueous buffer and the Cs + ions form a

The term satellite DNA is used for a DNA component that gives a sharp band in a density gradient and can be resolved from the broader main band of DNA in the gradient. The usual gradient material is CsCl in aqueous buffer and the Cs + ions form a density gradient in a centrifugal field. DNA in the solution sediments to its isopycnic point. The density of DNA is a function of base composition and sequence and so a homogeneous or highly repeated DNA sequence will form a sharp band in CsCl density gradients at a characteristic density. The resolution of this procedure may be enhanced or modified by binding ligands to the DNA. For example, netropsin binds specifically to A + T-rich regions of DNA and reduces their density (1,2). Another useful ligand is Ag + , which must then be centrifuged in Cs 2 SO 4 gradients to avoid precipitation of AgCl (3). Pharmacia has recently introduced CsCF 3 COO as a gradient material.

Figures ( 0 ) & Videos ( 0 )

Experimental specifications, other keywords.

Related articles

Opus-dsd: deep structural disentanglement for cryo-em single-particle analysis, cesium chloride-bisbenzimide gradients for separation of phytoplasma and plant dna.

Construction of Phage Mutants

In vitro base excision repair assay using mammalian cell extracts, purification of peripheral blood natural killer cells, the wheat germ cell-free expression system, dendritic cell generation from highly purified cd14+monocytes, methods for growth and purification of enteric adenovirus type 40, nucleotide excision repair assay in drosophila melanogaster using established cell lines.

• Matthews, H. R., Johnson, E. M, Steer, W. M., Bradbury, E. M., and Allfrey, V. G. (1978) The use of netropsin with CsCl gradients for the analysis of DNA and its application to restriction nuclease fragments of ribosomal DNA from Physarum polycephalum . Eur. J. Biochem. 82 , 569–576.
• Matthews, H. R., Pearson, M. D., and Maclean, N. (1980) Cat satellite DNA: isolation using netropsin with CsCl gradients. Biochim. Biophys. Acta 606 , 228–235.
• Jensen, R. H., and Davidson, N. (1966) Spectrophotometric, potentiometric and density gradient ultracentrifugation studies of the binding of silver ion by DNA. Biopolymers 4 , 17–32.
• Seebeck, T., Stalder, J., and Braun, R. (1979) Isolation of a minichromosome containing the ribosomal genes from Physarum polycephalum . Biochemistry 18 , 484–490.
• Hudspeth, M. E. S., Shumard, D. S., Tatti, K. M., and Grossman, L. I. (1980) Rapid purification of yeast mitochondrial DNA in high yield. Biochim. Biophys. Acta 610 , 221–228.
• Gould, H., and Matthews, H. R. (1976) Separation Methods for Nucleic Acids and Oligonudeotides , Elsevier/North Holland, Amsterdam.
• Judelson, H. S., and Vogt, V. M. (1982) Accessibility of ribosomal genes to trimethyl psoralen in nuclei of Physarum polycephalum . Mol. Cell. Biol. 2 , 211–220.

Advertisement

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

• View all journals
• Explore content
• About the journal
• Publish with us
• Sign up for alerts
• Open access
• Published: 22 August 2024

Structure of anellovirus-like particles reveal a mechanism for immune evasion

• Shu-hao Liou 1   na1   nAff2 ,
• Rajendra Boggavarapu 1   na1 ,
• Noah R. Cohen 1   nAff3 ,
• Yue Zhang   ORCID: orcid.org/0000-0003-4367-4470 1 ,
• Ishwari Sharma 1 ,
• Lynn Zeheb 1 ,
• Nidhi Mukund Acharekar 1 ,
• Hillary D. Rodgers 1 ,
• Saadman Islam 1   nAff4 ,
• Jared Pitts 1 ,
• Cesar Arze   ORCID: orcid.org/0000-0002-3606-4183 1 ,
• Harish Swaminathan 1   nAff5 ,
• Nathan Yozwiak 1   nAff6 ,
• Tuyen Ong 1 ,
• Roger J. Hajjar 1   nAff6 ,
• Yong Chang 1 ,
• Kurt A. Swanson   ORCID: orcid.org/0000-0003-4749-8355 1 &
• Simon Delagrave   ORCID: orcid.org/0000-0003-0865-2196 1   nAff7  

Nature Communications volume  15 , Article number:  7219 ( 2024 ) Cite this article

Metrics details

• Cryoelectron microscopy
• Drug delivery
• Virus structures

Anelloviruses are nonpathogenic viruses that comprise a major portion of the human virome. Despite being ubiquitous in the human population, anelloviruses (ANVs) remain poorly understood. Basic features of the virus, such as the identity of its capsid protein and the structure of the viral particle, have been unclear until now. Here, we use cryogenic electron microscopy to describe the first structure of an ANV-like particle. The particle, formed by 60 jelly roll domain-containing ANV capsid proteins, forms an icosahedral particle core from which spike domains extend to form a salient part of the particle surface. The spike domains come together around the 5-fold symmetry axis to form crown-like features. The base of the spike domain, the P1 subdomain, shares some sequence conservation between ANV strains while a hypervariable region, forming the P2 subdomain, is at the spike domain apex. We propose that this structure renders the particle less susceptible to antibody neutralization by hiding vulnerable conserved domains while exposing highly diverse epitopes as immunological decoys, thereby contributing to the immune evasion properties of anelloviruses. These results shed light on the structure of anelloviruses and provide a framework to understand their interactions with the immune system.

Introduction

Anelloviruses, the principal constituent of the commensal human virome, are universally acquired in infancy and found throughout the body 1 , 2 , 3 . Since the discovery of the original torque teno virus in 1997 4 , three highly diverse genera of human Anelloviridae have been described. These viruses elicit weak immune responses that permit multiple strains to co-exist and persist for years in a typical individual 5 . However, because they do not cause disease 6 and due to the lack of an in vitro culture system, the structural elements of the anellovirus capsid remain poorly understood 7 , 8 .

Anelloviruses comprise a single-stranded, negative-sense circular DNA genome ranging from ~2.0–3.9 kilobases (kb). Among the three genera of human anelloviruses, the Alphatorquevirus (torque teno virus; TTV) genome is ~3.8 kb, the Gammatorquevirus (torque teno midi virus; TTMDV) genome is ~3.2 kb, and the Betatorquevirus (torque teno mini virus; TTMV) genome is ~2.9 kb. Through alternative splicing, at least six different open reading frames (ORFs) are encoded by these genomes, the longest one being ORF1 9 . A distinguishing feature of ORF1 is the presence at the N-terminus of several arginine residues and some lysines which are consistently observed in anelloviruses of all genera. This arginine-rich motif (ARM) is seen in several viral families and is thought to function as a nuclear localization signal as well as to bind the viral genome 10 , 11 . ARM-containing viruses typically contain a jelly roll domain characteristic of viral capsid proteins 12 . Therefore, we hypothesized that the ORF1 protein also contains a jelly roll domain. Recently, Butkovic et al. sequenced >250 ANV genomes from Weddell seal (Leptonychotes weddellii) and grizzly bear (Ursus arctos horribilis) samples and used sequence similarity detection and structural prediction to suggest ORF1 would form a canonical jelly roll structure like those found in circular single-stranded DNA (ssDNA) viruses such as circoviruses 13 . However, no experimental evidence was available to confirm the presence of such a domain in the anellovirus ORF1. Indeed, a few failed attempts at expressing full-length ORF1 have been reported by others 9 , limiting characterization of this protein.

Here we describe the expression, purification, and structural characterization of the anellovirus virus-like particle generated in two different cell lines. Initial expression of the ORF1 protein in insect cells (Sf9) suggested a proteolytically susceptible region of ORF1 that resulted in the removal of the C-terminus. With the C-terminus of the ORF1 polypeptide removed, the protein spontaneously formed symmetric 60-mer particles, as determined by cryogenic electron microscopy (cryo-EM). To generate a higher resolution structure, we generated an ORF1 polypeptide lacking the C-terminus (ΔC-Term) in human cells (Expi 293) and determined its structure. Particles from both expression systems adopt comparable folds, including an icosahedral particle core. Extending from the particle are two subdomains, herein referred to as P1 and P2, that form a spike structure around the particle’s 5-fold axis. The spike domains place the P2 subdomain, comprised primarily of the hypervariable region (HVR), on the outermost surface of the particles. We therefore propose that anelloviruses have evolved this structural orientation to display their highly diverse sequences on the particle surface as a mechanism for immune evasion that potentially explains the weak immune responses to anelloviruses detected in humans.

Anellovirus LY1 ORF1 forms virus-like particles when the C-terminus is proteolyzed in insect cells

Our initial efforts to study the structure of an ANV particle, derived from a Betatorquevirus isolate called LY1 14 , were performed with the full-length ORF1 (residues 1–672) expressed in insect cells (Sf9) using a baculovirus system (pFastbac, Thermo Fischer) encoding a codon-optimized ORF1 construct (Fig.  1 ). Full-length ORF1 assembled into particles ~32 nm in diameter (Fig.  1D ), similar to previously reported estimates of anellovirus size 15 . However, full-length ORF1 particles lacked the structural homogeneity and symmetry expected of viral particles. Furthermore, we noticed that ORF1 tended to degrade in the cells or during purification. Genome packaging by ARM-containing viruses is believed to occur when the positive arginine residues bind a negatively charged genome, overcoming electrostatic repulsion between ARM domains to permit capsid assembly. To prevent potential proteolysis from occurring in the ARM by trypsin-like proteases, we designed the LY1 ΔARM construct wherein residues 2–45 are deleted (Fig.  1A, B ). LY1 ΔARM expression was confirmed using a rabbit polyclonal antibody raised against a spike P1 domain peptide (residues 485–502; Fig.  1B ). We observed a band consistent with LY1 ΔARM by western blot above the 62 kDa marker (expected mass of 73.3 kDa; Fig.  1C ). However, we continued to observe a smaller ORF1 band accumulating below the 62 kDa marker during purification, confirming the proteolysis was not limited to the ARM region. To identify the region of proteolysis, we generated polyclonal antibodies to peptides from the jelly roll domain at the extreme N-terminus (residues 46–58) and the C-terminal domain at the C-terminus (residues 635–672) of the LY1 ΔARM construct (Fig.  1B ) and confirmed the presence of the N-terminal peptide of the LY1 ΔARM fragment and the absence of the C-terminal peptide by western blot analysis (Fig.  1C ).

A A schematic representation of full-length LY1 ORF1 is shown as a cartoon labeled and colored by domains. The arginine-rich motif (ARM) is shown in purple, the jelly roll (JR) domain is shown in red, the spike P1 domain is shown in blue, the spike P2 domain is shown in green, and the C-terminal domain is shown in cyan. Residue numbers beginning each domain and the C-terminus are indicated above. B The sequence of full-length LY1 ORF1 colored as in A with residue numbers indicated above. A dashed line above the sequence indicates residues removed in either the insect or mammalian ORF1 constructs as indicated. Secondary structure elements are indicated above with β-strands as arrows and α-helices as zig-zag lines. The jelly roll β-strands are labeled B-I per convention while additional secondary structures are numbered by their domain. Three peptides used to generate polyclonal antibodies are underlined. C Western blot analysis of LY1 ΔARM after expression (Expression) and after purification (Purification) in insect cells. A molecular weight marker is labeled to the left of the gels, while arrows on the right indicate the band of LY1 ΔARM before (LY1 ΔARM) and after proteolysis (LY1 ΔARM Fragment). Polyclonal antibodies used to probe the western blots, colored and named by the peptides used to generate them ( A ), are indicated below. D – F A cartoon representing the primary structure of the construct is above each micrograph, showing full-length ORF1 contains a complete ARM and C-terminal region whereas the ΔARM constructs lack the ARM region and have truncated C-terminal regions. D Full-length LY1 shows highly heterogeneous particles by negative staining EM. E LY1 ΔARM virus-like particles have a more symmetrical appearance after C-terminal proteolysis. F Genetic truncation of the C-terminal (Δ552-672) preserves a symmetrically structured particle.

EM analysis of the LY1 ΔARM fragment showed that the particles formed were more homogeneous and symmetric in morphology relative to the particles formed by full-length ORF1 (Fig.  1D, E ). The formation of homogeneous particles following genetic removal of the N-terminus and proteolysis of the C-terminus has been observed in another jelly roll-containing virus, hepatitis E (HEV) 16 . To determine if the C-terminus of LY1 ORF1 is required for particle formation, we generated virus-like particles from construct LY1 ORF1 ΔARM ΔC-Term (Δ2-45 and Δ552–672). LY1 ORF1 ΔARM ΔC-Term produced particles of similar symmetry to LY1 ORF1 ΔARM fragment (Fig.  1F ). These results suggest that proteolysis or removal of the ORF1 C-terminus improved particle formation.

ΔC-Term virus-like particle from human cells

To generate an increased number of particles to permit higher resolution structure determination, we generated an ARM-containing ΔC-Term ORF1 construct (residues 1-609) for expression in Expi 293 cells (Fig.  2 ). Western blot and Coomassie-stained sodium dodecyl-sulfate polyacrylamide gel electrophoresis (SDS-PAGE) analysis of the ΔC-Term ORF1 construct purified from mammalian cells showed a principal band consistent with the expected size (Fig.  2B–D ). Despite containing the N-terminal ARM region, the ΔC-Term ORF1 formed symmetrical particles similar to the ΔARM fragments produced in insect cells, suggesting the presence of the ARM region does not hinder particle formation (Fig.  2E ).

A Schematic representation of ΔC-term LY1 ORF1 is shown as a cartoon labeled and colored by domain, colored as in Fig.  1 . Residue numbers beginning each domain and indicating the C-terminal deletion indicated above. B Western blot analysis of ΔC-term LY1 ORF1 expressed in mammalian cells. C Coomassie-stained SDS-PAGE displaying the high purity of the ΔC-term LY1 ORF1. D Western blot analysis shows the purified ΔC-term LY1 ORF1. E Negatively stained electron micrograph shows the homogeneous population of ΔC-term LY1 ORF1 particles.

Cryo-EM structure of the ANV particle

We determined the cryo-EM structure of the virus-like particles formed by the insect cell-derived ΔARM fragment to 3.9 Å resolution from 6271 particles (Supplementary Figs.  1 and 2 , Supplementary Table  1 ) and of the mammalian cell-derived LY1 ΔC-Term ORF1 particle using cryo-EM to 2.69 Å resolution from 22,743 particles (Fig.  3 , Supplementary Figs.  3 , 4 and 5 , Supplementary Table  2 ). The ORF1 particle is formed by sixty ORF1 fragments organized in an icosahedral T  = 1 symmetry. Despite being generated in two different cell lines and the presence of the ARM region in the mammalian virus-like particle (which is not visible in the density, suggesting flexibility) the two protomer structures are very similar (1.179 Å root mean square deviation when aligned over residues 48–562; Fig.  3E ). 7 β−strands within residues 46–228 (named β−strands B to H by convention) form the majority of the canonical 8-β−strand jelly roll domain, and the eighth β−strand (strand I) is formed by residues 531-542 (Fig.  1A, B , Fig.  3F ). Residues 229-530 are inserted between the H and I β−strands (herein called the H-I insertion) and form the exterior of the particle surface. Insertions between jelly roll β−strands are present in other viruses such as adeno-associated virus (AAV) and canine parvovirus (CPV) 17 , 18 . However, the insertions for AAV2 (228 residues) and CPV (227 residues) are between G and H β−strands while the 298-residue H-I insertion in LY1 is significantly larger. The H-I insertion extends from the jelly roll domain to form a structure herein referred to as a spike domain. The spike domain is formed by two subdomains; the spike P1 subdomain (residues 229–250 and 386–530) and the spike P2 subdomain (residues 251–385).

A The cryo-EM map demonstrating the local resolution of the SF9 cell expression purified LY1 ΔARM ORF1 particle colored by its resolution (shown on the scale to the right). B . The cryo-EM map marking the local resolution of expi293-expression system purified LY1 ΔC-Term particle colored by its resolution (same scale as in A ). C A ribbon representation of the LY1 ΔARM 60-mer particle atomic model. D A ribbon representation of the LY1 ΔC-Term atomic model shown as in ( C ). E An overlay of the SF9 cell expression purified (red) and mammalian cell-derived (green) protomers. The observable N- and C-termini and jelly roll, P1 and P2 domains are labeled. F One ΔC-Term ORF1 protomer shown in its electron density with domains labeled and colored as above.

Jelly roll domains form particle core

Sixty LY1 jelly roll domains form the core of the virus particle (Fig.  4 ). The β−strands form β −sheets which are characterized by a C-H-E-F pattern on the core’s exterior and B-I-D-G pattern on the core’s interior. The N-terminus of strand B is oriented to place the ARM on the interior of the core, where it is positioned to bind the viral genome. The observed C-terminal residues (545-562) extend from the C-terminus of β-strand I on the interior of the particle and thread through jelly roll domains on the particle exterior near the 2-fold axis contacting the neighboring jelly roll domain (Fig.  4B ). In several jelly roll-containing viruses, including the avian circovirus beak and feather disease virus (BFDV) 11 , positively charged residues (arginine and lysine) oriented internally on strands B, I, D, and G are expected to bind the negatively charged viral genome (Fig.  4C ). In LY1, basic residues Arg61, Lys62, Arg64 (β−strand B), Arg133 (β−strand D), Lys197 (β−strand G), Lys535, and Lys541 (β−strand I) are all oriented toward the particle interior and are likely responsible, together with the ARM motif, for binding the negatively charged viral genome. Alignment of 2,201 Betatorquevirus ORF1 sequences reveals that several of these putative DNA-binding residues are conserved within ORF1 (deposited https://doi.org/10.5281/zenodo.11099132 ). In addition, alignment of LY1 and LY2 with 3 randomly selected Betatorquevirus ORF1s, 5 Gammatorquevirus ORF1s, and 5 Alphatorquevirus ORF1s shows >50% conservation of these basic residues, suggesting a general conservation for basic residues lining the interior of the ANV jelly roll particles (Supplementary Fig.  6 ).

A The LY1 core structure comprised of 60 jelly roll domains pack in icosahedral symmetry with one domain uniquely colored in red. B Two jelly roll domains are shown in red with the observed C-terminal domain backbone colored in cyan. The jelly roll domains are arbitrarily labeled JR 1 and JR 2 with the first (K48) and last (V562) observed residues for each protomer labeled with the corresponding number for clarity. C A single jelly roll domain is oriented to show the β-sheet on the interior of the particle core. Sidechains of basic residues in position to contact with the viral genome are shown and labeled. D The sequence of LY1 is shown for the ARM, jelly roll, and C-terminal regions. Basic residues of LY1 positioned to potentially contact the viral genome and shown in ( C ) are indicated with asterisks above the sequence along with the secondary structure elements.

Spike domain

Residues 229–530 of the H-I insertion form the spike domain that extends ~6 nm from the jelly roll core (Fig.  5 ). A β-strand (residues 245–250) extending from jelly roll β−strand H is the first component of the spike P1 subdomain and is N-terminal to the spike P2 subdomain (residues 251–385). The previously described HVR of ORF1 comprises the majority of the spike P2 subdomain. The remaining residues of the spike domain (residues 386–530) form five additional β−strands and eight helices, which, together with the residues 245–250 strand, fold into the spike P1 subdomain. The local resolution of P1 is only slightly lower than the jelly roll domain (~3–3.5 Å), while the resolution of P2 is within 4.5–5 Å. This is likely a consequence of both being further from the radius of gyration and some flexibility of the HVR residues.

A The anellovirus particle structure shown as a surface rendering. The particle is shown in gray with 5 spikes forming a crown structure, numbered for clarity, shown in surface rendering and colored as in Fig.  1 . B The exterior of the crown structure shown from the side. Five spike domains are colored as in ( A ). C The spike domain sequence of LY1 is shown colored as in Fig.  1 . Secondary structure elements are indicated above the sequence.

Neighboring spike domains pack together around the five-fold symmetry axis to form a ringed structure of five spike domains henceforth called the crown (Fig.  5A , B). A receptor for anelloviruses has not been identified to date. Given the diverse tissue distribution of different anellovirus strains, it is possible that residues of the spike P2 HVR, which are the most surface-exposed, are involved in viral attachment and infection. However, we hypothesize that the hypervariable sequence of the spike P2 subdomain serves to aid in immune evasion rather than harboring a receptor-binding motif. If this is the case, a receptor-binding motif on the better-conserved spike P1 subdomain, or even on the surface of the jelly roll core, may exist.

Absent from both structures is the majority of the C-terminal region (in the higher resolution ΔC-Term particle, we see density up to residue 562). The portion of the C-terminus we can see extends from the jelly roll core toward the surface near the 2-fold axis between the spike domain crowns. This suggests that if the C-terminus is present on the full-length viral particle, it could be located on the particle surface. The C-terminal region is predicted to be helical in nature and has a glutamate-rich region (residues 636–640 on LY1) which could N-terminally cap a conserved helical motif. To determine if the C-terminal residues of LY1 would form a helical structure, we performed circular dichroism experiments on the C-terminal peptide used to generate the aforementioned antibodies. Indeed, circular dichroism shows the C-terminal region (residues 635-672), which has a leucine-zipper motif, is helical in solution (Supplementary Fig.  7 ). Given that anelloviruses have evolved to evade the immune system, and that the C-terminal region is the immunodominant region of ORF1 and may be antagonistic to particle formation, we hypothesize that the helical C-terminal region is required for wild-type particle assembly but somehow modified or removed during the viral life cycle.

Despite anelloviruses constituting the majority of the human virome, their capsid structure was unknown and the identity of the capsid protein itself had not been established experimentally 19 . Our determination of the LY1 structure demonstrates that ORF1 encodes the capsid protein and that anelloviruses evolved a novel spike domain that extends around the 5-fold axis to form a crown structure. These crowns are capped with hypervariable P2 regions, which likely inhibits the development of antibodies against the better-conserved spike P1 subdomain or the jelly roll core via steric hindrance. Butkovic et al. recently performed structure prediction, corroborating the insertion of P1 and P2 subdomains between strands H and I and suggesting that anelloviruses evolved from an ancestral circovirus-like genome into which progressively larger insertions were selected 13 . Given that this insertion, containing the HVR, is the region of highest variation between the genera, it suggests the virus adopted a mechanism of constantly mutating its crown structure as a way of evading immunodetection.

ANVs have an elongated C-terminal region, predicted to be helical, which is absent in circoviruses such as BFDV. It is unclear what role the C-terminus would play in the viral life cycle, however given that the recombinant ORF1 protein does not seem to form symmetrical particles with the C-terminal region attached, one might hypothesize it has a role in particle assembly. It is possible that manipulation of the C-terminal region through post-translational modification or interaction with some other protein factor could stimulate particle formation in the wild-type virus. Alternatively, given that the portion of the C-terminal region in the ANV-like particles is on the viral surface, perhaps the C-terminus plays a role in viral attachment and infection. Whatever the function, it is interesting that ANVs evolved this novel feature unique from the structurally similar circovirus, and this could be an area of investigation for future studies.

While the spike domain and C-terminal regions represent novel features distinct from the otherwise structurally similar circoviruses, the ARM and jelly roll particle core are very similar. Sarker et al. have previously solved two different structural forms of BFDV capsid protein oligomers by X-ray crystallography (a decamer form and a 60 mer icosahedral capsid structure proposed to represent stages of viral maturation). When ssDNA primers were added to the decamer species, the capsomeres reassembled around the nucleic acid to form the particle, and density for ssDNA bound to basic residues of the jelly roll could be observed. ANV ORF1 particles share similar DNA contact residues with BFDV capsids (Fig.  6 ). The B-I-D-G β -sheet represents the majority of the particle interior surface for both viruses. In BFDV, residues of jelly roll strand F are also oriented toward the particle interior (Fig.  6A, B ). Two basic residues (Lys 154 and Lys 155) are oriented toward the neighboring jelly roll, placing these residues in close proximity to the basic residues of strand I and making contact with the DNA strand. While there is no visible DNA in the ANV particle structure (and no detectable DNA in the purified ORF1 sample) it is interesting to note that the jelly roll domains have a similar packing with regards to the jelly roll domains forming the pentamer (Fig.  6C, D ). On LY1, basic residue Lys 180 (just N-terminal to the jelly roll F-strand) is in close proximity to basic residue Lys 535 on the I-strand of the neighboring jelly roll domain. This orientation is very similar to what is observed in BFDV. In the alignment of 2201 Betatorquevirus ORF1s, residues aligned to L1 Lys 180, as well as residues aligned to LY1 residues 178 and 179 are predominantly basic (lysine, arginine or histidine) or polar (asparagine; Supplementary Fig.  8 ). This observation would suggest that this region could have conserved DNA-binding residues similar to the Lys residues 154 and 155 of BFDV.

A A BFDV capsid pentamer showing basic residues (magenta and blue) interacting with single-stranded DNA (green). Two jelly roll domains are labeled JR1 and JR2. B An expanded view of ( A ) showing the F-strand of the jelly roll (blue) orients towards the neighboring jelly roll in such a way as to place basic residues Lys 154 and 155 in contact with the DNA strand bound by basic residues on the neighboring jelly roll (magenta). C An ORF1 pentamer from the ANV-like particle shown in the same orientation and labeled as the BFDV pentamer in ( A ). D Expanded view of ( C ) showing that residue Lys 180 and the F-strand (labeled and shown in blue) are in close proximity to basic residues of the neighboring jelly roll (magenta).

The structure of the LY1 Betatorquevirus can now be used to guide future anellovirus research, exploring the semi-conserved regions of the P1 subdomain, C-terminal region, and jelly roll regions and how they relate to receptor binding and tissue tropism. In addition, further exploration of whether and how the ARM motif and basic residues of the jelly roll interior coordinate DNA binding and whether, like BFDV, ANVs undergo a decamer intermediate state prior to genomic encapsidation would be of interest. Lastly, the diversity and immunological stealth of anelloviruses invite continued study of their application to the delivery of therapeutic genes because they could potentially target cell types not currently addressed by existing vectors and they may be less susceptible to pre-existing immunity or to the development of neutralizing antibodies following initial treatment 20 . In sum, the availability of the present structures will shed light on the anellovirus lifecycle and will help guide the study of anellovirus-based gene therapy vectors.

Antibody generation and western blot analysis

LY1 antibodies were generated by immunizing rabbits with a synthetic peptide representing one of three portions of the ORF1 protein (jelly roll residues 46–58: KIKRLNIVEWQPK, spike domain residues 485–502: SPSDTHEPDEEDQNRWYP, C-terminal domain residues 635–672: SEEEEESNLFERLLRQRTKQLQLKRRIIQTLKDLQKLE) conjugated to a carrier protein by an engineered N-terminal cysteine. Antibodies used in this study were generated and purified by Thermo Fisher Scientific, Custom Antibody Production. Briefly, rabbits were immunized twice with the indicated peptide and polyclonal antibodies were purified using protein A purification. Western blot analysis was performed using NUPAGE 4–12% gels (ThermoFisher) transferred to nitrocellulose membranes using the Transblot Turbo system (BioRad). Membranes were blocked with blocking buffer (LI-COR Biosciences), probed ~16 h with primary antibody at 1 ug/ml, and detected using anti-rabbit IRDye infrared secondary antibodies and imaging system (LI-COR Biosciences).

Sf9 construct design, cell culture, and protein expression/purification

The LY1 ORF2 and ORF1 sequences 14 were codon-optimized for insect cells and with different ORF1 construct length variations (full-length ORF1, ΔARM with a deletion of residue 2–45, and ΔARM/ΔC-term with deletions of residues 2–45 and 552–672). The ORF2 and ORF1 constructs were cloned into the pFastBac Dual plasmid which was used to generate baculoviruses (Genscript and Medigen). To express the ORF1 proteins, Sf9 cells (Gibco 11496015) were infected by baculovirus with multiplicity of infection = 1 and the cells were cultured for 3 days at 27 °C and harvested by centrifugation. The cells were suspended in 50 mM Tris pH 8.0, 100 mM NaCl, and 2 mM MgCl 2 . Cells were lysed by treatment with 0.01% Triton X-100 (Sigma-Aldrich 11332481001) and subjected to micro-fluidization, and treated with protease inhibitors (Thermo Scientific Halt Protease Inhibitor Cocktail PI78438) and DNase (Benzonase®; Sigma). Cell lysate was subsequently purified using HiTrap Heparin affinity chromatography (Cytiva) followed by size-exclusion chromatography (HiPrep 16/60 Sephacryl S-500 HR; Cytiva).

Mammalian expression construct design, cell culture, and protein expression/purification

The LY1 ORF1 ΔC-term (1-609) sequences were codon-optimized for mammalian cells and cloned into a cytomegalovirus promoter-expression vector (Genscript). The ORF1 construct was transiently transfected into Expi293 cells (Gibco A14527) using polyethylenimine and the cells were harvested after 3 days at 37 °C by centrifugation. The cells were suspended in 50 mM Tris pH 8.0, 100 mM NaCl, and 2 mM MgCl 2 . Cells were lysed by micro-fluidization and treatment with 0.01% Triton X-100 (Sigma-Aldrich 11332481001) and then treated with protease inhibitors (Thermo Scientific Halt Protease Inhibitor Cocktail PI78438) and DNase (Benzonase®; Sigma). Cell lysate was subsequently purified using HiTrap Heparin chromatography (Cytiva) followed by size-exclusion chromatography (HiPrep 16/60 Sephacryl S-500 HR; Cytiva). Lysate from 5 liters of expression was loaded onto a 100 mL Heparin Sepharose column (Cytiva, #90100192) and eluted with a 2 molar NaCl gradient followed by Capto Core 400 purification. Protein and DNA amounts were quantified by using bicinchoninic acid (Pierce™ BCA Protein Assay Kit Cat# 23225) and picogreen (Quant-iT™ PicoGreen™; the limit of detection is 50 pg) methods, respectively. Protein amount was quantified as 0.7 mg/mL and DNA estimation results show no detectable amounts of DNA. Concentrated protein quality was evaluated by SDS-PAGE, western blot, and negative stain EM.

Negative-stained EM data collection and analysis

Purified ΔARM with a deletion of residue 2–45, and ΔARM/ΔC-term (expressed in SF9 cells) and LY1 ORF1 ΔC-term (1–609) (expressed in mammalian cells) were adsorbed for ∼ 1 min to parlodion carbon-coated 400 mesh copper grids which were rendered hydrophilic by glow discharge at low pressure in air. Grids were washed with three drops of MilliQ water and stained with two drops of 1% uranyl acetate. Electron micrographs were recorded at a nominal magnification of 30,000 × with JEOL 1400 Flash Transmission Electron Microscope operated at 100 kV and equipped with a NanoSprint43L-MkII CMOS camera system with 43 megapixel sampling region (Advanced Microscopy Techniques).

Cryo-EM data collection and data analysis and molecular refinement

For the SF9 fragment virus-like particle structure determination, 3 µl of 0.3 mg/ml particle sample was applied to a 1.2 × 1.3 graphene oxide grid. A total of 11,083 micrographs were collected from Glacios cryo-TEM (Thermo Fisher Scientific) operated at 200 kV with a Falcon 4 direct electron detector, at a nominal defocus range of -1.0 – -2.5 μm and accumulated dose of 19.59 e - /Å for a total of 15 frames in 3 min. The pixel size was 0.923 Å, with the magnification 150,000×. Automated data collection was carried out by Leginon 21 software. All micrographs were motion-corrected by Relion-4.0 22 , 23 , 24 implemented MotionCor2 24 , and the contrast transfer function (CTF) parameters were estimated by Gctf 25 . With manually picked particles from 20 micrographs to train the network, SPHIRE-crYOLO 26 automatically picked 58,391 particles along with the PhosaurusNet network. All particles were extracted by Relion-4.0 and rescaled to 2-fold (pixel size 1.846 Å), followed by subsequence 2D classifications with 350 Å mask diameter to remove any junk particles. Two iterations of 2D classifications resulted in 11,185 particles, which were merged and reextracted to generate a de novo 3D initial model by Relion-4.0 with I1 symmetry. Notably, several similar initial models without symmetry imposed were obtained by Relion-4.0 and cisTEM 27 . To obtain a better classification result, all particles were first subjected to a Refine3D with initial angular sampling 3.7° and local angular search 0.9° per step. The alignment parameters of each particle were transferred to a 3D classification with angular sampling interval 0.9° and local angular search 5° per step. 3D classification of the entire particle set attributed most of the particles into a single class. After CtfRefine and Bayesian polishing, the post-processing results in 3.98 Å resolution under the gold standard (with Fourier shell correlation [FSC] = 0.143).

For the expi293 cell-derived ΔC-term structure determination 3 grid preparation, 3 µl of 0.4 mg/ml particle sample was applied to a 1.2 × 1.3 CCFQ and performed triple blot. A total of 19,697 micrographs were collected using a Titan Krios microscope operated at 300 kV equipped with a Gatan quantum 967 LS imaging filter, a Volta phase plate, and Gatan K3 Direct Detection camera at a nominal defocus range of –1.0 to 1.6 μm 50.53 e-/Å2 for a total of 35 frames. The pixel size was 0.834 Å with magnification 105,000×.

All micrographs were motion-corrected by Motioncor2 24 implemented in Cryosparc 3.3 and the CTF parameters were estimated by Gctf 25 . Particles were picked by automated particle picking using Cryosparc 3.3, and particles were extracted and subjected to 2D classification in Cryosparc 3.3. After iterative rounds of reference free 2D classification, 23,193 particles were selected for further processing. Initial models were generated ab initio from all selected particles, and the best 3D class was submitted to homogeneous 3D refinement that included dynamic masking. The final reported resolution of 2.69 Å was based on the gold standard FSC = 0.143 criterion 28 . Maps were visualized using Chimera 29 .

The initial anellovirus TTMV-LY1 capsid monomer structure was predicted by TrRosetta 30 , and the TTMV-LY1 structure was further predicted by RosettaCM 31 . Structural refinement was performed by Rosetta 32 and Phenix 33 , and fine-adjusted by COOT 34 , 35 . Detailed definitions of the secondary elements are shown in Supplementary Fig.  9 .

Statistics and reproducibility

The negative staining-electron micrograph of full-length LY1 ORF1 (Fig.  1D ) is a representative of over 100 micrographs recorded while the micrograph of C-terminal truncation (Δ552-672) is a representative of over 20 micrographs recorded (Fig.  1F ). The negative staining-electron micrographs of LY1 fragment (Fig.  1E and Supplementary Fig.  1A ) are representative micrographs of over 100 micrographs taken, while the cryo-electron micrograph of the same sample (Supplementary Fig.  1B ) is a representative of 11,000 micrographs recorded. The negative staining-electron micrograph of ΔC-term LY1 ORF1 (Fig.  2E ) is a representative of over 100 micrographs recorded while the cyro-electron micrograph of the same sample (Supplementary Fig.  3A ) is a representative of 19,000 micrographs recorded.

Sequence alignments

Sequences of select anellovirus ORF1s were taken from Genbank and aligned using Clustal Omega 36 . ORF1 sequences used in the alignment are indicated by their accession numbers in Supplementary Fig.  6 . ORF1 sequences used in the alignment are indicated by their accession numbers in Supplementary Fig.  4 , except LY1 (YP 007518450.1), LY2 (AGG91484.1), SAfiA (MN779270.1), and JA20 (AF122914.3). The consensus sequence is a composite of residues conserved (identity >30%) or similar residues (>50%; hydrophobic and bulky aromatic (Tyr) residues indicated with ϕ, and hydrophilic and basic residues indicated with γ). Alignments and consensus sequences are shown in Fig.  2 and Fig.  3 and Supplementary Fig.  4 . Bulk alignments of 2201 Betatorquevirus ORF1 amino acid sequences retrieved from Genbank were aligned utilizing MAFFT 36 with the ‘auto’ parameter. Amino acid frequencies were computed for each position of the alignment and visualized in Supplementary Fig.  6 .

Circular dichroism

To determine the secondary structure of the LY1 C-terminal domain, the peptides used to generate C-terminal antibodies (residues 635–672) were analyzed at 25.8 µM in PBS on a Jasco J-815 circular dichroism spectropolarimeter using 2 mm path length cell at ambient temperature. Each data set is the average of three consecutive scans. The secondary structure was determined by using the CDPro software package 37 , compared with a reference set containing 56 proteins (IBasis = 10). The final secondary structure fractions were averaged over the results from three programs (SELCON3, CDSSTR, CONTINLL) in CDPro.

Reporting summary

Further information on research design is available in the  Nature Portfolio Reporting Summary linked to this article.

Data availability

The Betatorquevirus amino acid alignment referenced in this paper is available via Zenodo utilizing DOI number https://doi.org/10.5281/zenodo.11099132 ( https://doi.org/10.5281/zenodo.11099133 ). Both particle structures from SF9- and Expi 293-expressed ORF1 were deposited to the Protein Data Bank and EMDB, and the accession codes are PDB: 8CYG, EMDB: 27077 and PDB: 8V7X, EMDB:43009 respectively.  Source data are provided with this paper.

Shulman, L. M. & Davidson, I. D. Viruses with circular single-stranded DNA genomes are everywhere! Ann. Rev. Virol. 4 , 1–22 (2016).

Google Scholar  

Takahashi, K., Iwasa, Y., Hijikata, M. & Mishiro, S. Identification of a new human DNA virus (TTV-like mini virus, TLMV) intermediately related to TT virus and chicken anemia virus. Arch. Virol. 145 , 979–993 (2000).

Article   CAS   PubMed   Google Scholar  

Kaczorowska, J. & Hoek, Lvander Human anelloviruses: diverse, omnipresent and commensal members of the virome. FEMS Microbiol. Rev. 44 , 305–313 (2020).

Article   CAS   PubMed   PubMed Central   Google Scholar  

Nishizawa, T. et al. Quasispecies of TT virus (TTV) with sequence divergence in hypervariable regions of the capsid protein in chronic TTV infection. J. Virol. 73 , 9604–9608 (1999).

Kaczorowska, J. et al. Diversity and long-term dynamics of human blood anelloviruses. J. Virol. 96 , e00109–22 (2022).

Koonin, E. V., Dolja, V. V. & Krupovic, M. The healthy human virome: from virus–host symbiosis to disease. Curr. Opin. Virol. 47 , 86–94 (2021).

Stulberg, E. et al. An assessment of US microbiome research. Nat. Microbiol. 1 , 15015 (2016).

Bal, A. et al. Metagenomic next-generation sequencing reveals individual composition and dynamics of anelloviruses during autologous stem cell transplant recipient management. Viruses 10 , 633 (2018).

Kakkola, L. et al. Expression of all six human Torque teno virus (TTV) proteins in bacteria and in insect cells, and analysis of their IgG responses. Virology 382 , 182–189 (2008).

Requião, R. D. et al. Viruses with different genome types adopt a similar strategy to pack nucleic acids based on positively charged protein domains. Sci. Rep. 10 , 5470 (2020).

Article   ADS   PubMed   PubMed Central   Google Scholar  

Sarker, S. et al. Structural insights into the assembly and regulation of distinct viral capsid complexes. Nat. Commun. 7 , 13014 (2016).

Article   ADS   CAS   PubMed   PubMed Central   Google Scholar  

Nath, B. K. et al. Structural perspectives of beak and feather disease virus and porcine circovirus proteins. Viral Immunol. 34 , 49–59 (2021).

Butkovic, A. et al. Evolution of anelloviruses from a circovirus-like ancestor through gradual augmentation of the jelly-roll capsid protein. Virus Evol. 9 , vead035 (2023).

Article   PubMed   PubMed Central   Google Scholar  

Galmès, J. et al. Potential implication of new torque teno mini viruses in parapneumonic empyema in children. Eur. Respir. J. 42 , 470–479 (2013).

Article   PubMed   Google Scholar  

Itoh, Y. et al. Visualization of TT vrus particles recovered from the sera and feces of infected humans. Biochem. Biophys. Res. Commun. 279 , 718–724 (2000).

Li, T.-C. et al. Essential elements of the capsid protein for self-assembly into empty virus-like particles of hepatitis E virus. J. Virol. 79 , 12999–13006 (2005).

Xie, Q. et al. The atomic structure of adeno-associated virus (AAV-2), a vector for human gene therapy. Proc. Natl Acad. Sci. USA 99 , 10405–10410 (2002).

Tsao, J. et al. The three-dimensional structure of canine parvovirus and its functional implications. Science 251 , 1456–1464 (1991).

Article   ADS   CAS   PubMed   Google Scholar  

Desingu, P. A., Nagarajan, K. & Dhama, K. Can a Torque teno virus (TTV) be a naked DNA particle without a virion structure? Front. Virol . 1–5 https://doi.org/10.3389/fviro.2022.821298 (2022).

Arze, C. A. et al. Global genome analysis reveals a vast and dynamic anellovirus landscape within the human virome. Cell Host Microbe. 29 , 1305–1315.e6 (2021).

Suloway, C. et al. Automated molecular microscopy: the new Leginon. Syst. J. Struct. Biol. 151 , 41–60 (2005).

Article   CAS   Google Scholar  

Zivanov, J., Nakane, T. & Scheres, S. H. W. Estimation of high-order aberrations and anisotropic magnification from cryo-EM data sets in RELION-3.1. IUCrJ. 7 , 253–267 (2020).

Kimanius, D., Dong, L., Sharov, G., Nakane, T. & Scheres, S. H. W. New tools for automated cryo-EM single-particle analysis in RELION-4.0. Biochem. J. 478 , 4169–4185 (2021).

Zheng, S. Q. et al. MotionCor2: anisotropic correction of beam-induced motion for improved cryo-electron microscopy. Nat. Methods 14 , 331–332 (2017).

Zhang, K. Gctf: Real-time CTF determination and correction. J. Struct. Biol. 193 , 1–12 (2016).

Wagner, T. et al. SPHIRE-crYOLO is a fast and accurate fully automated particle picker for cryo-EM. Commun . Biol. 2 , 218 (2019).

Grant, T., Rohou, A. & Grigorieff, N. cisTEM, user-friendly software for single-particle image processing. eLife 7 , e35383 (2018).

Scheres, S. H. W. & Chen, S. Prevention of overfitting in cryo-EM structure determination. Nat. Methods 9 , 853–854 (2012).

Goddard, T. D., Huang, C. C. & Ferrin, T. E. Visualizing density maps with UCSF. Chimera. J. Struct. Biol. 157 , 281–287 (2007).

Yang, J. et al. Improved protein structure prediction using predicted interresidue orientations. Proc. Natl Acad. Sci. USA 117 , 1496–1503 (2020).

Song, Y. et al. High-resolution comparative modeling with RosettaCM. Struct. (Lond., Engl. 1993) 21 , 1735–1742 (2013).

Wang, R. Y.-R. et al. Automated structure refinement of macromolecular assemblies from cryo-EM maps using Rosetta. eLife 5 , e17219 (2016).

Liebschner, D. et al. Macromolecular structure determination using X‐rays, neutrons and electrons: recent developments in Phenix. Acta Crystallogr. Sect. D. 75 , 861–877 (2019).

Article   ADS   CAS   Google Scholar  

Emsley, P., Lohkamp, B., Scott, W. G. & Cowtan, K. Features and development of coot. Acta Crystallogr. Sect. D. 66 , 486–501 (2010).

Sievers, F. et al. Fast, scalable generation of high‐quality protein multiple sequence alignments using Clustal Omega. Mol. Syst. Biol. 7 , 539–539 (2011).

Katoh, K. & Standley, D. M. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol. 30 , 772–780 (2013).

Sreerama, N. & Woody, R. W. Estimation of protein secondary structure from circular dichroism spectra: comparison of CONTIN, SELCON, and CDSSTR methods with an expanded reference set. Anal. Biochem. 287 , 252–260 (2000).

Download references

Acknowledgements

We thank Dr. Joseph Che-Yen Wang for manuscript discussions and suggestions, and Maria Ericsson and the team at Harvard Medical School Electron Microscopy facility for their technical support. Cryo-EM data were collected at NanoImaging Services under the leadership of Dr. Giovanna Scapin and Dr. Phat Dip. Expi 293 produced virus-like particle cryo-EM data collection and image processing at NanoImaging Services were supported by Dr. Weili Zheng. Baculovirus was generated by Dr. Peter Pushko at Medigen. We thank Peter Riebling for proof reading and suggestions.

Author information

Shu-hao Liou

Present address: Carbon Biosciences, Waltham, MA, 02451, USA

Noah R. Cohen

Present address: AbbVie Bioresearch Center, Worcester, MA, 01605, USA

Saadman Islam

Present address: GSK, Cambridge, MA, 02139, USA

Harish Swaminathan

Present address: DaCapo Brainscience, Cambridge, MA, 02139, USA

Nathan Yozwiak & Roger J. Hajjar

Present address: Gene and Cell Therapy Institute, Mass General Brigham, Cambridge, MA, 02139, USA

Simon Delagrave

Present address: Delagrave Life Sciences, LLC, Sudbury, MA, 01776, USA

These authors contributed equally: Shu-hao Liou, Rajendra Boggavarapu.

Authors and Affiliations

Ring Therapeutics, 140 First Street Suite 300, Cambridge, MA, 02139, USA

Shu-hao Liou, Rajendra Boggavarapu, Noah R. Cohen, Yue Zhang, Ishwari Sharma, Lynn Zeheb, Nidhi Mukund Acharekar, Hillary D. Rodgers, Saadman Islam, Jared Pitts, Cesar Arze, Harish Swaminathan, Nathan Yozwiak, Tuyen Ong, Roger J. Hajjar, Yong Chang, Kurt A. Swanson & Simon Delagrave

You can also search for this author in PubMed   Google Scholar

Contributions

S.-H.L. and R.B. purified samples, conducted negative EM, and processed the data with structural refinement. N.C. aided in structural model refinement. S.-H.L., J.P., and N.C. designed Sf9 ORF1 constructs and R.B. and Y.Z. designed HEK ORF1 constructs. Y.Z. and N.C. performed circular dichroism measurement. C.A. and H.S. performed multiple sequence alignments and analysis. I.S., N.M.A., H.D.R., S.I., and L.Z. performed ORF1 purification and construct screening. S.D., R.H., N.Y., T.O., and Y.C. oversaw the project. K.S. aided in ORF1 construct design and purification, aided in the design of the experiments and oversaw the structural determination. S.-H. L., R.B., S.D., and K.S. wrote the manuscript.

Corresponding author

Correspondence to Kurt A. Swanson .

Ethics declarations

Competing interests.

All authors are or were employed by and hold equity interests in Ring Therapeutics. S.-H. L., N.C., R.B., Y.Z., S.D., and K.S. are the inventors on a patent application related to this work.

Peer review

Peer review information.

Nature Communications thanks Mart Krupovic, and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Supplementary information, peer review file, source data, rights and permissions.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/ .

Reprints and permissions

About this article

Cite this article.

Liou, Sh., Boggavarapu, R., Cohen, N.R. et al. Structure of anellovirus-like particles reveal a mechanism for immune evasion. Nat Commun 15 , 7219 (2024). https://doi.org/10.1038/s41467-024-51064-8

Download citation

Received : 18 December 2023

Accepted : 26 July 2024

Published : 22 August 2024

DOI : https://doi.org/10.1038/s41467-024-51064-8

Share this article

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

By submitting a comment you agree to abide by our Terms and Community Guidelines . If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.

Quick links

• Explore articles by subject
• Guide to authors
• Editorial policies

Sign up for the Nature Briefing: Translational Research newsletter — top stories in biotechnology, drug discovery and pharma.

Experimental and Numerical Study of Taylor Bubble in Counter-Current Turbulent Flow

• RESEARCH ARTICLE - Special Issue - Challenges and Recent Advancements in Nuclear Energy Systems
• Open access
• Published: 21 August 2024

Cite this article

You have full access to this open access article

• Iztok Tiselj   ORCID: orcid.org/0000-0001-9340-5397 1 , 2 ,
• Jan Kren 1 , 2 ,
• Blaž Mikuž 1 ,
• Raksmy Nop 3 ,
• Alan Burlot 3 &
• Grégoire Hamrit 3  

The stagnant Taylor bubble in vertical isothermal turbulent counter-current flow was analyzed using 2D shadowgraphy experiments and two distinct high-fidelity numerical simulations. One simulation employed the geometrical VOF interface tracking method within the OpenFOAM code, while the other utilized the explicit front tracking method of the TrioCFD code. Interface recognition algorithms were applied to the photographs and compared with the results of 3D simulations performed with LES and pseudo-DNS accuracy in OpenFOAM and TrioCFD, respectively. The measured Taylor bubbles exhibited an asymmetric bullet-train shape and a specific speed, which were compared with the predictions of both numerical approaches. Reproducing the experiment proved challenging for both otherwise well-established methods frequently used in interface tracking simulations of two-phase flows. Grid resolution and subgrid turbulent models, known for their success in single-phase turbulence, were less accurate near the water–air interface. Additional experimental parameters compared with simulations were related to the dynamics of tiny disturbance waves with amplitudes ranging from 10 to 100 µm along the interface of the Taylor bubbles. The speed and spectra of the surface disturbance waves were reproduced numerically with moderate success despite detailed grid refinement in the relevant region of the computational domain.

Avoid common mistakes on your manuscript.

1 Introduction

Gas–liquid mixture flows exhibit various two-phase flow patterns, with vertical pipes commonly experiencing bubbly, slug, churn, annular, and droplet flow regimes [ 1 ]. The specific flow regime depends on factors such as flow velocities, phase volume fractions, fluid properties, pipe size, and orientation. This study is focused on Taylor bubble flow, which falls under the slug flow regime. Taylor bubbles are bullet-shaped bubbles that move at different speeds from the bulk liquid, occupying almost the entire pipe cross section. Slug flows are relevant for chemical, nuclear, petroleum, and other types of processing engineering. They are encountered in a wide range of practical applications, including vaporizers, boilers, filtration and membrane processes [ 2 ], as well as extreme events in the petroleum industry [ 3 ] or steam generators in nuclear power plants. The most recent review paper discussing vertical gas–liquid slug flows by Holagh & Ahmed from 2024 [ 4 ] demonstrates the extensive scope of this research area. This remarkable review is citing 470 references.

The type of slug flows relevant to the present study occurs in the inertia-dominant regime, where the influence of viscosity and surface tension is minimal [ 1 ]. In this regime, the drift velocity of Taylor bubbles \({U}_{0}\) in pipe with diameter D is given by the correlation \({U}_{0}=k\sqrt{gD}\) ( g acceleration of gravity). Based on the constant value of k ≈0.35, this correlation predicts a drift velocity of approximately 0.18 m/s for Taylor bubbles in our experiments [ 5 ]. This value is close to the average measured liquid velocity, \({U}_{L}\) , which is in the downward direction (negative sign) and keeps the bubble fixed in position. Dumitrescu [ 6 ] demonstrated that soon after the liquid flow was directed downward, the Taylor bubble became unstable. One of the earliest detailed experiments on counter-current turbulent flow was performed by Martin [ 7 ], who studied air–water mixtures in circular pipes with diameters of 2.6, 10.16, and 14.0 cm. Martin found that the bubble velocity in counter-current slug flow could not be adequately explained by existing theories for co-current background flow or stagnant liquid conditions. This discrepancy is due to bubble instability, which increases bubble velocity when the bubble is displaced from the pipe axis. Lu and Prosperetti [ 8 ] performed a stability analysis and showed that the breakup of Taylor bubble symmetry occurs at liquid velocities below a critical negative velocity of \({U}_{c}=-0.13\sqrt{gD}\) . Figueroa-Espinoza and Fabre [ 9 ] performed numerical analyses of symmetry breakup at different surface tension values. They found that asymmetry leads to increased bubble velocity and a decreased curvature radius at the stagnation point of the bubble nose. Fabre and Figueroa-Espinoza [ 10 ] further investigated symmetry breakup experimentally and determined that asymmetry is largely independent of whether the flow regime is turbulent or laminar. They identified the vorticity-to-radius ratio at the stagnation point as a crucial parameter for symmetry breakup. Fershtman et al. [ 11 ] also studied counter-current slug flow, measuring a liquid velocity that exactly balances buoyancy \({U}_{L}=0.35\sqrt{gD}\) = 0.178 m/s, which was also observed in experiments in the present study. The latest study by Abubakar and Matar [ 12 ] provided a detailed numerical and parametric analysis of the effects of downward liquid velocity, viscosity, and surface tension on bubble shape and motion. Their linear stability analysis identified regions of dimensionless parameters where the bubble is unstable and assumes an asymmetric shape, and they explained the mechanisms behind symmetry breakup. Our experiments were conducted in the unstable region with asymmetric bubble shapes, necessitating dynamic flow rate control, as described in the following section.

The interactions between Taylor bubbles and turbulent liquid flow have been the subject of various studies. Unlike in laminar liquid flow, the tail of the bubble starts to break in the turbulent background flow. The breakup and recoalescence processes in the bubble wake region have been observed. The studies [ 13 , 14 ] measured the gas loss from a stationary Taylor bubble in a counter-current liquid flow using a special spherical Teflon cap to hold the bubble in a fixed position. More recent experiments involving the turbulent counter-current regime have utilized high-speed cameras in visible light to measure the bubble’s disintegration rate [ 15 ]. They have shown that disintegration by the breakup stops when the Taylor bubble is sufficiently short. This result is important for the present study focused on the Taylor bubbles, which do not lose their mass due to the breakup of the bubble’s tail. As shown in [ 15 , 37 ], even shorter bubbles are gradually losing their mass due to the dissolving of the gas in the liquid. However, this mechanism is much slower than the physical breakup and can be neglected over the time intervals relevant for the present study. Consequently, dynamical control of the liquid flow rate is used to trap the bubble in an equilibrium position for hours, and allowing for studies over several minutes.

Two fundamentally different approaches are used for numerical simulations of turbulent two-phase flows:

One approach is the Euler-Euler method, wherein the Navier–Stokes equations are solved independently for each phase [ 16 ].

Another approach is the one-fluid formulation, which involves applying a single set of governing equations across the entire domain, encompassing the interface [ 17 ]. Various methods for interface advection exist, with the Volume of Fluid (VOF) method [ 19 , 20 ], Front tracking method [ 21 ], and Level-Set method [ 30 ] being among the most widely used. The present research is focused on the (VOF) method and on the Front tracking method, which were tested with the stagnant Taylor bubble in turbulent flow.

The specific version of the VOF method used in this study is implemented in the open-source OpenFOAM computer code [ 18 ]. One of the main advantages of the VOF method over alternative methods is its well-established framework and guaranteed volume conservation. Central to the VOF method is the concept of a marker function, which represents the volume fraction of one fluid within each computational cell of the domain. A key challenge in advecting a marker function is the numerical diffusion that arises from using a cell-averaged marker function [ 21 ]. To mitigate this diffusion, the VOF method reconstructs the interface so that the marker does not move into a new cell until the current cell is completely filled. Reconstruction models of the VOF method are categorized into algebraic and geometric types. Significant effort has been focused on geometric methods because they produce better results than algebraic reconstruction methods [ 23 , 24 ]. One such geometric method is the piecewise linear interface calculation (PLIC) method, which has been investigated for large eddy simulations (LES) by Kren et al. [ 22 ] and is tested in the present paper. Large eddy simulation (LES) is a very useful compromise between the high accuracy and cost of direct numerical simulation (DNS) and the lower cost but reduced accuracy of Reynolds-averaged Navier–Stokes (RANS) simulations. The development of LES methods for multiphase flows is still in its early stages, largely due to the lack of experimental and DNS benchmarks. Klein et al. [ 25 ] have laid out a framework for developing LES in multiphase flows. In single-phase flows, modeling small scales is required only for the convective term, typically achieved by adding eddy viscosity to the equations. However, multiphase flows require modeling several terms, with at least two being particularly significant. While the convective term can be approached in a similar manner to single-phase flows [ 26 ], new models are needed for other closures. Recent developments highlight that the closure of the sub-grid term for surface tension is the most critical [ 22 ].

Numerical simulations of multiphase flows using the front-tracking method have been presented by Unverdi in 1992 [ 27 ] and upgraded by Tryggvason et al. [ 28 ]. Like in the VOF method, interfacial terms are incorporated by adding the appropriate sources as delta functions at the phase boundaries. The unsteady Navier–Stokes equations are solved using a conventional finite volume or finite element method on a fixed, structured grid. The interface, or front, is explicitly tracked using connected marker points in Lagrangian coordinates, with surfaces in the 3D domain. This type of method is implemented in TrioCFD code [ 35 ], which was the second type of numerical model tested with the stagnant Taylor bubble in turbulent flow. Interfacial source terms such as surface tension are computed on the front and transferred to the fixed grid. The advection of fluid properties, like density, is handled by following the motion of the front. When large topology changes occur, the distance between marker points can increase or decrease, leading to reduced accuracy. To address this, new marker points might be added where marker density is low, and some marker points are removed where density is high. Although the method is complex to implement, it provides highly accurate tracking of the interface position. However, interface breakup does not occur well unless a special model is implemented.

Existing numerical simulations of Taylor bubbles have predominantly focused on stagnant or co-current background liquid flows across various setups, ranging from simplistic 2D and Euler-Euler simulations [ 31 , 32 ] to comprehensive 3D simulations with interface tracking [ 33 , 34 ]. LES investigations of Taylor bubbles within co-current turbulent regimes, without special sub-grid scale models for bubble coalescence or breakup, revealed disintegration rates one to two orders of magnitude faster than observed in experiments of Taylor bubbles within counter-current turbulent flow [ 22 , 36 ].

The present study aims to address details of the Taylor bubble’s interface dynamics in the region of the liquid film and the capabilities of state-of-the art numerical schemes to describe the flow of the thin liquid film and the velocity fields in the liquid and air phase around the interface. An important part of the current work is the study of interface waves: as discussed below, accurate measurements of these waves were performed despite their tiny amplitudes. Consequently, we have examined properties of the interface waves predicted by the OpenFOAM and TrioCFD codes. All instances of waves are traveling fluid oscillations sustained by the surface tension force. For experimental Taylor bubble the interface waves were studied by Kren et al. [ 24 ]. The work on the interfacial waves analysis spans beyond the films of Taylor bubbles: closely related interfacial waves are observed in vertical annular flows, where a recent review is available [ 39 ] and a more specific example by Tekavčič et. al. [ 38 ] who analyzed interfacial waves in water–air churn flow regime. Another slightly less relevant area are studies of horizontally stratified flows with the presence of capillary waves, where some of the recent activities can be found in [ 40 , 41 ].

The goal of the present study is a detailed analysis of one specific experimental case performed within the experimental campaign [ 37 ], where stagnant Taylor bubbles were observed in the counter-current turbulent flow of water. The selected experimental case is reproduced with two state-of-the-art interface tracking approaches capable to operate in turbulent flow regime. Comparison of all measured parameters is performed: the main integral parameter is mass flow rate of the water that is balancing the buoyancy; flow rate is closely related to the axially asymmetric shape of the bubble, which is influencing the bubble drag. The last test is a comparison of the measured and computed interfacial waves that are traveling over the body of the Taylor bubble.

2 Experimental Setup

Experiments were conducted in a loop, as depicted in Fig.  1 . The test section consisted of a 1.5 m long glass pipe with an internal diameter of D = 26 mm and wall thickness 2 mm. Water is entering the test section from the top. All experimental cases, including the one considered in the present paper, were performed in the turbulent flow regime of the liquid above the bubble, with a Reynolds number of approximately 5600. The straight section of the pipe above the bubble spanned around 40 pipe diameters, ensuring statistically uniform turbulence impinging on the bubble. To maintain a constant water temperature of 30 °C, a heat exchanger was utilized in the tank. The Taylor bubble was injected into the test section with a syringe through a small dedicated connection beneath the test section. The flow through the test section was regulated using a control valve, which adjusted the flow distribution between the main loop and a bypass loop.

Schematics of the test loop

Observations of the Taylor bubble were carried out using a high-speed camera (camera Phantom v 1212, objective: TOKINA AT-X PRO D 100 mm), set up at a distance of around 30 cm between the objective and the pipe axis, with a field of view covering a 14 cm (~ 5 diameters) section of the pipe. The pipe was immersed in a rectangular glass section filled with water to minimize optical distortion. The observed Taylor bubbles typically had lengths ranging from 1.5 to 5 pipe diameters. Measurements were performed over different time intervals of 8, 4, 2, and 1 min, with camera frequencies of 100, 200, 400, and 800 Hz, respectively. Around 90 mm long bubble filmed over two minutes interval at frequency 400 Hz was used for detailed analyses in the present study.

The camera’s useful resolution for the measurements was approximately 1280 × 240 pixels, corresponding to around 9 pixels per millimeter. The measurement of absolute liquid film thickness is achieved with a precision ranging between 0.5 and 1 pixel. However, this level of precision introduces relative errors exceeding 40% for very thin films below 3 pixels. The estimated optical distortion, due to the light refraction, results in a maximum enlargement of the liquid film thickness by up to 2%. This value is significantly lower than the uncertainty associated with interface reconstruction.

Each experimental run was performed in the following steps:

Preparation of the water loop, establishment of the water circulation with an expected mass flow rate.

Set-up of the camera, illumination, pressure, temperature and mass flow rate sensors.

Injection of the air bubble.

Fine-tuning of the water mass flow rate to bring the bubble into the camera’s view.

Start of measurements and active fine-tuning of the mass flow rate during the measurement.

Data (image) processing.

In the counter-current flow configuration, the instability of the Taylor bubble requires dynamic adjustments of the mass flow rates during the experiment to ensure that the bubble remains within the camera’s field of view. Minor corrections of the valve position are made every few seconds, leading to fluctuations in the bulk liquid velocity within the test section. The bulk liquid velocity is based on the readings from the Coriolis flow meter, which measures mass flow rate through the section at a frequency of 1 Hz. The variations in the mean velocity measurement due to manual mass flow rate corrections range between 3 and 10% of the bulk velocity across various experimental cases. In analysis of the results these minor changes in the bubble position were considered with a sliding coordinate system which was fixed to the tip of the Taylor bubble’s nose, while the dynamics of the Taylor bubble was practically unaffected by the mass flow rate changes. As shown in [ 37 ] the bubbles moved up and down with vertical velocities below 0.01 m/s, which were considerably lower than the upstream mean liquid velocity of 0.18 m/s and the velocities around 1 m/s observed on the liquid–air interface of the bubble. While the system of bubble position control may present challenges when comparing results with similar experiments or numerical simulations, it closely resembles the numerical technique employed in high-fidelity simulations of co-current Taylor bubble flow [ 45 , 46 ]. In these simulations, the Taylor bubble is modeled within a moving frame of reference to ensure the bubble remains inside the computational domain. The same approach was used in TrioCFD simulation in this paper, while the adaptive mass flow rate boundary condition was used in the OpenFOAM.

The processing of each recording involved analyzing a set of 50,000 photographs using a dedicated in-house software. The best scaling distance on the photographs turned out to be the outer diameter of the glass pipe, which is not affected by the optical distortion. This software utilized widely used libraries for tasks such as fitting two-dimensional surfaces and one-dimensional lines, performing Fourier transformations, and cross-correlating one-dimensional functions. The algorithms and techniques employed in the software were based on methods described in the Numerical Recipes book [ 47 ]. The main focus of the computer codes was image processing, specifically the extraction of the Taylor bubble surface from the images. Given the constraints of automated analysis on a large dataset, we used established image processing methods found in the open literature [ 48 ]. In order to handle the large number of photographs, manual corrections and artifact removal were limited. To address this, a robust procedure was developed that could identify potential failures in bubble interface reconstruction. A description of the algorithm can be found in [ 37 ], while the intermediate results of the particular step of silhouette reconstruction are shown in Fig.  2 :

a =  > b: conversion of image density matrix into gradient matrix.

b =  > c: identification of the bubble outer surface

c =  > d: sub-pixel interface position refinement.

a original image, b magnitude of the gradients field (step 1), c extracted bubble interface and pipe inner walls at pixel level (step 2), d) refinement of the interface position at subpixel level (step 3) with pixel grid in the background: “ + ”—pixel level interface, “x”—sub-pixel level interface. All units in pixels

Distinguishing between absolute and relative accuracy is crucial when considering interface recognition. The absolute uncertainty of the interface position on a single photograph ranges from half a pixel to one pixel. However, when analyzing a time series or spatial profiles of the interface, the relative uncertainty of the interface motion between neighboring pixels in space or time is reduced by a factor of approximately 5 to around ± 0.1 pixel. This improvement in relative accuracy allows very precise characterization of interface movements.

3 Numerical methods

Modeling of the Taylor bubble in the counter-current turbulent flow was performed with two highly accurate but fundamentally rather different approaches (Table  1 ):

The Front tracking method is implemented in TrioCFD code [ 35 ]. The interface position is tracked with a Lagrangian mesh which is advected by the flow. Such an approach allows very accurate interface recognition and dynamics as the interface has a zero thickness; however, it means more expensive numerical algorithm.

The geometric VOF method implemented in the finite volume code OpenFOAM [22] introduces high-order interface capturing scheme, which consists of two parts—interface advection using void fraction property \(\alpha\) and interface reconstruction with different submodels of the VOF method. Computing time of this approach turns out to be shorter, but slightly less accurate in modeling the interface dynamics and surface tension effects.

3.1 Front Tracking in TrioCFD

TRUST/TrioCFD is an open-source CFD code developed by the CEA (the French Atomic Energy and Alternative Energies Commission). Massively parallel, it can handle various physical situations: single or two-phase flow, chemistry, fluid–structure interaction… To model a two-phase flow, one of the possible approaches of TrioCFD is to use its front-tracking method. The latter uses an Euler–Lagrange approach coupled with a Volume-of-Fluid like method. It uses the one-fluid formulation for the fluid problem and an explicit interface tracking by considering:

The indicator function χ k equals 1 in the phase k and 0 otherwise,

The variable \({\phi }_{k}\) (velocity, pressure, density…) having its given value in the phase k .

One can define a one-fluid variable as \(\phi = {\sum }_{k}{\chi }_{k}{\phi }_{k}\) . By summing the continuity equation for incompressible flows, the Navier–Stokes equation for each phase and the Laplace pressure law at the interface, one can derive the equation system:

with \(\overrightarrow{u}\) the velocity field, ρ the density, μ the dynamic viscosity, \(\overrightarrow{g}\) the gravity field, \(\sigma\) the surface tension at the interface, \(\kappa\) the local curvature of the interface \({\delta }_{I}\) the indicator function of the interface and \(\overrightarrow{n}\) the vector normal to the interface. The system of Eq. ( 1 ) is the one-fluid problem solved for the fluid mixture in the front-tracking of TrioCFD. The choice was made to use a non-conservative discrete surface tension rather than a classic Continuum Surface Force model as the latter introduces parasitic currents. The non-conservative effect has been studied and is found negligible [ 29 ].

Regarding the interface (Fig.  3 ), the Lagrangian mesh is advected between two time steps by the velocity field with simple advection equation applied to all Lagrangian markers on the interface \(i=1,N\) :

where \({\overrightarrow{s}}_{\text{i}}\) denotes position of the Lagrangian marker i, and \(\overrightarrow{u}({\overrightarrow{s}}_{\text{i}})\) represents the velocity at the marker position. To ensure the mass conservation, the transport of the phase indicator function with a calculation of the volume of gas is performed to refine the position of the Lagrangian markers. The interface treatment ends by a smoothing and a remeshing to regularize the markers position.

Illustration of the Lagrangian mesh tracking the interface in TrioCFD

In this investigation, an explicit Euler scheme was used for the time integration. Regarding the space discretization, we used the Finite Element Volume method, which is a hybrid between the finite volume and the finite element methods.

As illustrated in Fig.  4 the computational domain is a vertical circular pipe with a length of L = 23 cm and a diameter corresponding to the experimental pipe. At the outlet, a free pressure condition was implemented and at the pipe wall, a no-slip condition was imposed. Fixed mass flow rate was prescribed at the inlet ( v 0  = 0.17 m/s) and the moving frame of reference approach was used in TrioCFD simulation.

Left: scheme of the computational model. Right: cross section of the mesh

As the inlet flow is turbulent, a synthetic turbulence developing an isotropic and homogeneous turbulence was chosen as a boundary condition. The values of turbulent kinetic energy \(k\) and turbulent dissipation rate \(\varepsilon\) on inlet were computed with a preliminary RANS standard k-ε computation. Even though this does not exactly mimic the conditions in a pipe, this approach combined with a certain development length above the bubble head provides satisfactory results and a reasonable computational cost. Let us note that this branching between the synthetic turbulence and the two-phase domain significantly increases the computational time compared to simulations of Taylor bubbles performed with a laminar inlet condition showing excellent results [ 42 ]. Future development to optimize the process is ongoing.

The bubble is initiated as a hemisphere at the head, contiguous to a cylinder of a same radius letting a liquid film of a thickness δ chosen as 2.6 mm (10% of the diameter). This is significantly thicker than the final film thickness but it eases the initialization of the computation. To avoid the bubble to break during the initialization because of the fluid forces, two techniques were used:

The surface tension of the bubble was initiated with a value ten times higher than the physical one, then linearly decreased to its physical value during 0.1 s.

Three zones of uniform velocity were set to avoid excessive initial fluid force: no velocity in the bubble, higher velocity around the bubble, and nominal velocity in the rest of the domain.

The computational domain was meshed using the internal mesh generator of TRUST/TrioCFD coupled with GMSH. The process starts with meshing a quarter of disk in two different sections:

a corona filled with prisms whose size follows a geometrical expansion law. The latter are then cut to triangles.

the bulk domain homogeneously filled with triangles.

The resulted mesh is then extended to generate the mesh of the circular cross section (see Fig.  4 ) and then extruded to produce the entire pipe. The final mesh for TrioCFD simulation constituted of about 2.8 million tetrahedral elements. Details of the mesh characteristics can be found in Table  2 , η being the Kolmogorov length scale computed as η = ν 3/4 ε −1/4 with ν the liquid kinematic viscosity and ε the turbulent dissipation rate, itself estimated with the turbulent length scale approach. The subgrid turbulent model was not used, which means that the code worked with a so-called quasi-DNS approach: mesh was too coarse for DNS but sufficiently fine to allow simulations without turbulent diffusivity.

The length of the bubble after the initial non-physical interval of 0.1 s, was established at around 70 mm. Time interval of the observation was 1 s and the computation took 45,000 CPU-hours on 128 CPU cores.

3.2 Geometric VOF interface tracking in OpenFOAM

A two-phase gas–liquid system has been modeled using the one-fluid formulation of the Navier–Stokes Eqs. ( 1 ) and the geometric VOF approach for interface capturing. Within the VOF framework, a void fraction, denoted as \(\alpha\) , is defined. Its advection equation is given as:

Solution of this equation represents the starting point for the VOF reconstruction of the interface. It is important to emphasize that interface treatment with Lagrangian markers in Eq. ( 2 ) and through the volume fraction advection Eq. ( 3 ) is the key difference between both numerical models used in the present study.

In this computational study, the OpenFOAM v10 software, a widely recognized tool in the field of computational fluid dynamics (CFD), is employed to solve the relevant equations. The core of the simulation is a highly sophisticated and modified interFoam solver. This solver is notable for its capability to utilize diagonally implicit Runge–Kutta (DIRK) time integration schemes, seamlessly integrated with the piecewise linear interface calculation (PLIC) for geometric reconstruction of interfaces. This solver is an advanced iteration of an original version developed in OpenFOAM v4 by Frederix et al. [ 45 ], which was later refined by Kren et al. [ 22 ].

For modeling the subgrid-scale phenomena, the study adopts the Wall-Adapting Local Eddy-viscosity (WALE) model. This eddy viscosity model is particularly effective in capturing the dynamics of turbulent flows at smaller scales. The computational strategy includes the use of the Pressure-Implicit with Splitting of Operators (PISO) algorithm. This algorithm plays a pivotal role in the computational process, as it adeptly decouples the pressure and velocity equations, allowing for their segregated solution. A notable feature of this methodology is the incorporation of two inner corrector loops. This design implies that the pressure equation is reformulated and resolved twice within each stage of the Runge–Kutta (RK) time-stepping process.

Surface tension, a crucial aspect in multi-phase flow simulations, is computed using the continuum surface force (CSF) model [ 51 ]. This model effectively distributes the surface tension force across several computational cells, utilizing a Dirac delta function. To enhance the accuracy of this approach, the Dirac-delta function in the surface tension term is smoothened via the α function. The primary objective of this modeling technique is to precisely equilibrate the forces due to pressure gradients and surface tension, ensuring accurate representation of the physical phenomena in the simulated fluid system. For the spatial discretization of the divergence terms, the finite volume framework was used, which eventually reduces to a simple summation of all the face-normal fluxes across all the faces enclosing each control volume. Similarly to Kren et al. [ 22 ], a blended scheme was used for momentum convection term that stabilizes the artificial breakup and does not have a detrimental effect on the single-phase turbulence far away from the bubble. All other interpolations and gradients are discretized using linear schemes, which are second-order accurate. The modified solver is able to use any Runge– Kutta scheme. In the present computation, a Diagonally Implicit Runge–Kutta scheme of second order (DIRK2) was used with the CFL number in the simulations below 0.4.

To simulate a Taylor bubble in a counter-current flow under turbulent conditions, a recycling boundary condition was used at the inlet, situated upstream from the Taylor bubble. This recycling process occurs five hydraulic diameters ( D ) from the inlet, allowing sufficient space for the velocity field to develop fully. Beyond the recycling point, a distance of two to three diameters is maintained before the Taylor bubble’s nose, ensuring it doesn’t affect the boundary condition. The flow rate is continuously adjusted at each time step to balance the bubble’s buoyancy against hydrodynamic drag, keeping the bubble’s position relatively stable in the simulation. This method is depicted in Fig.  5 . To counter minor fluctuations, a gentle relaxation factor of 0.01 was employed, ensuring the bubble remains steady.

Schematic of recycling boundary condition inside the computational domain (note the direction of the gravity)

The OpenFOAM model described above was successfully used for simulation of a stagnant Taylor bubble in the counter-current water flow in a thinner pipe at Re = 1400 [ 22 ]. These results were compared with experiments performed under the same conditions. The Re = 1400 case represents laminar flow of water above the bubble, laminar flow in the liquid film and chaotic flow with developing turbulence under the tail of the bubble. The model successfully described the laminar region of the water flow and the turbulent region under the tail of the Taylor bubble. The laminar-turbulent case at Re = 1400, where bubble retains axial symmetry, was a starting point for the present, fully turbulent model at Re = 5600, and the corresponding mesh density.

In the present investigation, characterized by a Reynolds number of 5600, three distinct mesh resolutions were used, all of them prepared with “Salome Meca” tool [ 49 ], with two of them depicted in Fig.  6 . Each mesh shared an identical cylindrical shape, with a length of 0.52 m and a diameter of 26 mm. The G30 mesh consisted of around 700,000 hexahedral cells and the G15 mesh was composed of about 4.1 million cells. The naming convention of the meshes corresponds to the dimensionless spanwise cell size in the bulk flow.

Mesh G15 ( left), G30 (right)

Near the wall regions, the spanwise cell size was maintained at less than one wall unit. This wall unit, denoted as \(d{x}^{+}\) , is defined by the formula \(d{x}^{+} = dxU/\nu\) , where \(dx\) represents the cell width in actual units, \(U\) signifies the friction velocity, and ν is the kinematic viscosity. Along the streamwise direction, the cell sizes were generally kept comparable to those in the spanwise direction, albeit with additional refinement in the vicinity of the bubble.

The OpenFOAM set-up described in this section was used and verified in simulations of Taylor bubble flow in Kren et al. [ 22 ], where laminar liquid flow at Re = 1400 was prescribed at the inlet. The problem considered in [ 22 ] had simpler nature of the flow at the Taylor bubble’s body, but more complex tail behavior than in the present work: Taylor bubble in [ 22 ] was longer and experienced break-up of tiny bubbles at the flapping tail, where laminar to turbulent transition occurred. In the present work, the liquid flow is fully turbulent everywhere, but the bubble is shorter and does not exhibit a significant break-up. The intensity of the bubble breakup was actually the most important parameter that was affected by the mesh density as shown in Fig.  7 ; coarser mesh (G30) bubble has lost larger amount of mass in the same time interval in comparison with the finer mesh (G15). Further comparison of grid sensitivity study shown in Fig.  7 shows different azimuthal orientation of both bubbles, which is randomly established very early in the simulations. Other properties of the Taylor bubble, like the mean mass flow rate and shape of the nose, were rather similar on both meshes. Consequently, only the fine mesh (G15) results are presented in the next sections.

Comparison of instantaneous liquid velocity magnitude fields on G15 (fine—top) and G30 (coarse—bottom) mesh. Bright red denotes air phase. Both drawings are given at the same time t = 10 s in the simulation

The length of the bubble in the OpenFOAM simulation was around 120 mm. Time interval of the simulation was 12 s, which took 110,000 h of CPU time on 192 computer cores.

In this section, we compare one experimental and two simulated bubbles with their main properties collected in Table  3 . One can see that the bubbles are of different length. This was not planned: the initial intention was to have bubbles of the same length in both simulations performed by both groups of researchers. When it was found that bubbles in TrioCFD and OpenFOAM simulations have different lengths, both simulations were running already for a month and a decision was made to continue without changing the length. The decision was based on experimental results described in [ 37 ], which have shown that for the bubbles of the length between two and six diameters, the main properties, like the bubble velocity, shape, and water mass flow rate, do not depend on the length.

As demonstrated in [ 37 ] time-averaged bubble shape is not axisymmetric. Instead the bubble exhibits a quasi-stable asymmetric shape: the bubble is always inclined toward one side of the pipe wall. The azimuthal direction of the inclination is determined during the injection of the bubble. As discussed in [ 37 ], this behavior did not allow predictions of the time-averaged 3D shape of the Taylor bubble in turbulent counter-current flow based on the 2D shadowgraphy measurements. This asymmetry of the bubble can be easily analyzed in 3D simulations, in the 2D photographs however, it is important to perform measurements in the plane, where the bubble asymmetry is maximal. Out of around 10 different experimental cases, where asymmetry was observed at different azimuthal angles, we have selected one with the most pronounced asymmetry in the camera’s field-of-view. Even this case is not exactly perpendicular to the 2D plane of the photograph; however, it is sufficiently close and is chosen as a representative case in the present paper.

The second important parameter for comparison of the Taylor bubbles is the time interval of the observation. Our experience shows that it is ideal to have a time interval of around one minute. As seen in Table  3 , such time intervals were not achievable in simulations due to the high computational costs and/or large numerically induced breakup of the bubble.

When time interval is considered, the influence of the initial conditions in the simulations must be emphasized: both simulations start with symmetric and non-physical shape of the bubble, which is followed by a rather vigorous semi-physical transient. Consequently, a short time interval of a couple of tenths of a second must be neglected in analyses of the results. The second characteristic time in the simulations is the time interval where the initial symmetry of the bubble is broken and a quasi-stable asymmetric bubble shape is obtained. In OpenFOAM simulation the asymmetry was established after 3–5 s. After that time, the comparison of the bubble shape in experiments and simulations is feasible and the total drag of the bubbles can be compared.

As seen in SubSect.  4.1 and in Table  3 , attaining the proper steady-state bubble shape, which is ultimately responsible for the bubble’s drag, was a challenge in OpenFOAM simulation. The ultimate bubble shape predicted by the OpenFOAM had lower drag than the actual experimental bubble. Consequently, 17% higher mass flow rate was needed to keep the bubble stagnant in simulation (Table  3 ).

In the TrioCFD simulation, a fixed mass flow rate (mean water velocity v 0  = 0.18 m/s) and a moving frame of reference were used to keep the bubble in place. Since the TrioCFD simulation did not achieve the steady-state asymmetric shape, the bubble experienced stronger drag force and was moving down with velocity around 0.04 m/s. This motion was compensated with the moving frame of reference. As a rough estimate, the mean downward velocity of water, which would keep the symmetric bubble stagnant, would be around 0.13 m/s. This is about 75% of the measured downward velocity (Table  3 ).

4.1 Taylor Bubble Shape

Typical instantaneous shapes and sizes of the bubbles taken at times several seconds apart are shown in Fig.  8 , while Fig.  9 shows time-averaged bubble shapes. In the OpenFOAM simulation, the interface was defined as an isoline with values of the gas (and liquid) volume fraction equal to 0.5. In TrioCFD this definition is not needed: exact position of the interface markers is available at any time. Time averaging in the experiment was performed over a 2-min interval. In the OpenFOAM simulation, time averaging was performed in the interval from 6 to 12 s, where the bubble lost the symmetry imposed by the initial conditions and developed a roughly steady-state asymmetric shape. In the TrioCFD simulation, the analyzed time interval was only 1 s long and that was not enough for development of the asymmetric shape. Consequently, instantaneous TrioCFD bubble is the same as the time-averaged bubble, only one profile is given in Fig.  8 , and there is no TrioCFD profile in Fig.  9 .

Instantaneous silhouettes of measured Taylor bubble (left), OpenFOAM simulation (center), and TrioCFD simulation (right). Pipe walls are drawn in experimental image. Walls in simulations are left–right edges of the image

Time-averaged Taylor bubble interface position: experiment—2 min time interval (left), OpenFOAM average Taylor bubble position over 6 to 12 s time interval (right)

Time averaging takes into account slow vertical motion of the Taylor bubble in experiment and in simulations. The time-averaged bubble shapes in Fig.  9 are obtained in a moving frame of reference: origin of the coordinate system in axial direction is attached to the axial position of the bubble nose tip at every instantaneous snapshot before averaging.

The time-averaged measured bubble in Fig.  9 is inclined to the right side of the image; however, one of the three instantaneous profiles in Fig.  8 shows a bubble inclined slightly to the left. One can also note that the tail of the bubble is not resolved in all instantaneous experimental cases. The dynamics of the bottom surface is a strong 3D phenomenon, and capturing the tail position from 2D shadowgraphy does not give a particularly useful information. Consequently, the software for interface identification was not forced to work in this region.

Three instantaneous shapes of the bubble from OpenFOAM simulation are given in five seconds intervals in Fig.  8 , center: the first silhouette shown 1 s after the start of simulation is showing nearly symmetric bubble, which has not achieved the quasi-steady asymmetric shape yet. Symmetry is broken and developed at times 6 s and 11 s. These two silhouettes exhibit higher asymmetry than their experimental counterparts. Consequently, such a shape of the bubble exhibits lower drag and requires a higher mass flow rate for dynamical balance, as reflected in Table  3 .

As seen from the experimental silhouettes in Fig.  8 the bubble’s nose occasionally crosses the axis; however, on average, it remains attached to the same azimuthal angle of the pipe throughout the several minute time interval. This phenomenon is further demonstrated in Fig.  10 : 12 s time interval allows direct comparison with the OpenFOAM simulation. Wobbling of the nose is stronger in the experiment and less intensive in the OpenFOAM simulation. The important information from Fig.  10 is more intensive fluctuations of the experimental bubble in comparison with the simulated bubble from OpenFOAM simulation. This is a clear sign that the existing OpenFOAM model cannot provide a perfect description of the phenomena.

Radial position of Taylor bubble’s nose

Transition from symmetric OpenFOAM bubble to quasi-stable asymmetric one is also rather clearly seen in Fig.  10 . The simulations started with an axially symmetric Taylor bubble, and the bubble nose position has moved to the wall in about three to five seconds. After that the bubble nose stays close to the wall and the final asymmetry of OpenFOAM bubble is much more pronounced than in the measurements.

Deviations from the experiment are not without consequences: next to the higher liquid mass flow rate needed to counter-balance the OpenFOAM bubble given in Table  3 , simulated bubble creates a very thin liquid film seen on the right side of the 2D image in Fig.  9 (center). The film eventually becomes so thin that it ruptures and the air comes into the direct contact with the pipe wall. This phenomenon does not have a physical background since it was never observed in experiments. Consequently, the simulations are stopped at such instances and subgrid models of interface friction and surface tension are being investigated to improve that behavior.

The asymmetry in the bubble of the TrioCFD simulation was not achieved in the observed time interval, thus its long-time behavior cannot be predicted at this point. Nonetheless, during that time frame, radial motion of the bubble nose was observed, but it was not relevant to compare as it was only about one percent of the pipe diameter. The axial velocity of the TrioCFD bubble was 0.04 m/s downward and this was taken into account in Table  3 mass flow rate ratio calculation. This seemingly larger difference is due to the higher drag coefficient of the symmetric bubble compared to the asymmetric one. This ratio cannot be used to extrapolate the final mass flow rate that would balance the TrioCFD bubble once the quasi-stable asymmetric shape is reached.

4.2 Liquid Film Thickness and Interface Axial Velocity

The asymmetry observed in the time-averaged Taylor bubble, as depicted in Sect.  4.1 , presents challenges when it comes to independently verifying our measurements. However, a way to mitigate this issue was averaging the bubble shape and corresponding liquid film thickness over both sides of the photographs. These findings are presented in Fig.  11 , which illustrates the liquid film thickness along the bubble. The axial distances along the z -axis of Fig.  11 are measured from the bubble nose. Solid lines represent the time-averaged and left–right spatial averaged profile derived from the measurements and simulations. Dashed lines represent time-averaged left and right (thick and thin film) profiles separately, i.e., magnification of the near-wall region in Fig.  9 . The magnification of the film thickness in Fig.  11 shows that the average thickness is very similar in experiment and OpenFOAM simulation. Figure  11 is clearly exposing the feature mentioned above: bubble asymmetry is more pronounced in OpenFOAM simulation than in the experiment and the difference between the thick and thin film side is larger in OpenFOAM simulation.

Average liquid film thickness: solid lines. Dashed lines: thick and thin films. Measurement uncertainty is below 0.1 mm and is not plotted

The approximately 25% thinner film in TrioCFD simulation is due to the lower drag of the symmetric bubble, which is proportional to the reduced effective mass flow rate (effective mass flow rate is equal to imposed mass flow rate minus moving frame of reference contribution). Since the thick and thin films in TrioCFD are equal to the average, only one profile is shown in Fig.  11 .

By utilizing the averaged film thickness profiles from the measurements, an additional curve related to the mean downward liquid velocity in the film region can be derived. Based on the continuity equation and known upstream liquid velocity v 0 (0.179 m/s), the mean liquid velocity can be calculated as \(v(z)={v}_{0}{R}^{2}/[2Rh\left(z\right)-{h(z)}^{2}]\) , where R represents the pipe radius (13 mm) and h(z) represents the liquid film thickness at a given axial position z . By assuming the measured film thickness, one can calculate the mean liquid velocity within the film region. However, the quantity of interest is interface velocity and with the presented measurement techniques it cannot be measured. Nevertheless, a reasonably good approximation is available: the velocity of the interface can be estimated from the mean liquid velocity based on the liquid film thickness. The interface velocity was obtained from the mean liquid film velocity by multiplication with a factor of 1.15. This semi-empirical factor stems from the DNS simulations of the turbulent flume [ 50 ] and contains relative error around 2–3% [ 37 ]. It applies to the fully developed free liquid surface near an infinite flat wall and disregards the air shear force. As shown in [ 37 ], this approximation can be used, because the interface velocity is very similar on all sides of the bubble and does not depend on the thickness of the liquid film.

In simulations, the interface velocity can be obtained directly from the data. In the OpenFOAM results, the interface velocity is defined as a velocity at the isoline with the values of the gas (and liquid) volume fraction equal to 0.5. In the TrioCFD results, this velocity is directly obtained as a velocity of the Lagrangian marker points on the interface.

Figure  12 presents the experimental interface velocity profile for the time-averaged and spatially left–right averaged film and interface velocities from simulations. The relative uncertainty of the measured time-averaged velocity profile is similar in magnitude to the uncertainty of the mean film thickness measurement, approximately 10% at distances greater one diameter D from the bubble nose.

Taylor bubble interface velocity (m/s): OpenFOAM and TrioCFD—obtained directly from simulations. Experiment—obtained from film thickness and continuity equation. OpenFOAM: left and right interface velocity

Numerical interface velocity profiles in Fig.  12 can be easily extracted for each side of the Taylor bubble separately and are also plotted separately for OpenFOAM simulation. Only one profile is given for TrioCFD simulation in Fig.  12 due to the symmetry of the bubble.

As further explained in [ 37 ] and in Sect.  4.3 , where interface velocity is analyzed with a different approach, very similar interface velocities are expected on both sides of the Taylor bubble. Nevertheless, the OpenFOAM simulation shows very good agreement of the interface velocity only on the thick side of the liquid film. Significantly slower interface velocity is observed in the region of the thin film. Mesh resolution in OpenFOAM G15 mesh simulation describes 0.5 mm liquid film with around 10 mesh points in radial direction. This is actually coarse for an accurate description of a turbulent film with LES. If we add that the interface is smeared over 2 to 3 points, we can conclude that the higher radial resolution or a more elaborated subgrid model is needed around such interface.

The velocity profile of the TrioCFD simulation is close to the measurements and within the uncertainty of the measurements, which are characterized with a single uncertainty bar at around z  = 70 mm. It is clear that the markers are accurate in specifying the location and the velocity of the interface. Nevertheless, the high similarity between the TrioCFD and measured interface velocities is due to the fact that the walls in the TrioCFD simulation are moving up at around 0.04 m/s along with the moving frame of reference, while the bubble is fixed in space.

4.3 Axial Velocities of Disturbance Waves on the Interface

The presented measurement techniques and image processing algorithms allow us to track small disturbance waves traveling along the Taylor bubble interface [ 43 , 44 ]. This technique was used and described by Kren et al. [ 37 ]. The specific mechanisms generating these waves are not entirely clear, but the waves are believed to be induced by the turbulence. However, assuming that most of the waves are produced by random disturbances, they are expected to travel in all directions parallel to the air–water interface. The velocities of these waves are governed by the capillary wave equations, as described in [ 5 ]. The dispersion relation of capillary waves can be expressed as:

\(\omega^{2} = \frac{{{ }\sigma k^{3} }}{\rho }{\text{tanh}}\left( {k d} \right)\) ,

where \(k=2\pi /\lambda\) is the wavenumber, \(\omega =\) \(2\pi \nu\) angular frequency, and \(c=\) \(\lambda \nu\) is the phase velocity. For typical “thick” 2 mm liquid film waves, the characteristic frequencies, wavelengths, and phase velocities are approximately 1 Hz, 50 mm, and 0.05 m/s, respectively. For the thinner 0.5 mm liquid film, these values are approximately 40 Hz, 5 mm, and 0.2 m/s. These estimates indicate that the characteristic phase velocities of the waves are lower than the interface velocities shown in Fig.  12 , implying that practically all waves on the interface travel downward.

To estimate the axial velocities of the disturbance waves traveling over the interface, measurements of the axial disturbance wave velocity w are performed using cross-correlations of the time signals at various axial positions along the pipe. By selecting a distance H between specific points in space, for example, H  = 200 pixels, the velocity w can be obtained from the measured time lag τ of the signals as w  =  H/τ . For example, the time lag at the point 400 pixels downstream of the bubble nose and at a distance H  = 200 pixels, is computed from the cross-correlation of time signals at points 400 −  H /2 = 300 pixels and 400 +  H /2 = 500 pixels. The procedure is described in [ 37 ], where it was applied to the measurements.

The same procedure was used also for analyses of the disturbance wave velocities in the OpenFOAM simulation results and the results are presented in Fig.  13 together with the interface velocities. Velocities are given for one real and two simulated Taylor bubbles listed in Table  3 . The disturbance wave speed profiles are derived from the time lags observed in the right-hand side of the silhouettes, except for TrioCFD results, where both sides are symmetric and the results are very similar and not repeated twice.

Disturbance velocity (dashed) vs. interface velocity (solid) in the experiment (top-left), TrioCFD (top-right), and OpenFOAM (bottom-left [0.3 s:3.6 s] interval and bottom-right [8.4 s:11.8 s] interval)

The complete time history consisting of 50,000 frames over 125 s interval is analyzed for measurement. In the OpenFOAM results two intervals of the same length 3.4 s at the beginning and at the end of simulation (50–2050 frames and 5000–7000 out of total 7000 frames) were used. In the TrioCFD simulation around 1700 frames were analyzed over a 1.7 s interval. Cross-correlations are compared at a fixed distance of around 11 mm in measurements and in OpenFOAM, and 20 mm in TrioCFD. The discrete values of cross-correlation time lags are smoothed using parabolic interpolation.

Figure  13 shows disturbance velocities and the corresponding interface velocities in experiment and in simulations. They are separated into four separated graphs plotted on the same spatial scale and with the same velocity range in order to make the differences and similarities clear. Both distinct types of velocity profiles are obtained from the same measurements but through entirely different analysis. The interface velocities are obtained directly from the simulations and from the measurement of the liquid film thickness in the experiment. On the other hand, the disturbance velocities are determined based on the relative motions of the liquid–air interface. Notably, both types of velocities exhibit remarkable similarity in most of the graphs in Fig.  13 . This observation confirms the hypothesis given in [ 37 ] that the time-averaged velocity of the disturbance waves on the water–air interface effectively represents the velocity of the interface itself.

Before proceeding to the further discussion it is important to address the discrepancies seen in the graphs of Fig.  13 . Disturbance velocity profiles obtained in experiment show very similar disturbance wave velocities on both sides of the photograph, despite significant difference in the liquid film thickness seen in Fig.  8 and in Fig.  11 . Similar observation was reported in other experimental cases in [ 37 ]. Further observation of experimental profiles shows increasing discrepancies on the right side of the photographs in the region beyond 70 mm from the bubble’s nose. This discrepancy appears in the region of poorer spatial resolution in the thin film, which is resolved on 2–5 pixels. Very thin liquid film filters the long wavelengths, which are the most relevant for accurate evaluation of the disturbance velocity measurements. Short wavelengths, which remain on the thin films, are more difficult to resolve with the limited resolution of the presented experiment. The only solution is to increase the spatial resolution of the photographs.

The disturbance velocities in TrioCFD were obtained in 8 points along the bubble, which are marked with symbols in Fig.  13 . Computation of cross-correlation functions is difficult and sensitive to the noise even on the fixed meshes of the photographs and OpenFOAM results. The problem becomes even more difficult for TrioCFD, where cross-correlation functions must be obtained from the moving markers on the interface. Consequently, for the given spatial resolution, only very rough estimates of disturbance velocities are available in TrioCFD (Fig.  13 ).

Lastly, we need to comment on the disturbance wave velocities computed from the OpenFOAM results using the same procedure as in the experiment. Two separate graphs are shown in Fig.  13 for OpenFOAM, one for the interval at the beginning of the simulation during the time interval [0.3 s: 3.6 s] where the bubble is losing its symmetry, and the other for the time interval [8.4 s: 11.9 s] where the bubble is close to quasi-steady-state and with a very thin liquid film on the left side plane shown in Fig.  8 . Disturbance velocity profiles computed in the first time interval where the Taylor bubble is close symmetric are not smooth, but reasonably close to the computed interface velocities. Coarse, ~ 0.05 mm, radial resolution in the near wall region can only provide a very rough approximations of the disturbances with amplitudes 0.01–0.05 mm.

The problem is more exaggerated in the time interval at the end of the OpenFOAM simulation. At that time, the liquid film becomes very thin on one side and rather thick on the other side. As shown in Fig.  13 the cross-correlation technique still works for the thin film (not with a great precision), but fails on the thick side. The reason for this failure is film position, which is not within the finely resolved boundary layer visible in G15 mesh of Fig.  6 , but in the coarser central region where small disturbance waves cannot be captured anymore.

The equivalence between the time-averaged velocities of the interface waves and the convective velocity of the interface observed in experiments, is thus roughly confirmed also in simulations. The main reason for the equivalence lies in the fact that the characteristic velocities of the dominant disturbance waves are significantly lower, at least by an order of magnitude, compared to the mean velocities of the liquid film (approximately 1 m/s). As a result, the time averaging process predicts the final disturbance wave velocity equal to the time-averaged convective velocity of the water–air interface.

The final comparison of the experiment and simulations is focused on spectra of the disturbance waves traveling in the axial direction over the body of the Taylor bubble. Figure  14 shows power spectra of the disturbance waves analyzed in a point approximately 50 mm below the Taylor bubble’s nose. The common property of measured and TrioCFD spectra is very sharp drop at frequencies above 10–20 Hz. An exception is seen in the OpenFOAM results and we assume that the difference is due to the implemented numerical scheme in combination with the LES-WALE mode, which does not sufficiently suppress the fluctuations in the range of frequencies between 10 and 70 Hz.

Power spectra of interface disturbance waves at point 50 mm downstream the bubble’s nose

The highest frequencies of the turbulent fluctuations in the single-phase flow of water above the bubble can be estimated from the DNS database of Kasagi (Fukata and Kasagi, 2002) pipe flow case at Re = 5300, which is close to the present experimental conditions. Their highest frequencies of Kolmogorov scale vortices are between 10 and 70 Hz in the axis of the pipe and in the near-wall region, respectively. Since the turbulent kinetic energy in the vortices at the Kolmogorov scales is low, frequencies between 3 and 20 Hz, which correspond to the Taylor microscales, seem to be more relevant as the upper limits.

Spectral analysis of the waves traveling over the Taylor bubble’s body shown in Fig.  14 is representative also at other distances from the bubble’s nose.

5 Conclusions

This paper summarizes the studies of the Taylor bubble in a vertical turbulent counter-current air–water flow, excluding the bubble’s tail region. The analyses are based on results available from high-precision 2D shadowgraphy observations which are compared with high-fidelity simulations with two different computer codes and numerical schemes: OpenFOAM is using geometrical VOF interface tracking, while the TrioCFD code is based on explicit front tracking method. Considered Taylor bubble was observed in the inertia dominant regime, where the influence of viscosity and surface tension are minor.

The primary goal of this study was to compare the time-averaged shape of the Taylor bubble’s interface. From the comparison between the experimental data and the OpenFOAM simulation—which used a model verified and validated on the case of a stagnant Taylor bubble in a laminar counter-current flow [ 22 ] —two main findings emerged:

The bubble asymmetry in the simulation was more pronounced than that observed in the experiments, resulting in lower bubble drag and, consequently, a roughly 20% higher liquid mass flow rate needed to keep the bubble stationary.

A less critical, but still relevant issue: fluctuations around the time-averaged bubble shape were weaker in the simulation than in the experiment.

The simulation of TrioCFD was performed with a more expensive numerical approach and only one second of transient was analyzed. This was not enough to develop an asymmetric bubble shape. Consequently, only qualitative bubble shape comparison was performed.

The second part of the study was focused on dynamics of interfacial waves traveling over the body of the Taylor bubble. Our analyses of various experiments in [ 37 ] have shown that the disturbance wave velocity, measured over a sufficiently long interval of several tens of seconds, becomes equal to the axial water–air interface velocity. The cross-correlation measurements primarily capture low-frequency waves, which are slower than the interface velocity. Therefore, tracking these waves provides a technique for measuring the time-averaged interface velocity.

The same analysis of disturbance waves on the interface was performed in both numerical simulations, where the accuracy of the analyses was severely limited with the spatial discretization of both simulations. Refined mesh in the near-wall region was barely sufficient to capture the disturbance waves and to reconstruct their propagation velocities in OpenFOAM. The specific numerical approach of TrioCFD was even less appropriate for disturbance velocity measurements and allowed only rough approximation. When the interface fell out of the refined mesh in the near-wall region into the coarse meshing in the center of the pipe in OpenFOAM simulation, disturbance waves were not recognized anymore.

Spectra of the interfacial waves were compared at a point fixed from the bubble’s nose. TrioCFD showed greater precision than OpenFOAM, accurately reflecting the sharp frequency decline above 10–20 Hz observed in experiments. In contrast, OpenFOAM’s spectra erroneously displayed significant frequencies up to 70 Hz, likely misrepresenting the physical phenomena.

The stagnant Taylor bubble in turbulent low-Reynolds counter-current flow was identified as a challenging test case for two advanced interface tracking models in the TrioCFD and OpenFOAM codes, even though they were both well validated in laminar conditions [ 22 , 41 ]. We demonstrated that fine spatial resolution is necessary not only in the near-wall region, where the liquid boundary layer forms, but also at the interface itself. Future simulations with both codes will aim to enhance the subgrid models for interface friction and surface tension.

Wallis G.B.: One-Dimensional Two-Phase Flow. McGraw Hill, (1969)

Morgado, A.O.; Miranda, J.M.; Araújo, J.D.P.; Campos, J.B.L.M.: Review on vertical gas-liquid slug flow. Int. J. Multiph. Flow 85 , 348–368 (2016)

MathSciNet   Google Scholar  

Zhou, G.; Prosperetti, A.: Violent expansion of a rising Taylor bubble. Phys. Rev. Fluids 4 , 073903 (2019)

Google Scholar  

Holagh, S.G.; Ahmed, W.H.: Critical review of vertical gas-liquid slug flow: an insight to better understand flow hydrodynamics’ effect on heat and mass transfer characteristics. Int. J. Heat Mass Transf. 225 , 125422 (2024)

Liberzon, D.; Shemer, L.; Barnea, D.: Upward-propagating capillary waves on the surface of short Taylor bubbles. Phys. Fluids 18 , 048103 (2006)

Dumitrescu, D.T.: Strömung an einer Luftblase im senkrechten Rohr. Z. angew. Math. Mech. 23 , 139 (1943)

Martin, C.S.: Vertically downward two-phase slug flow. ASME J. Fluids Eng. 98 (4), 715 (1976)

Lu, X.; Prosperetti, A.: Axial stability of Taylor bubbles. J. Fluid Mech. 568 , 173–192 (2006)

Figueroa-Espinoza, B.; Fabre, J.: Taylor bubble moving in a flowing liquid in vertical channel: transition from symmetric to asymmetric shape. J. Fluid Mech. 679 , 432–454 (2011). https://doi.org/10.1017/jfm.2011.159

Fabre, J.; Figueroa-Espinoza, B.: Taylor bubble rising in a vertical pipe against laminar or turbulent downward flow: symmetric to asymmetric shape transition. J. Fluid Mech. 755 , 485–502 (2014). https://doi.org/10.1017/jfm.2014.429

Fershtman, A.; Babin, V.; Barnea, D.; Shemer, L.: On shapes and motion of an elongated bubble in downward liquid pipe flow. Phys. Fluids 29 , 112103 (2017). https://doi.org/10.1063/1.4996444

Abubakar, H.; Matar, O.: Linear stability analysis of Taylor bubble motion in downward flowing liquids in vertical tubes. J. Fluid Mech. 941 , A2 (2022). https://doi.org/10.1017/jfm.2022.261

Delfos, R.; Wisse, C.J.; Oliemans, R.V.A.: Measurement of air-entrainment from a stationary Taylor bubble in a vertical tube. Int. J. Multiph. Flow 27 , 1769–1787 (2001)

Kockx, J.P.; Nieuwstadt, F.T.M.; Oliemans, R.V.A.; Delfos, R.: Gas entrainment by a liquid film falling around a stationary Taylor bubble in a vertical tube. Int. J. Multiph. Flow 31 , 1–24 (2005)

Mikuž B., Kamnikar, J., Prošek, J., Tiselj, I.: Experimental observation of Taylor bubble disintegration in turbulent flow. Proc. In: 28th Int. Conf. Nuclear Energy for New Europe 9, (2019).

Crowe, C.T.; Troutt, T.R.; Chung, J.N.: Numerical Models for Two-Phase Turbulent Flows. Annu. Rev. Fluid Mech. 28 , 11–43 (1996)

Trujillo, M.F.: Reexamining the one-fluid formulation for two-phase flows. Int. J. Multiph. Flow 141 , 103672 (2021)

The OpenFOAM Foundation, OpenFOAM | Free CFD Software, (2022).

Hirt, C.; Nichols, B.: Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 39 , 201–225 (1981)

Sint, A.M.; Deen, N.; Kuipers, J.: Numerical simulation of gas bubbles behaviour using a three-dimensional volume of fluid method. Chem. Eng. Sci. 60 , 2999–3011 (2005)

Tryggvason, G.; Scardovelli, R.; Zaleski, S.: Direct Numerical Simulations Of Gas-Liquid Multi-Phase Flows. Cambridge University Press, Cambridge, New York (2011)

Kren, J.; Frederix, E.M.A.; Tiselj, I.; Mikuž, B.: Numerical study of Taylor bubble breakup in countercurrent flow using large eddy simulation. Phys. Fluids 36 (2), 023311 (2024). https://doi.org/10.1063/5.0186236

Cifani, P.; Michalek, W.; Priems, G.; Kuerten, J.; Geld, C.; Geurts, B.: A comparison between the surface compression method and an interface reconstruction method for the VOF approach. Comput. Fluids 136 , 421–435 (2016)

Dai, D.; Tong, A.Y.: Analytical interface reconstruction algorithms in the PLIC-VOF method for 3D polyhedral unstructured meshes. Int. J. Numer. Meth. Fluids 91 , 213–227 (2019)

Klein, M.; Ketterl, S.; Hasslberger, J.: Large eddy simulation of multiphase flows using thevolume of fluid method: part 1—governing equations and a priori analysis. Exp. Comput. Multiph. Flow 1 , 130–144 (2019)

Vreman, A.W.: An eddy-viscosity subgrid-scale model for turbulent shear flow: algebraic theory and applications. Phys. Fluids 16 , 3670–3681 (2004)

Unverdi, S.O.: A front-tracking method for viscous, incompressible, multi fluid flows. J. Comput. Phys. 100 (1), 25–37 (1992)

Tryggvason, G.; Bunner, B.; Esmaeeli, A.; Juric, D.; Al-Rawahi, D.; Tauber, W.; Han, J.; Nas, S.; Jan, Y.: A front-tracking method for the computations of multiphase flow. J. Comput. Phys. 169 (2), 708–759 (2001)

Mathieu, B.: Physical, experimental and numerical study of fundamental mechanisms involved in two-phase flows, Ph.D. Thesis, Université Aix-Marseille 1, 2003

Osher, S.; Fedkiw, R.P.: Level Set Methods: An Overview and Some Recent Results. J. Comput. Phys. 169 , 463–502 (2001)

Araújo, J.; Miranda, J.; Pinto, A.; Campos, J.: Wide-ranging survey on the laminar flow of individual Taylor bubbles rising through stagnant Newtonian liquids. Int. J. Multiph. Flow 43 , 131–148 (2012)

Morgado, A.; Miranda, J.; Araújo, J.; Campos, J.: Review on vertical gas–liquid slug flow. Int. J. Multiph. Flow 85 , 348–368 (2016)

Cerqueira, R.F.; Paladino, E.E.; Evrard, F.; Denner, F.; Wachem, B.: Multiscale modeling and validation of the flow around Taylor bubbles surrounded with small dispersed bubbles using a coupled VOF-DBM approach. Int. J. Multiph. Flow 141 , 103673 (2021)

Gutiérrez, E.; Balcázar, N.; Bartrons, N.; Rigola, J.: Numerical study of Taylor bubbles rising in a stagnant liquid using a level-set/moving-mesh method. Chem. Eng. Sci. 164 , 158–177 (2017)

Angeli P.E., Puscas M.A., Fauchet G., Cartalade A.: FVCA8 Benchmark for the stokes and navier-stokes equations with the TrioCFD Code. finite volumes for complex applications VIII - methods and theoretical aspects, (2017)

Mikuž B., Frederix E.M.A., Komen E.M.J., Tiselj I.: Taylor bubble behaviour in turbulent flow regime. Proceedings of the conference Computational Fluid Dynamics for Nuclear Reactor Safety (CFD4NRS-8), 12 (2020)

Kren, J.; Zajec, B.; Tiselj, I.; El Shawish, S.; Perne, Ž; Tekavčič, M.; Mikuž, B.: Dynamics of Taylor bubble interface in vertical turbulent counter-current flow. Int. J. Multiphase Flow 165 , 104482 (2023)

Tekavčič, M.; Končar, B.; Kljenak, I.: The concept of liquid inlet model and its effect on the flooding wave frequency in vertical air-water churn flow. Chem. Eng. Sci. 175 , 231–242 (2018). https://doi.org/10.1016/j.ces.2017.09.050

Xue, Y.; Stewart, C.; Kelly, D.; Campbell, D.; Gormley, M.: Two-phase annular flow in vertical pipes: a critical review of current research techniques and progress. Water 14 , 3496 (2022). https://doi.org/10.3390/w14213496

Slavchov, R.I.; Peychev, B.; Ismail, A.S.: Characterization of capillary waves: a review and a new optical method. Phys. Fluids 33 (10), 101303 (2021)

Giamagas, G.; Zonta, F.; Roccon, A.; Soldati, A.: Propagation of capillary waves in two-layer oil–water turbulent flow. J. Fluid Mech. 960 , A5 (2023). https://doi.org/10.1017/jfm.2023.189

Nop, R., Hamrit, G., Burlot, A., Bois, G., Mikuž, B., Tiselj, I.: The 3D DNS of a Taylor bubble in counter-current flow with a turbulent wake using the Front-Tracking method in TrioCFD. In: 32nd International conference : nuclear energy for new europe, Portorož, Slovenia, September 11–14, (2023).

Pan, L.; He, H.; Ju, P.; Hibiki, T.; Ishii, M.: Experimental study and modeling of disturbance wave height of vertical annular flow. Int. J. Heat Mass Transf. 89 , 165–175 (2015)

Lin, R.; Wang, K.; Liu, L.; Zhang, Y.; Dong, S.: Study on the characteristics of interfacial waves in annular flow by image analysis. Chem. Eng. Sci. 212 , 115336 (2020)

Frederix, E.M.A.; Komen, E.M.J.; Tiselj, I.; Mikuž, B.: LES of turbulent co-current Taylor Bubble flow. Flow Turbulence Combust. 105 , 471–495 (2020)

Taha, T.; Cui, Z.F.: CFD modelling of slug flow in vertical tubes. Chem. Eng. Sci. 61 (2), 676–687 (2006). https://doi.org/10.1016/j.ces.2005.07.022

Press W.H., Teukolsky, S.A., Vetterling, W.T., Flannery B.P.: Numerical recipes 3rd edition: the art of scientific computing. Cambridge Press, (2007).

Grishchenko, D.: KROTOS image analysis for water-corium interactions (KIWI). OECD SERENA project report DEN/DTN/STRI/LMA/NT/2011/009/0, CEA, France, (2011).

Consortium, Open source. “Salome Meca”. Version 9.*, http://www.salome-platform.org/ (2023). Accessed 2023

Bergant, R.; Tiselj, I.: Near-wall passive scalar transport at high Prandtl numbers. Phys. Fluids 19 , 065105 (2007)

Brackbill, J.; Kothe, D.; Zemach, C.: A continuum method for modeling surface tension. J. Comput. Phys. 2 , 335–354 (1992)

Download references

Acknowledgements

The authors gratefully acknowledge financial support provided by Slovenian Research Agency, grant P2-0026 and Slovenia-CEA grant NC-0026. The TrioCFD computations were made possible by the granted access to the HPC resources of IDRIS under the allocation 2023-R0131010339 made by GENCI.

Funding was provided by Javna Agencija za Raziskovalno Dejavnost RS, NC-0026, Iztok Tiselj, P2-0026

Author information

Authors and affiliations.

Jožef Stefan Institute, Jamova 39, 1000, Ljubljana, Slovenia

Iztok Tiselj, Jan Kren & Blaž Mikuž

Faculty of Mathematics and Physics, University of Ljubljana, Jadranska Ulica 19, 1000, Ljubljana, Slovenia

Iztok Tiselj & Jan Kren

Service de Thermohydraulique Et de Mécanique Des Fluides, Université Paris-Saclay, CEA, 91191, Gif-Sur-Yvette, France

Raksmy Nop, Alan Burlot & Grégoire Hamrit

You can also search for this author in PubMed   Google Scholar

Corresponding author

Correspondence to Iztok Tiselj .

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/ .

Reprints and permissions

About this article

Tiselj, I., Kren, J., Mikuž, B. et al. Experimental and Numerical Study of Taylor Bubble in Counter-Current Turbulent Flow. Arab J Sci Eng (2024). https://doi.org/10.1007/s13369-024-09489-2

Download citation

Received : 11 March 2024

Accepted : 28 July 2024

Published : 21 August 2024

DOI : https://doi.org/10.1007/s13369-024-09489-2

Share this article

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

• Taylor bubble
• Shadowgraphy
• VOF interface capturing
• Front tracking
• Find a journal
• Publish with us
• Track your research