# Which numerical method can I use to solve this system of hyperbolic PDEs?

Backround

The mathematical model I am trying to numerically solve models wave propagation inside a cylinder with specific material properties suited for dynamic loading. The cylinder's upper base is under the impact pressure of a rigid body, while the lower base is fixed (as presented in the picture below):

Model

In its dimensionless form, the model consists of the following three hyperbolic PDEs in which the given functions depend on time, as well as the radial and longitudinal coordinates:

$${\partial ^2 u \over \partial r^2} + {1 \over r}{\partial u \over \partial r} +{\partial ^2 u \over \partial z^2}-f={1\over c^2}{\partial ^2 u \over \partial t^2} \tag 1$$

$${\partial ^2 f \over \partial r^2} + {1 \over r}{\partial f \over \partial r} +{\partial ^2 f \over \partial z^2}={\partial ^2 f \over \partial t^2} \tag 2$$

$${\partial ^2 v \over \partial r^2} + {1 \over r}{\partial v \over \partial r} -{1\over r^2}v+c^2{\partial ^2 v \over \partial z^2}+\big(1-c^2\big){\partial ^3 u \over \partial r\partial z^2}={\partial ^2 v \over \partial t^2} \tag 3$$

where $$u$$ is the longitudinal displacement, $$f$$ is a function defined by $$(1)$$, $$v$$ is the longitudinal rate of change of radial displacement and $$c$$ is a material constant.

The initial conditions of these functions are zero:

$$u(r,z,0)=f(r,z,0)=v(r,z,0)=0 \tag 4$$ $${\partial u \over \partial t} \Bigg |_{t=0}={\partial f \over \partial t} \Bigg |_{t=0}={\partial v \over \partial t} \Bigg |_{t=0}=0 \tag 5$$

The cylinder is fixed at its bottom base, so I have the following zero boundary conditions at $$z = 1$$:

$$u(r,1,t)=f(r,1,t)=v(r,1,t)=0 \tag 6$$

Also, I have the following zero boundary conditions at the cylinder's axis:

$${\partial u \over \partial r} \Bigg |_{r=0}={\partial f \over \partial r} \Bigg |_{r=0} = v(0,z,t) = 0 \tag 7$$

Note that instead of equations $$(1)$$ and $$(2)$$ which contain singularities, the following two equations can be used at the cylinder's axis:

$$2{\partial ^2 u \over \partial r^2} + {\partial ^2 u \over \partial z^2}-f={1\over c^2}{\partial ^2 u \over \partial t^2} \tag 8$$

$$2{\partial ^2 f \over \partial r^2} + {\partial ^2 f \over \partial z^2}={\partial ^2 f \over \partial t^2} \tag 9$$

The boundary conditions of the cylinder's curved surface are the following:

$${\partial u \over \partial r} \Bigg |_{r=r_{max}}=-v(r_{max},z,t) \tag {10}$$ $${\partial v \over \partial r} \Bigg |_{r=r_{max}}={2c^2-1 \over 2\big(c^2-1\big)}f(r_{max},z,t) \tag {11}$$

$${\partial f \over \partial r} \Bigg |_{r=r_{max}}= {c^2-1 \over c^2} \Bigg( {\partial ^2 v \over \partial r^2}\Bigg |_{r=r_{max}} + {1 \over r}{2c^2-1 \over 2\big(c^2-1\big)}f(r_{max},z,t) -{1\over r^2}v(r_{max},z,t)-{\partial ^2 v \over \partial z^2}\Bigg |_{r=r_{max}} \Bigg) \tag {12}$$

For the upper base of the cylinder, I have the following two zero boundary conditions:

$${\partial f \over \partial z} \Bigg |_{z=0} = v(r,0,t)=0 \tag {13}$$

As for the boundary condition for the longitudinal displacement, it satisfies the position continuity condition:

$$u(r,0,t) = p(t) \tag {14}$$

where $$p(t)$$ is the position of the rigid body. This function satisfies the following equation (derived by using Newtons and Hooks laws):

$${d^2 p \over dt^2}=-m {\partial u \over \partial z} \Bigg |_{z=0}\tag {15}$$ where $$m$$ represents the elastic body and rigid body mass ratio. The following two initial conditions are known:

$$p(0)=0 \tag {16}$$ $${dp\over dt}\Bigg |_{t=0}=V_0 \tag {17}$$

where the initial rigid body velocity is under the restriction $$V_0<1$$. Also, at $$z=0$$, I have the following boundary conditions:

$${\partial u \over \partial r} \Bigg |_{z=0}=0 \tag {18}$$ $${\partial^2 u \over \partial z^2} \Bigg |_{z=0}={c^2 \over c^2-1}f(r,0,t) \tag {19}$$

The model I wrote down has the following problems. The first is the corner singularity problem at the corner $$(r=r_{max},z=0)$$. Namely, the condition $$(11)$$ and $$v(r_{max},0,t)=0$$ can not be satisfied both at the same time. For my simulation, I chose the condition $$v(r_{max},0,t)=0$$ to be satisfied. However, this is a weak singularity, meaning the solution can still be approximated numerically, but with reduced accuracy. More on this can be found in J.P. Boyd's book "Cheybeshev and Fourier Spectral Methods"".

The other problem is that the rigid body and the elastic body contact surfaces do not have the same velocities at the initial moment of the impact. This is why I would restrict equation $$(1)$$ to be solved for $$z>0$$ and I would calculate the longitudinal displacement at $$z=0$$ by solving equation $$(15)$$ first and then applying the condition $$(14)$$.

I tried to solve these equations using the finite difference method. As expected, it didn't give a stable solution. I used the center difference scheme which has the problem of numerical diffusion. Also, the solution for $$v$$ has a big spike near the corner $$(z=0,r=r_{max})$$ which throws all of the other calculations to spike at the same place out of bounds. This is shown in the pictures below.

However, I also noticed another problem. It takes five time steps for the function $$v$$ to finally be non-zero at $$(r, z=\Delta z)$$, and six time steps for the function $$f$$ to be non-zero at $$(r,z=0)$$. I fear this is physically impossible, especially since the function $$f$$ needs to satisfy the compatibility condition $$(19)$$. At the sixth time step, the function $$u$$ was already non-zero at $$(r, z=\{0,\ \Delta z,\ 2\Delta z,\ 3\Delta z,\ 4\Delta z,\ 5\Delta z\})$$

This means that choosing a numerical method that can actually solve this problem is not trivial.

I have a good background in finite differences and basic knowledge of the spectral method. Apart from that, I don't know how to use some other popular methods such as finite elements. But I would like to study and learn any numerical method, as long as it can solve my model.

So my question is, which numerical method can I use to solve my problem? How can I compensate the corner singularity, numerical diffusion and the compatibility-delay problem I described?

A reference to literature to learn the appropriate method would also be helpful. Thank you for your time.