Stegotons: solitary waves arising in non-dispersive periodic media

A fascinating new class of solitary waves was discovered in a 2003 paper by R.J. LeVeque and D. Yong. Solitary waves usually appear as solutions of nonlinear dispersive wave equations, like the KdV or NLS equations. But the waves discovered by LeVeque and Yong arise in a system of nonlinear wave equations with no dispersion! The nonlinear elasticity equations they investigated are: \begin{align} \epsilon_t(x,t) - u_x(x,t) & = 0 \\ (\rho(x) u(x,t))_t - \sigma(\epsilon(x,t),x)_x & = 0. \end{align} They took the density $\rho(x)$ and bulk modulus $K(x)$ to be periodic functions and the stress strain relation $\sigma(\epsilon,x)$ nonlinear (they used an exponential function, but any nonlinear function will work).

If $\rho$ and $K$ are chosen so that the impedance $Z = \sqrt{\rho K}$ varies in space, then waves undergo reflection on the fine scale. Remarkably, the effect of these reflections is an effective behavior that mimics dispersion -- as predicted already by Santosa & Symes in 1993.

Here we reproduce some of their original experiments in PyClaw. If you're interested in the simulations, all of the code is provided here and you can run it yourself. If you're only interested in the results, feel free to just skip over the code.

In [1]:
from clawpack import riemann
from clawpack import pyclaw
import numpy as np

riemann_solver = riemann.nonlinear_elasticity_fwave_1D
solver = pyclaw.ClawSolver1D(riemann_solver)
solver.fwave = True

# Boundary conditions
solver.bc_lower[0] = pyclaw.BC.extrap
solver.bc_upper[0] = pyclaw.BC.extrap
solver.aux_bc_lower[0] = pyclaw.BC.extrap
solver.aux_bc_upper[0] = pyclaw.BC.extrap

xlower=0.0; xupper=1000.0
cells_per_layer=12; mx=int(round(xupper-xlower))*cells_per_layer
x = pyclaw.Dimension('x',xlower,xupper,mx)
domain = pyclaw.Domain(x)
state = pyclaw.State(domain,solver.num_eqn,3)

#Initialize q and aux
KA    = 1.0; rhoA  = 1.0
KB    = 4.0; rhoB  = 4.0
xfrac = xc-np.floor(xc)

state.aux[0,:] = rhoA*(xfrac<0.5)+rhoB*(xfrac>=0.5) #Density
state.aux[1,:] = KA  *(xfrac<0.5)+KB  *(xfrac>=0.5) #Bulk modulus
state.aux[2,:] = 0. # not used

sigma = 0.5*np.exp(-((xc-500.)/5.)**2.)
state.q[0,:] = np.log(sigma+1.)/state.aux[1,:]  # Strain
state.q[1,:] = 0.                               # Momentum

claw = pyclaw.Controller()
claw.solution = pyclaw.Solution(state,domain)

claw.output_style = 1
claw.num_output_times = 100
claw.tfinal =  550.
claw.solver = solver
claw.keep_copy = True
claw.output_format = None

Before running the simulation, let's take a quick look at the setup. The initial condition (stored in state.q) is a Gaussian stress perturbation with zero velocity, while the impedance $Z(x)$ is piecewise constant and periodic. Notice in the plot below that, though the stress is a Gaussian, the strain is discontinuous since the bulk modulus $K(x)$ is discontinuous.

In each plot, the inset is a closeup of the central region.

In [2]:
import plotly
py = plotly.plotly(username_or_email='DavidKetcheson', key='mgs2lgb203')

data = [{'x':xc, 'y':state.q[0,:]},{'x':xc[5800:6200],'y':state.q[0,5800:6200],"xaxis":"x2","yaxis":"y2"}]
layout = {"title": "$\\text{Strain }(\epsilon)$",'showlegend':False,    "xaxis2": {"domain": [0.6, 0.95],"anchor": "y2"},
    "yaxis2":{"domain": [0.2, 0.6],"anchor": "x2"}}
In [3]:
Z = np.sqrt(state.aux[0,:]*state.aux[1,:])
data = [{'x':xc, 'y':Z},{'x':xc[5800:6200],'y':Z[5800:6200],"xaxis":"x2","yaxis":"y2"}]
layout = {"title": "$\\text{Impedance }(Z)$",'showlegend':False,    "xaxis2": {"domain": [0.6, 0.95],"anchor": "y2"},
    "yaxis2":{"domain": [0.2, 0.6],"anchor": "x2"}}
py.iplot(data,layout = layout)
In [4]:
{'cflmax': 0.9014990065747429,
 'dtmax': 0.06465763584693647,
 'dtmin': 0.06436880097776293,
 'numsteps': 86}
In [5]:
%pylab inline
from matplotlib import animation
import matplotlib.pyplot as plt
from clawpack.visclaw.JSAnimation import IPython_display

fig = plt.figure(figsize=[8,4])
ax = plt.axes(xlim=(xc[0], xc[-1]), ylim=(0, 0.4))
line, = ax.plot([], [], lw=1)

def fplot(i):
    frame = claw.frames[i]
    strain = frame.q[0,:]
    line.set_data(xc, strain)
    ax.set_title('Strain at t='+str(frame.t))
    return line,

animation.FuncAnimation(fig, fplot, frames=len(claw.frames))
Populating the interactive namespace from numpy and matplotlib

Once Loop Reflect