Circular Inclusion Benchmark

Schmid and Podladchikov [Clast] derived a simple analytic solution for the pressure and velocity fields for a circular inclusion under simple shear as in Figure 1.

Figure 1: Schematic for the circular inclusion benchmark

The file input/benchmarks/circular_inclusion/README has instructions on how to run this benchmark.

Because of the symmetry of the problem, we only have to solve over the top-right quarter of the domain. For the velocity boundary conditions, the analytic solution is a bit complicated. So we used the simple relation

[;\begin{align*} v_x &= -\dot{\epsilon}y, \\ v_y &= \dot{\epsilon}x, \end{align*};]

for the boundaries, where [;\dot{\epsilon}=1;] is the magnitude of the shear and [;x;] and [;y;] are the coordinates. This induces an error of order [;r_i^2/r^2;], where [;r_i=0.1;] is the radius of the inclusion, and [;r;] is the radius. We have the boundaries at 80 times the radius of the inclusion, giving an error of about 0.01%, which is much smaller than the other errors we were looking at. Just to make sure, we did runs with boundaries at 40 times the radius of the inclusion and got very similar results.

A characteristic of the analytic solution is that the pressure is zero inside the inclusion, while outside it follows the relation

[;p_m=4\dot{\epsilon}\frac{\mu_m(\mu_i-\mu_m)}{\mu_i+\mu_m}\frac{r_i^2}{r^2}\cos(2\theta),;]

where [;\mu_i=2;] is the viscosity of the inclusion and [;\mu_m=1;] is the viscosity of the background media. Many numerical codes that solve Stokes flow, including Gale, assume that pressure, velocity, and viscosity are continuous. The pressure discontinuity at the surface of the inclusion violates that assumption, so the error tends to concentrate near the surface of the inclusion.

Figure 2 plots the error in the pressure along the line [;y=x/2;] for different resolutions. Inside the inclusion near the surface, the pressure is consistently wrong. The pressure does not converge with higher resolution, giving us a clue that the numerical scheme is not completely accurate.

Figure 2: Pressure along the line [;y=x/2;] for resolutions of 128 × 128 (blue), 256 × 256 (red), and 512 × 512 (black). The inclusion has a radius [;r_i=0.1;]. Note that the pressure should be zero inside the inclusion, but the numerical solutions consistently underestimate the pressure.

Outside the inclusion, the error is better behaved. Figure 3 plots the error in the pressure along the line [;y=x/2;] outside the inclusion for different resolutions. While there are still problems near the surface, away from the surface the solutions are quite good. Figure 4 plots the error scaled with resolution, and we can see that the error scales linearly with resolution. This gives us confidence that, at least away from the inclusion, the code is giving the right answer. This kind of result, where the solution is bad close to the surface, but good otherwise, is typical for numerical solutions of this problem.

Figure 3: Error in the pressure outside the inclusion along the line [;y=x/2;] for resolutions of 128 × 128 (blue), 256 × 256 (red), and 512 × 512 (black). The inclusion has a radius [;r_i=0.1;].

Figure 4: As in Figure 3, but with the error scaled with [;h;]. So the medium-resolution error is multiplied by 2 and the high-resolution error is multiplied by 4.