Modeling a Sphere Falling on a Water Surface

A sphere slightly lighter than water is dropped into a small liquid-filled beaker. How do you model the resulting shape of the water surface and the motion of the sphere as it floats on the surface? In this blog post, we demonstrate how to model this system in the COMSOL Multiphysics® software. Although we are considering the case of a small sphere, the technique discussed here can be applied to other larger shapes as well.

Oscillating Motion of a Buoyant Sphere: A CFD Problem

Let’s consider a buoyant sphere that is dropped into a water-filled beaker. The sphere, with a radius of 0.15 cm and density of 800 kg/m³, is initially suspended in air. The cylindrical beaker has a radius of 0.65 cm, a height of 1.7 cm, and is partially filled with water. The size of the system is small enough to assume laminar flow and also allows for a faster solution process on a regular PC. A larger system can be modeled in exactly the same way, but would require more computational resources.

To model the motion of the sphere through the water surface, we need to model the fluid flow in both the air and water phases. If we consider the sphere to be rigid, we can compute the motion of it by solving the ordinary differential equations (ODEs) of motion that account for the force experienced by the sphere falling through the air, hitting the water surface, and finally interacting with both air and water. The total fluid stress exerted by the fluids on solid boundaries is given by

(1)

where u_fluid denotes the fluid velocity field, p the fluid pressure, μ the dynamic viscosity of the fluid, n the outward normal to the boundary, and I the identity matrix.

Buoyant sphere suspended in air above the water surface in a water-filled beaker. Part of the beaker wall has been hidden to show the sphere.

Assuming that the solid only moves in the vertical direction (z), we can write

(2)

where m denotes the mass, v_s the velocity in the z direction, and u_s the displacement in the z direction of the sphere.

The force F_z can be calculated by integrating the z-component of \mathbf{f} from Eq. 1 over the surface of the sphere. The above equations of motion (ODEs) are solved along with the finite element model to compute the displacement and velocity of the sphere.

Since the sphere moves through both air and water, a two-phase flow model is formulated on a spatial frame with time-dependent spatial coordinates to model both fluids. The spatial frame is the usual fixed global Euclidean system with spatial coordinates (x, y). These spatial coordinates, corresponding to the mesh nodes of the fluid domains, are functions of time and represent the current location of a point in space. This is achieved by solving the fluid flow problem on a Moving Mesh interface.

The mesh, on which the fluid flow problem is solved, deforms with a mesh deformation that is equal to the displacement of the sphere at the interfacing boundaries between the sphere and the fluid domains. Inside the fluid domains, the mesh deformation is computed using an appropriate smoothing method that is available in the Moving Mesh interface.

(3)

On the fluid boundaries interfacing the sphere, the fluid velocity is given by

(4)

Setting Up the Model in COMSOL Multiphysics®

We can set up the above problem as a 2D axisymmetric analysis. To model this problem, we couple a Two-Phase Flow interface with a Moving Mesh interface and Global Equations by defining the coupling expressions illustrated in Eqs. 1–4.

Geometry of the 2D axisymmetric model showing the initial position of the water surface represented by the horizontal line. The sphere is represented as a boundary in the model domain because we are modeling it as a rigid object using ODEs.

The fluid flow in this case is laminar, in both the air and water phases, so we can use the Laminar Two-Phase Flow, Phase Field interface. (Read this blog post on how to choose a multiphase flow interface for more information.) This multiphysics interface automatically defines a coupling between the laminar flow and phase field method, which is used to track the detailed shape of the water surface.

The Laminar Two-Phase Flow, Phase Field interface enables us to model the boundaries of the sphere as moving wetted walls. The Wetted Wall boundary condition of the Phase Field interface lets us specify the contact angle between the water surface and the solid walls, while the Wall Movement setting in the Wall boundary condition of the Laminar Flow interface lets us specify the velocity of the moving solid walls. An Outlet boundary condition is used for the boundary corresponding to the top of the beaker.

The structure of the model setup and the various interfaces — namely Laminar Flow, Phase Field, and Moving Mesh — is shown in the screenshot below:

In the Phase Field interface, we specify the value of the phase field variable as Φ = 1 (air for Fluid 1) and Φ = -1 (water for Fluid 2) for the initial configuration, along with the boundary representing the initial position of the water surface with the air above and the water below the horizontal line. The initial configuration is shown in the figure above.

We use a Wetted Wall boundary condition to specify the contact angle between the air-water interface, solid walls of the beaker, and sphere. For simplicity, we assume the same contact angle for both the sphere and walls of the beaker.

The velocity field in the fluids is automatically coupled to the phase field method to account for advection. The properties of the fluids and surface tension coefficient are specified in the multiphysics coupling between the Laminar Flow and Phase Field interfaces.

The Moving Mesh interface solves for the mesh deformation within the fluid domains due to movement of the sphere. We also use boundary conditions to describe how the mesh at the fluid boundaries should deform. The important boundary condition here is Mesh Displacement, which is equal to the sphere’s displacement on the spherical boundaries (Eq. 3).

The movement of the sphere is purely due to the gravitational load. However, the fluid exerts a drag force as defined by Eq. 1. The z-component of this force can be computed by integrating the built-in variable for total fluid stress in the z direction over the sphere boundaries. The built-in variable spf.T_stressz computes the total stress in the z direction that is exerted on the sphere boundary by the surrounding fluid.

Learn more in this blog post on how to compute lift and drag.

The ODEs (Eq. 2) are solved using the Global Equations interface. The screenshot below shows one of the equations implemented in the Global Equations interface, where vs is the velocity of the sphere as noted in Eq. 2, vst is the time derivative of velocity (or the acceleration), and intop1 is an integration operator defined over the solid boundaries.

Settings of the Global Equations interface used to solve for the velocity of the sphere.

Since there is a large movement of the sphere, the mesh deformation will be large. To help keep the element quality high and avoid inverted elements, we can use the Automatic Remeshing feature, which stops a transient simulation based on a mesh quality metric and remeshes the current deformed shape before it continues.

Analyzing the Simulation Results

The plot below shows the sphere falling on the water surface, the water-air interface (Φ = 0.5), and a cross-sectional plot of the fluid velocity magnitude.

Animation of the movement of the buoyant sphere in the water-filled beaker, the water surface, and the fluid velocity magnitude plotted on a cross-sectional slice. To see the variation in velocity magnitude at later stages of the simulation, the maximum value of the color scale is capped at 0.1[m/s].

The figure below shows the plot of the displacement of the sphere versus time. As expected, its motion is damped by the fluid surrounding it and the displacement amplitude decreases exponentially with time.

Concluding Thoughts

In this blog post, we have demonstrated how to model the oscillating motion of a buoyant sphere. Here, the sphere and its motion are considered axisymmetric with no rotational motion. In general, should the situation require, the same strategy discussed here can be applied to model both rotation and translation. Modeling such a case would require solving for additional rotational motion of the rigid object using ODEs. It would be much easier to instead use the built-in functionality of the Solid Mechanics interface to automatically solve for the rigid body motions. This would also allow for modeling a larger system where the flow becomes turbulent and the motion can take place in all directions.

Access the file for the model featured in this blog post by clicking the button below: