Santa Barbara simulations
This page describes some preliminary simulation results carried out at Santa Barbara, to model a freely moving and deforming object in two dimensions, where the boundary is captured using a level set description. The equations which are being simulated are
Here, u and v are the horizontal and vertical components of velocity. The remaining three components describe the stress tensor,
where p is the pressure, and s and τ are the deviatoric components. At the boundary of the object, the normal/normal and normal/tangential components of the stress tensor are zero. For the current simulations, all parameters in the equations were set to 1.
The computational implementation
The code was written in C++ and makes use of Frederic Gibou's level set library. In addition to the fields above, a function φ is defined, which is negative in the interior of the shape, and positive in the exterior. The boundary of the object is given by the zero contour of φ. As the simulation progresses, the boundary moves through different grid points, and a number of computaional issues need to be addressed:
- Storing of the field values: The values of the five fields are stored at all points in the interior the object. In addition, they are stored at all the first neighbors of the interior points. This allows for the specification of the fields on the exact edge of the object by linear interpolation.
- Computation of derivatives: In the interior of the shape, derivatives are computed using finite differences. The advection terms in the equations are calculated using one-sided derivatives, while all the rest are calculated using centered differences. At the boundary of the shape, the field values may not exist, so this cannot be employed. Instead, the code locally solves a 2D linear regression problem on all the available field values in a 3 × 3 grid surrouding each point.
- Passing the velocity field to the level set library: The velocity fields are only defined within the shape, so before the level set routine is called, the code extrapolates their values to all the grid points which are within three hops of the defined values. The extrapolation is done by solving a local linear regression problem on all available field values in a 5 × 5 grid surrounding each point.
- Updating of grid points: When the level set is updated, some points may be deleted from the object, and some new grid points may be created. At the created points, the code initializes the field values by solving a local regression problem as described above.
- Boundary conditions: The boundary condition creates a coupling between the three stress variables p, s, and τ which has to be well-captured by the simulation in order to produce valid results. Previously, simulations have tried to impose the boundary condition at all grid points which are at the edge of the object, but this can create large errors, since often the boundary of the shape (as described by the level set) passes between grid points. For these simulations, and alternative approach was taken. The code finds all pairs of neighboring points where one lies in the object, and one lies outside, so that the boundary of the object passes between them. The stress values are linearly interpolated to the boundary, and the boundary condition is calculated there. The field values at the two grid points are then modified so that the condition is exactly upheld at the boundary.
For the simulations presented below, a 64 × 64 grid was employed. In all three cases, the initial stress tensor was set to zero, and a velocity field was initialized to create motion.
Demo 1 – A collapsing and expanding circle
The first test consists of a circle of radius 0.6, where the velocity field intially corresponds to uniform contraction, so that u(x,y)=-x and v(x,y)=-y. The image below shows the pressure that is set up in the body to oppose the contraction.
The gradient above eventually causes the object to rebound. In the movies below, the object undergoes several oscillations. By the end of the simulation, the boundary of the object is no longer perfectly circular. A longer simulation or a smaller timestep may alleviate this.
Demo 2 – A rotating ellipse
The second test consists of an ellipse, where the initial velocity field is a rotation, so that u(x,y)=-y and v(x,y)=x. The image below shows the pressure gradient during the motion, which becomes negative to oppose the rotational forces:
As the simulation progresses, the object elongates and loses its elliptical shape. A longer simulation with more grid points may be able to capture this shape more accurately.
Demo 3 – A localized push
In this simulation, an ellipse experiences an initial localized vertical velocity, so that v(x,y)=0.3-|x| if |x|<0.3, and v(x,y)=0 otherwise. The image below shows the pressure, after the push has taken place. As expected, a negative pressure is set up in the center-top of the object, and a positive pressure as set up at the center-bottom:
As the simulation progresses, the pressure differences become smaller, and the body begins to move uniformly upwards, as would be expected by conservation of momentum. Some small oscillations are visible in the pressure at late times which may be a numerical artefact; this needs to be studied in more detail.