Deflection of a membrane

Authors: Hans Petter Langtangen and Anders Logg.

Modified for DOLFINx by Jørgen S. Dokken

In the first FEniCSx program, we solved a simple problem which we could easily use to verify the implementation. In this section, we will turn our attentition to a physically more relevant problem with solutions of a somewhat more exciting shape.

We would like to compute the deflection \(D(x,y)\) of a two-dimensional, circular membrane of radius \(R\), subject to a load \(p\) over the membrane. The appropriate PDE model is

(11)\[\begin{align} -T \nabla^2D&=p \quad\text{in }\quad \Omega=\{(x,y)\vert x^2+y^2\leq R^2 \}. \end{align}\]

Here, \(T\) is the tension in the membrane (constant), and \(p\) is the external pressure load. The boundary of the membrane has no deflection. This implies that \(D=0\) is the boundary condition. We model a localized load as a Gaussian function:

(12)\[\begin{align} p(x,y)&=\frac{A}{2\pi\sigma}e^{-\frac{1}{2}\left(\frac{x-x_0}{\sigma}\right)^2-\frac{1}{2}\left(\frac{y-y_0}{\sigma}\right)^2}. \end{align}\]

The parameter \(A\) is the amplitude of the pressure, \((x_0, y_0)\) the location of the maximum point of the load, and \(\sigma\) the “width” of \(p\). We will take the center \((x_0,y_0)\) to be \((0,R_0)\) for some \(0<R_0<R\). Then we have

(13)\[\begin{align} p(x,y)&=\frac{A}{2\pi\sigma}e^{-\frac{1}{2}\left(\left(\frac{x}{\sigma}\right)^2 +\left(\frac{y-R_0}{\sigma}\right)^2\right)}. \end{align}\]

Scaling the equation

There are many physical parameters in this problem, and we can benefit from grouping them by means of scaling. Let us introduce dimensionless coordinates \(\bar{x}=\frac{x}{R}\), \(\bar{y}=\frac{y}{R}\), and a dimensionless deflection \(w=\frac{D}{D_e}\), where \(D_e\) is a characteristic size of the deflection. Introducing \(\bar{R}_0=\frac{R_0}{R}\), we obtain

(14)\[\begin{align} -\frac{\partial^2 w}{\partial \bar{x}^2} -\frac{\partial^2 w}{\partial \bar{y}^2} &=\frac{R^2A}{2\pi\sigma TD_e}e^{-\frac{R^2}{2\sigma^2}\left(\bar{x}^2+(\bar{y}-\bar{R}_0)^2\right)}\\ &=\alpha e^{-\beta^2(\bar{x}^2+(\bar{y}-\bar{R}_0)^2} \end{align}\]

for \(\bar{x}^2+\bar{y}^2<1\) where \(\alpha = \frac{R^2A}{2\pi\sigma TD_e}\) and \(\beta=\frac{R}{\sqrt{2}\sigma}\).

With an appropriate scaling, \(w\) and its derivatives are of size unity, so the left-hand side of the scaled PDE is about unity in size, while the right hand side has \(\alpha\) as its characteristic size. This suggests choosing alpha to be unity, or around unity. In this particular case, we choose \(\alpha=4\). (One can also find the analytical solution in scaled coordinates and show that the maximum deflection \(D(0,0)\) is \(D_e\) if we choose \(\alpha=4\) to determine \(D_e\).) With \(D_e=\frac{AR^2}{8\pi\sigma T}\) and dropping the bars we obtain the scaled problem

(15)\[\begin{align} -\nabla^2 w = 4e^{-\beta^2(x^2+(y-R_0)^2)} \end{align}\]

to be solved over the unit disc with \(w=0\) on the boundary. Now there are only two parameters which vary the dimensionless extent of the pressure, \(\beta\), and the location of the pressure peak, \(R_0\in[0,1]\). As \(\beta\to 0\), the solution will approach the special case \(w=1-x^2-y^2\). Given a computed scaled solution \(w\), the physical deflection can be computed by

(16)\[\begin{align} D=\frac{AR^2}{8\pi\sigma T}w. \end{align}\]