Positive disclination in a thin elastic sheet with boundary

We consider a positive disclination of strength $s$ located at a point $o$ within the 2D simply-connected plate domain $\omega$ with a piecewise smooth boundary $\partial\omega$. The equilibrium equations for a Föppl-von Kármán plate, in the absence of body forces, are written in terms of the in-plane stress $\boldsymbol{\sigma}$ (with components $\sigma_{\alpha\beta}$) and moment $\boldsymbol{m}$ (with components $m_{\alpha\beta}$) tensors, both symmetric, as [14]

where w represents the transverse displacement field of the plate. The stress and moment components are assumed to be related to the stretching and bending strains (given by $\varepsilon_{\alpha\beta}$ and $\text{w}_{,\alpha\beta}$, respectively) through the isotropic, materially uniform, linear elastic constitutive relations [14]

where $E$ is the Young’s modulus of the 2D sheet, $D$ is the bending modulus, and $\nu$ is the Poisson’s ratio. The equilibrium equation (1)₁ is identically satisfied if the stress components are expressed in terms of the scalar Airy stress function $\Phi$ such that $\sigma_{\alpha\beta}=e_{\alpha\mu}e_{\beta\nu}\Phi_{,\mu\nu}$. Due to the wedge disclination at $o\in\omega$, the strain field $\boldsymbol{\varepsilon}$ (with components $\varepsilon_{\alpha\beta}$) cannot be written in terms of a single-valued displacement field. More precisely, we can write $\varepsilon_{\alpha\beta}=\gamma_{\alpha\beta}+\frac{1}{2}\text{w}_{,\alpha}\text{w}_{,\beta}$, where $\gamma_{\alpha\beta}$ are the components of the infinitesimal stretching strain tensor which is not expressible in terms of a well defined in-plane displacement field [9, 6]. This incompatibility of the strain field is expressed in terms of the following equation [7, 6]:

where $\delta_{o}$ is the Dirac measure supported at point $o$. Note that $\frac{1}{2}\left[\text{w},\text{w}\right]=\text{det}(\text{w}_{,\alpha\beta})$ is the Gaussian curvature of the deformed plate. If $s=0$, i.e., there is no disclination, then there exists a single-valued in-plane displacement field $\text{u}_{\alpha}$ such that $\varepsilon_{\alpha\beta}=\frac{1}{2}(\text{u}_{\alpha,\beta}+\text{u}_{\beta,\alpha})+\frac{1}{2}\text{w}_{,\alpha}\text{w}_{,\beta}$; the converse is also true. The incompatibility associated with the positive disclination can be understood in terms of removal of wedges (of net angle $s$) centred at $o$, see Figure 1. Consider, for instance, removing of four symmetrical wedges (of angle $s/4$) from a square domain of side length $L_{0}$. The domain $\omega$ is then given by the square plate of side length $L=L_{0}(1-\tan(s/8))$ obtained from the plate of size $L_{0}$ after removing the four wedges and identifying each pair of the newly exposed edges with a straight line.

The Föppl-von Kármán equations with a disclination can be obtained by writing strain and moment components in terms of stress (and hence stress function) and w, using (2a) and (2b), respectively, and then substituting them into Equations (3) and (1)₂. We obtain

in $\omega$. These equations, combined with appropriate boundary conditions, can be used to obtain the stress field and the deformed shape of the plate due to the presence of a disclination of strength $s$ at $o\in\omega$. More detailed derivations of these equations are available elsewhere [7, 6]. In writing (4), and in interpreting Gaussian curvature as $\text{det}(\text{w}_{,\alpha\beta})$, we assume both $\Phi$ and w to be sufficiently smooth over $\omega$. Strictly speaking, this is overly restrictive and one alternative is to interpret these equations in the sense of distributions. This is however not immediate due to the nonlinear terms in the equations. In Appendix A.2 we have given the assumptions on $\Phi$ and w such that both (4) and the Gaussian curvature can be interpreted reasonably in a distributional form.

Considering a length parameter $L$ (representing the size of the finite plate), we can write the non-dimensional form of Föppl-von Kármán equations (4) as

where $\hat{\Phi}=\Phi/D$, $\hat{\text{w}}=\text{w}/L$, $\hat{x}_{i}=x_{i}/L$, and $\Gamma={EL^{2}}/{D}$. The dimensionless constant $\Gamma$ is known as the Föppl-von Kármán number [13]. The condition for inextensibility is $\boldsymbol{\varepsilon}=\boldsymbol{0}$ in $\omega$, which corresponds to $(1/\Gamma)\Delta^{2}\hat{\Phi}=0$. We note that $\Gamma\to\infty$, or $E\to\infty$ with $D$ and $L$ fixed, does not necessarily imply $(1/\Gamma)\Delta^{2}\hat{\Phi}\to 0$. In fact, as we demonstrate through our analytical and numerical results, the Gaussian curvature remains necessarily negative in finite sub-regions of $\omega$ even as $\Gamma\to\infty$. This in conjunction with (5a) requires $(1/\Gamma)\Delta^{2}\hat{\Phi}\neq 0$ away from $o$. Therefore we have a situation where $\Gamma\to\infty$ does not lead to inextensibility.

We state three types of boundary conditions that are most commonly used with Equations (4) to yield a well posed boundary value problem [15]. All of these can be derived as part of the stationarity conditions from the functional (17) with appropriate choice of test functions. The free boundary condition requires the plate edges to be free of forces and moments, i.e.,

on $\partial\omega$, where $\boldsymbol{t}$ is the unit tangent and $\boldsymbol{n}$ is the in-plane unit normal to the boundary. Whereas the first two conditions enforce that there are no net in-plane forces applied at any point of the boundary, the latter two ensure that there is no moment (about $\boldsymbol{t}$) and no transverse shear force, respectively, being applied at any point of the boundary. The simply supported boundary condition requires the in-plane traction, the moment about the edge tangent, and the out-of-plane displacement to all vanish at the plate boundaries, i.e.,

on $\partial\omega$. In the clamped boundary condition, the plate boundaries are free of in-plane traction and are clamped with respect to the out-of-plane displacement, i.e.,

on $\partial\omega$. The three boundary conditions are illustrated in Figure 2. We will be discussing the solution to the disclination problem (4) subjected to either (6), (7), or (8). It is evident that if w is a solution to any of these boundary value problems, then so is $-\text{w}$. The problems are, in general, analytically intractable and have to be attended numerically. There are however two scenarios, as discussed next, when we are able to obtain exact closed form solutions.

For any of the three boundary value problems stated above, the flat plate solution, with $\text{w}=0$ in $\omega$, always holds true. All three problems are reduced to

whose unique solution for a circular plate of radius $R$ is [16]

with the corresponding stress field

where $\boldsymbol{e}_{r}$ and $\boldsymbol{e}_{\theta}$ are the orthonormal basis vectors in the polar coordinate system $(r,\theta)$. The flat plate solution is not well defined for an unbounded plate. The solution is, in any case, unstable beyond a critical value of $R$ (for fixed $Es/D$) giving way to buckled solutions with $\text{w}\neq 0$ [6]. In this article we will always be working with parametric values where the buckled solution is the stable solution.

The simplest buckled solution is obtained assuming the plate to be elastically inextensible and with an unbounded domain. The former is tantamount to a priori imposing $\boldsymbol{\varepsilon}=\boldsymbol{0}$ in $\omega$. Substituting this into (3), we can derive the reformulated governing equations as

with the stress and moment fields vanishing identically as $r\to\infty$. The stress field, and hence $\Phi$, appears here as a Lagrange multiplier associated with the inextensibility constraint. The minimum energy solution to this problem is given by $\text{w}=\sqrt{\frac{s}{\pi}}r$, which represents a perfect cone, and $\Phi=-D\ln r$, whence we calculate

A rigorous verification of the claim, that $\text{w}=\sqrt{\frac{s}{\pi}}r$ and $\Phi=-D\ln r$ indeed solves the problem at hand, is not straightforward. We use distribution theory to establish the result in Appendix A.2, wherein we also derive stress field (13) from the stress function. Note that both stress and Gaussian curvature fields develop a Dirac singularity at the location of the defect in the plate. This should be compared with (11), where the stress is unbounded at $o$ but has no Dirac singularity. More importantly, the stress field in (13) is independent of the defect strength $s$. For $s=0$, and considering $\text{w}=0$ as a solution, any stress function field (including $\Phi=-D\ln r$) which yields a stress field vanishing at infinity is a solution. With $s\neq 0$, irrespective of the magnitude of $s$, the extent of non-uniqueness in stress is significantly reduced. We show, in Appendix A.3, that given $\text{w}=\sqrt{\frac{s}{\pi}}r$ the most general form of the solution for the stress function is $\Phi=-D\ln r+g_{0}(\theta)+rg_{1}(\theta)$, where $g_{0}(\theta)$ and $g_{1}(\theta)$ are arbitrary periodic functions (with period $2\pi$) which satisfy $\int_{0}^{2\pi}{g_{0}^{\prime}}\boldsymbol{e}_{\theta}\text{d}\theta=\boldsymbol{0}$ and $\int_{0}^{2\pi}{g_{1}}\text{d}\theta=0$. The corresponding non-uniqueness in the stress solution is given in Equation (51) in Appendix (A.3). Due to the inextensibility constraint the variations in the expression for the stress field have no bearing on the stored energy of the plate and hence all the solutions, with fixed w, are energetically equivalent.

A somewhat surprising result is that if we consider the inextensional equations (12) for a bounded plate, with boundary conditions of any form given in Section 2.2, then the ensuing boundary value problem has no solution. To establish this, we start by writing a loop condition

where $\boldsymbol{e}_{3}$ is the unit vector such that $\{\boldsymbol{e}_{r},\boldsymbol{e}_{\theta},\boldsymbol{e}_{3}\}$ form a right-handed orthonormal triad in $\mathbb{R}^{3}$ and $\text{d}L$ is an infinitesimal line element in $\omega$. Equation (14) can be derived by integrating Equation (37) from Appendix A.2, which is the distributional counterpart of (12)₁, over the plate domain $\omega$ and using the Stokes’ theorem. The loop condition in fact holds for any arbitrary loop enclosing the defect point $o$. Under certain regularity conditions on w, the loop condition over arbitrary loops is equivalent to (12)₁. According to (14), $\nabla^{2}\text{w}$ cannot vanish everywhere on $\partial\omega$. We first consider free and simply supported boundary conditions. The boundary condition $\langle\boldsymbol{m},\boldsymbol{n}\otimes\boldsymbol{n}\rangle=0$, on using the constitutive relationship, yields

Using this we can calculate the Gaussian curvature on the boundary as

Therefore, with $\langle\boldsymbol{m},\boldsymbol{n}\otimes\boldsymbol{n}\rangle=0$ on $\partial\omega$, $\nabla^{2}\text{w}\neq\boldsymbol{0}$ implies $[\text{w},\text{w}]\neq 0$. Since $\nabla^{2}\text{w}$ cannot vanish everywhere on $\partial\omega$, the same would follow for $[\text{w},\text{w}]$. This is inconsistent with (12)₁ which requires $[\text{w},\text{w}]=0$ at each point in $\omega-o$. In the case of clamped boundary condition, we have $\text{w}=0$ and $\langle\nabla\text{w},\boldsymbol{n}\rangle=0$ on $\partial\omega$, which together imply that $\nabla\text{w}=\boldsymbol{0}$ on $\partial\omega$. This, however, would trivialise the loop integral (14) and hence render it unequal to the right-hand side constant term.

We solve the boundary value problems using a finite element methodology. We have developed our own code based on a mixed variational principle, according to which the governing equations appear as the stationary conditions of the functional [17, p. 165]

where $\text{d}A$ is an infinitesimal area measure on $\omega$. The square plate domain is discretised using non-conforming C¹-continuous rectangular elements and the weak form of the variational principle is used to obtain a system of nonlinear algebraic equations. The algebraic equations are solved using an arc-length method which is able to trace the nonlinear equilibrium path through the limit point (including snap-back and snap-through). We note that the equations are nonlinear and hence the solutions obtained are not unique. Different solution paths can be traced depending on the initial guess of the parameters involved in the numerical procedure. All the solutions are stationary points of the functional $\Pi$ but not all are necessarily stable. The stable (metastable) solution corresponds to a point of global (local) minima in the strain energy landscape.

We state our results in terms of arbitrarily prescribed length ($l$) and force ($f$) units. The side length of the square plate $L$ and the deformation w both have units as $l$. The Gaussian curvature has a unit of $l^{-2}$. The constitutive parameters $E$ and $D$ have units of $l^{-1}f$ and $lf$, respectively. The stress and the stress function have units of $l^{-1}f$ and $lf$, respectively.

We use the numerical framework to solve a typical problem. We will use the results to motivate the central concerns of this work. We consider a square plate with free boundary condition and a positive disclination of strength $s=\pi/3$ at the centre of the plate. There are no external loads acting on the plate. We take $L=2$, $E/D=8000$, $D=0.01$, $\nu=0.3$, and a mesh of $48\times 48$ square elements. The plate axes are denoted as X and Y (both taking values from the interval $[0,2]$) with origin at one corner. The simulation results are given in Figure 3. In solving for w we fix three corners of the plate to avoid any rigid body motions. The plate appears to deform into a conical shape with a rounded vertex [18]. The smoothening of the cone tip is due to extensional elasticity; the fourth-order derivative term ($\frac{1}{E}\Delta^{2}\Phi$) acts as a regulariser for the nonlinear Monge-Ampere bracket term. Both the Gaussian curvature field and the scaled biharmonic of stress function ($\frac{1}{E}\Delta^{2}\Phi$) show singular behaviour at the defect location. However, unlike the inextensional case, it is not clear how the Dirac singularity in Equation (4a) is distributed between the two terms. The biharmonic plot also reveals an interesting cusp-like feature with the function decreasing sharply to a negative value, as one moves away from the defect, before rising again to a near zero magnitude. This feature is neither a numerical artefact nor a consequence of the boundary conditions, as has been checked rigorously through numerical experiments. The stresses again are singular but whether they have a Dirac singularity, or not, is unclear. The behaviour of the deformation and Gaussian curvature, away from the defect, seems uninteresting from the plots in Figure 3. This is however not so. Indeed, a simple conical solution for w away from the defect will not work at the boundary. It will violate all three sets of boundary conditions mentioned in Section 2.2.

Motivated by the discussion so far, we enumerate the questions that will be addressed in the rest of this article:

1. 1.

   What is the nature of solution close to the defect? More precisely, a) whether the Gaussian curvature field and the stress fields have a Dirac singularity at the defect location?, b) how the Dirac source term in (4a) is shared between the biharmonic and the Gaussian curvature terms?, c) are solution fields, in the close vicinity of the defect, invariant with respect to the type of boundary conditions considered? d) how do these solutions behave as $E/D$ is increased for a fixed $L$?

2. 2.

   What is the nature of solution away from the defect? We study this question with an emphasis on the behaviour of the Gaussian curvature field away from the defect location. In particular, a) how the field behaves for varying $E/D$ (keeping $L$ fixed) and varying plate sizes (keeping $E/D$ fixed)? and b) how the three boundary conditions affect the Gaussian curvature field away from the defect point?

3. 3.

   To what extent buckling is dependent on the three boundary conditions? In this we extend the previous work of Mitchell and Head [19] and Seung and Nelson [6].

In the rest of the paper the domain $\omega$ is taken to be a square plate of side length $L$ with the disclination of strength $s$ located at its centre (position denoted as $o$). We will fix $D=0.01$, $s=\pi/3$, and $\nu=0.3$, unless stated otherwise.

We assume both $\frac{1}{E}\Delta^{2}\Phi$ and $\frac{1}{2}\left[\text{w},\text{w}\right]$ to be measures like $\mu$, i.e.,

where $G_{1}$ and $G_{2}$ are integrable functions and $a_{1},a_{2}$ are constants. In other words, we posit both the scaled biharmonic term and the Gaussian curvature to be given in terms of an integrable function (possibly unbounded at $o$) and a Dirac concentration. Their sum, as it appears in (4a), is equal to $s\delta_{o}$. Consequently, $G_{1}=-G_{2}$ and $a_{1}+a_{2}=s$. The former can be proved by integrating (4a) over an arbitrary $\Omega\subset\omega$ with $o\not\in\Omega$. The latter result can then be established by integrating (4a) over any arbitrary $\Omega\subset\omega$ with $o\in\Omega$. We determine the values of $a_{1}$ and $a_{2}$ using a series of numerical experiments where, for definiteness, we take ${E}/{D}=8000,~{}L=2$, and the free boundary condition. We choose the mesh element containing $o$ as $\Omega_{n}$. For a sequence of mesh refinements we plot the variations in $\frac{1}{E}\Delta^{2}\Phi$ and $\frac{1}{2}[\text{w},\text{w}]$ at a section of the plate containing $o$, see Figures 4(a) and 4(b). In writing a mesh size as $2/24$, for instance, we refer to the case of discretising the plate domain of side length $L=2$ into $24\times 24$ elements. With increasing mesh refinement we expect $\int_{\Omega_{n}}\text{d}A\to 0$. For each instance of the mesh refinement we calculate two numbers: $V_{\frac{1}{E}\Delta^{2}\Phi}^{n}=\int_{\Omega_{n}}\frac{1}{E}\Delta^{2}\Phi\text{d}A$ and $V_{\frac{1}{2}[\text{w},\text{w}]}^{n}=\int_{\Omega_{n}}\frac{1}{2}[\text{w},\text{w}]\text{d}A$. We observe from Figure 4(c) that the former tends to $s$, while the latter tends to $0$, with increasing mesh refinement. This suggests that $a_{1}=s$ and $a_{2}=0$. Such a conclusion remains invariant irrespective of the choice of parameter values (as long as they remain finite) and boundary conditions, as has been verified through several numerical simulations. The term $\frac{1}{E}\Delta^{2}\phi$ therefore takes the whole of Dirac singularity. The Gaussian curvature at $o$ is unbounded but it does not have a Dirac concentration. This is contrary to what we observed in the inextensible case. The elastic extensibility of the plate, no matter how weak, alters the behaviour of the Gaussian curvature field at the defect location. Our result also explains the presence of the cusp like feature in the $\frac{1}{E}\Delta^{2}\Phi$ plots. Indeed, with $a_{1}=s$, $a_{2}=0$, and $G_{1}=-G_{2}$ in (21a), these plots can be interpreted in terms of a superposition of a Dirac onto the negative of the Gaussian curvature distribution.

For establishing that the solution close to $o$ is not significantly affected by our choice of the boundary condition as well as the plate size, we introduce an error

for a given size ($L>1$) and boundary condition, where $\mathcal{R}$ is a domain centred around $o$ of a size less than that of a unit square, $\text{w}_{1}$ is the deformation field corresponding to a plate of size $L=1$, $\text{w}_{L}$ is the deformation field for a plate of size $L$, $k_{1}$ is the Gaussian curvature field for a plate of size $L=1$, and $k_{L}$ is the Gaussian curvature field for a plate of size $L$. For a chosen boundary condition, and for a fixed region $\mathcal{R}$, error $e$ measures the deviation of the solution for a plate of size $L$ from that for a plate of size $L=1$. The results are reported in Figure 5, where each plot corresponds to a different boundary condition. Within each plot, we have reported errors for four plate sizes ($L=1.33,1.5,1.67,2$) and four choices of domain $\mathcal{R}$. For the latter, we have considered square domains, with centre at $o$, having one mesh element (${\mathcal{R}}_{1}$), nine elements (${\mathcal{R}}_{2}$), sixteen elements (${\mathcal{R}}_{3}$), and twenty-five elements (${\mathcal{R}}_{4}$). The error values are low for every case considered. The solution in small regions enclosing the defect therefore does not vary much for different plate sizes and boundary conditions. Moreover, for every boundary condition and plate size, the error values are the lowest when we compute them for solutions in the smallest neighbourhood ${\mathcal{R}}_{1}$ of $o$, and increasing only slightly as we move towards larger domain sizes of ${\mathcal{R}}$. The solution close to the defect therefore changes only minimally as we compare it for various boundary conditions and plate sizes.

We plot $\frac{1}{E}\Delta^{2}\Phi$ and $\frac{1}{2}[\text{w},\text{w}]$, at a section of the plate containing $o$, for increasing values of $E/D$, see Figures 6(a) and 6(b). We consider a fixed $\Omega_{0}$ containing $o$ and evaluate the integrals $V_{\frac{1}{E}\Delta^{2}\Phi}^{0}=\int_{\Omega_{0}}\frac{1}{E}\Delta^{2}\Phi\text{d}A$ and $V_{\frac{1}{2}[\text{w},\text{w}]}^{0}=\int_{\Omega_{0}}\frac{1}{2}[\text{w},\text{w}]\text{d}A$ for various values of $E/D$ while keeping the mesh refinement of $48\times 48$ elements, $L=2$, and the free boundary condition. The large values of the ratio $E/D$ are in fact equivalent to large values of the dimensionless parameter $EL^{2}/D$ for a fixed $L$. We identify $\Omega_{0}$ with the fixed domain of a single element containing $o$. According to Figure 6(c), $V_{\frac{1}{E}\Delta^{2}\Phi}^{0}$ decreases monotonically, possibly towards $0$, and $V_{\frac{1}{2}[\text{w},\text{w}]}^{0}$ increases monotonically, possibly towards $s(=\pi/3)$, as we approach large values of $E/D$. Their sum, as expected, is always close to $s$ in confirmation with the Föppl-von Kármán equation (4a) (the slight but persisting deviation of the sum from $s$ in Figure 6(c) is due to the limited numerical accuracy in calculating the biharmonic of $\Phi$ using finite difference method, especially close to the defect point). This indicates development of a concentration in the Gaussian curvature field. The monotonicity trend persists irrespective of the mesh refinement, plate size, and boundary condition, although with a different rate of convergence. Moreover, for any arbitrary domain (say $\hat{\Omega}$) in the vicinity of $o$ but not containing it, we observe the volume $V_{\frac{1}{2}[\text{w},\text{w}]}^{P}=\int_{\hat{\Omega}}\frac{1}{2}[\text{w},\text{w}]\text{d}A$ to decrease towards zero as we increase $E/D$; the results for one such patch in the form of an annular region (of sixteen elements) are given in Figure 6(d), with the patch shown in the inset. All together, this indicates that the scaled biharmonic term converges to $0$ while the Gaussian curvature converges to a Dirac at $o$ as $E/D\to\infty$. Indeed, a sequence of measures $f_{n}$ (of the type $\mu$) converges to $s\delta_{o}$ if, for any arbitrary open subset $\Omega\subset\omega$, $\int_{\Omega}f_{n}\text{d}A\to s\xi$, where $\xi=1$ if $o\in\Omega$ and $\xi=0$ otherwise. One should keep in mind that, as discussed in Section 3.3, the inextensible problem with boundary has no solution with Gaussian curvature field given only in terms of a Dirac $o$. Our results should not be seen contradictory, to those discussed in Section 3.3, for we are dealing with the solution only in the neighbourhood of the defect. As we shall see in a following section, the value of the Gaussian curvature, away from the defect closer to the boundary, indeed does not become vanishing small even for large values of $E/D$.

We now combine the arguments presented in the last two paragraphs. We showed that, for any finite $E$ (keeping $D$ and $L$ fixed), the Gaussian curvature behaves like an integrable function $G_{2}$ and the scaled biharmonic $\frac{1}{E}\Delta^{2}\Phi$ behaves like $G_{1}+s\delta_{o}$, in a neighbourhood of $o$, with $G_{2}=-G_{1}$. As $E\to\infty$, $G_{2}\to s\delta_{o}$ and $\frac{1}{E}\Delta^{2}\Phi\to 0$. However, as shown in Appendix A.4, $\Delta^{2}\Phi\to c\Delta\delta_{o}$, where $c\in\mathbb{R}$ is a constant. Such a behaviour of $\Delta^{2}\Phi$ would follow from a stress field which has a Dirac concentration at $o$. The latter is indeed the case, as verified numerically in the following section. The limiting behaviour of both the Gaussian curvature and the stress are in line with the solution for the unbounded plate with inextensional constraint (as obtained in Section 3.2). The corresponding solution in a bounded plate hence retains the essential aspects of the infinite plate solution close to the defect.

We now study the singular nature of the stress field around the defect. The stress distribution is observed to remain invariant with respect to the choice of the boundary conditions and plate size, while keeping all other parameters fixed. We establish the nature of singularity in the stress field for a fixed $E/D$. Following the framework developed in the preceding section, we assume all the three Cartesian components of the stress ($\sigma_{11}$, $\sigma_{22}$, and $\sigma_{12}$) to be measures like $\mu$, i.e., each of them is given in terms of an integrable (possibly unbounded) function and a Dirac concentration. As before, we take ${E}/{D}=8000,~{}L=2$, the free boundary condition, and choose the smallest mesh element containing $o$ as $\Omega_{n}$. For a sequence of mesh refinements we plot the variations in the stress components $\sigma_{\alpha\beta}$ at a section of the plate containing $o$, see Figure 7. For each instance of the mesh refinement we calculate three numbers: $V_{\sigma_{11}}^{n}=\int_{\Omega_{n}}\sigma_{11}\text{d}A$, $V_{\sigma_{22}}^{n}=\int_{\Omega_{n}}\sigma_{22}\text{d}A$, and $V_{\sigma_{12}}^{n}=\int_{\Omega_{n}}\sigma_{12}\text{d}A$. We observe that $V_{\sigma_{11}}^{n}=V_{\sigma_{22}}^{n}$ and $V_{\sigma_{12}}^{n}=0$, irrespective of the mesh element size. Moreover, with increasing mesh refinement $V_{\sigma_{11}}^{n}=V_{\sigma_{22}}^{n}$ tend towards $0$. Therefore, like the Gaussian curvature field, the stresses are unbounded at $o$ but without developing a Dirac concentration. This is again contrary to what was derived in the inextensible case (for an unbounded plate).

We consider a fixed $\Omega_{0}$ containing $o$ and evaluate the integrals for increasing $E/D$ values (for fixed $L$). The results are given in Figure 8, considering the free boundary condition and $\Omega_{0}$ as the single element centred at $o$ in a plate domain discretised with $48\times 48$ square elements. Whereas $\int_{\Omega_{0}}\sigma_{12}\text{d}A=0$, for all values of $E/D$, the integrals $V_{\sigma_{11}}^{0}=\int_{\Omega_{0}}\sigma_{11}\text{d}A$ and $V_{\sigma_{22}}^{0}=\int_{\Omega_{0}}\sigma_{22}\text{d}A$ are both equal and increase (in magnitude) with increasing $E/D$. There is always a bulk contribution to the integrals, as is clearly evident from the $\sigma_{22}$ plots in Figure 8. The limiting value of the integrals over an arbitrary open set in $\omega$, containing $o$, will therefore have contributions from the limiting concentration at $o$ and the limiting non-zero bulk value. If we conjecture that this limiting value of the stress is of the form given in (13) then, clearly, the limiting bulk field for stress is non-integrable and hence not well defined for an arbitrary measurable set. We can resolve this problem by interpreting the integrals as $\int_{\Omega}\sigma_{\alpha\beta}\text{d}A=\lim_{\epsilon\to 0}\int_{\Omega-B_{\epsilon}}\sigma_{\alpha\beta}\text{d}A$ for any open set $\Omega$ containing $o$, where $B_{\epsilon}$ is an open disc of radius $\epsilon$ centred at $o$.

We begin by determining the sign of the Gaussian curvature at the boundary points of the plate domain for various boundary conditions. We assume that the curvature $\nabla^{2}\text{w}$ remains non-zero in the considered regions. This is reasonable since we do not expect the solution to deviate far from the perfect cone. Recall, for the free boundary condition, that we require $\langle\boldsymbol{m},\boldsymbol{n}\otimes\boldsymbol{n}\rangle=0$ for all points on $\partial\omega$, which on using the constitutive relation can be rewritten as $\langle\nabla^{2}\text{w},\boldsymbol{n}\otimes\boldsymbol{n}\rangle=-\nu\langle\nabla^{2}\text{w},\boldsymbol{t}\otimes\boldsymbol{t}\rangle$, assuming $D\neq 0$. If $\nu>0$ then $\langle\nabla^{2}\text{w},\boldsymbol{n}\otimes\boldsymbol{n}\rangle$ and $\langle\nabla^{2}\text{w},\boldsymbol{t}\otimes\boldsymbol{t}\rangle$ are of opposite sign and if $\nu<0$ (auxetic materials) then they are of same sign on $\partial\omega$. The Gaussian curvature $\frac{1}{2}[\text{w},\text{w}]=\langle\nabla^{2}\text{w},\boldsymbol{n}\otimes\boldsymbol{n}\rangle\langle\nabla^{2}\text{w},\boldsymbol{t}\otimes\boldsymbol{t}\rangle-\langle\nabla^{2}\text{w},\boldsymbol{t}\otimes\boldsymbol{n}\rangle^{2}$ on $\partial\omega$ is therefore negative for plates with $\nu>0$ while its sign is undecided when $\nu<0$. If $\nu=0$ then $\langle\nabla^{2}\text{w},\boldsymbol{n}\otimes\boldsymbol{n}\rangle=0$ yielding a negative Gaussian curvature on $\partial\omega$. For the simply supported boundary condition, $\text{w}=0$ and $\langle\boldsymbol{m},\boldsymbol{n}\otimes\boldsymbol{n}\rangle=0$ on $\partial\omega$. If in addition the boundary is piecewise straight, as is the case for the square plate, we have $\nabla\boldsymbol{t}=\boldsymbol{0}$ and $\nabla\boldsymbol{n}=\boldsymbol{0}$ almost everywhere on $\partial\omega$. Under this simplification, the boundary condition $\text{w}=0$ can be differentiated twice to yield $\langle\nabla^{2}\text{w},\boldsymbol{t}\otimes\boldsymbol{t}\rangle=0$ which, on using the other boundary condition, gives $\langle\nabla^{2}\text{w},\boldsymbol{n}\otimes\boldsymbol{n}\rangle=0$ almost everywhere on $\partial\omega$, regardless of the value of $\nu$. Therefore, the Gaussian curvature for almost all the boundary points of a simply supported square plate is necessarily negative. For the clamped boundary condition, $\text{w}=0$ and $\langle\nabla\text{w},\boldsymbol{n}\rangle=0$ on $\partial\omega$. For a square plate these conditions yield $\langle\nabla^{2}\text{w},\boldsymbol{t}\otimes\boldsymbol{t}\rangle=0$ and $\langle\nabla^{2}\text{w},\boldsymbol{t}\otimes\boldsymbol{n}\rangle=0$ almost everywhere on $\partial\omega$. Consequently the Gaussian curvature is identically zero at almost all boundary points of the square plate with a clamped boundary. Following similar arguments we can show that the derivative of the Gaussian curvature, along $\boldsymbol{n}$, also vanishes at almost all boundary points for the square plate with clamped boundary condition. The results about the sign of Gaussian curvature for free and simply supported boundary conditions, and those about vanishing of the same (and its slope) for the clamped boundary condition, are in agreement with the numerical solutions in Figures 9(a) and 9(b).

We can obtain some further analytical understanding of the nature of the Gaussian curvature, and its slope, at the boundary if we restrict our attention to a circular plate (of radius $r_{0}$) with the free boundary condition. This allows us to consider a smooth axisymmetric solution for w, away from the defect, of the form $\text{w}=f(r)$. Hence $\nabla^{2}\text{w}=f^{\prime\prime}\boldsymbol{e}_{r}\otimes\boldsymbol{e}_{r}+({f^{\prime}}/{r})\boldsymbol{e}_{\theta}\otimes\boldsymbol{e}_{\theta}$, where the superscript prime denotes the derivative with respect to $r$. The Gaussian curvature and its slope along the radial direction can then be written as

respectively. The moment tensor takes the form

The boundary condition (6)₃ then implies $f^{\prime\prime}=-\nu({f^{\prime}}/{r_{0}})$ whereas (6)₄ yields $f^{\prime\prime\prime}+({f^{\prime\prime}}/{r_{0}})-({f^{\prime}}/{{r_{0}}^{2}})=0$. It is reasonable to assume that $f(r)$ is close to a cone like solution in the sense that $f^{\prime}\approx\sqrt{s/\pi}$, as is clear from our numerical simulations. Consequently, $f^{\prime\prime}\approx-\sqrt{s/\pi}({\nu}/{r_{0}})$ and $f^{\prime\prime\prime}\approx\sqrt{s/\pi}(1+\nu)/{r_{0}}^{2}$. Substituting these into (23) we obtain

The Gaussian curvature and its slope at the boundary therefore scale as $-({1}/{{r_{0}}^{2}})$ and $({1}/{{r_{0}}^{3}})$, respectively, with the size of the domain. Taking the effective radius of the square plate with side $L$ as $r_{0}=L/\sqrt{\pi}$, we superpose the analytically predicted behaviour of the end-values with those obtained from numerical solutions in Figure 9(c). The two solutions are in very good agreement except for the slope value at $L=1$.

The Gaussian curvature field oscillates as it moves away from the defect towards the boundary, irrespective of the material parameters, plate size, and the type of boundary condition, see Figures 9(a) and 9(b). In doing so, the curvature values become negative over large regions in the plate, see Figure 10. The existence of non-positive Gaussian curvature values at the boundary has been argued analytically in the preceding discussion for all the boundary conditions. In fact, as we demonstrate below, the average Gaussian curvature (over the plate) is necessarily zero both for the clamped and the simply supported case irrespective of the material and geometric parameters. Therefore the sheet will necessarily have regions of negative Gaussian curvature, large enough in their extent and in the magnitude of the curvature, so as to balance the substantial positive Gaussian curvature in the neighbourhood of the defect. Using identity (31) from Appendix A, we can write

For a clamped boundary $\nabla\text{w}=\boldsymbol{0}$ on $\partial\omega$. The net Gaussian curvature over $\omega$, given by the left-hand side of (26), is therefore zero. This is strictly a topological requirement for sheets with clamped boundary condition over boundaries of arbitrary shape. Indeed, clamping the sheet forces the integrated geodesic curvature on the boundary to be $2\pi$ which, on using the Gauss-Bonnet theorem, implies the vanishing of the net Gaussian curvature in the domain. For a simply supported boundary we can reach the same conclusion for piece-wise linear boundaries (like that of a square plate). Over the boundary, away from the corner points, $\langle\nabla^{2}\text{w},\boldsymbol{t}\otimes\boldsymbol{t}\rangle={0}$ and $\nabla\text{w}=\langle\nabla\text{w},\boldsymbol{n}\rangle\boldsymbol{n}$, as is immediate from (7)₄ and the constancy of $\boldsymbol{t}$. This leads to the vanishing of the integrand on the right-hand side of (26) over the boundary except at the corner points. At a corner point, a non-zero value of $\nabla\text{w}$ would necessarily lead to a non-trivial jump in the value of $\nabla\text{w}$ and therefore to a concentration in $\nabla^{2}\text{w}$, which is energetically unfavourable. Hence $\nabla\text{w}$ will necessarily vanish at the corner points, leading to our assertion. This result, again topological in nature, will hold for any simply supported plate with a polygonal shape. We have verified these claims, for the vanishing of the average Gaussian curvature, from our numerical simulations. Keeping them in mind, and recalling Figure 9(a), it is difficult to argue for the formation of a boundary layer as we move towards large $E/D$ values, contradictory to what one would intuitively expect in such a limiting solution.

The total strain energy $U$ stored in the plate due to a positive disclination is given in terms of stretching and bending energies, $U=U_{s}+U_{b}$, where

respectively. The identity (31) from Appendix A, when used for $\Phi$, yields

For any of the boundary conditions given in Section 2.2, $\nabla\Phi=\boldsymbol{0}$. The contribution from stretching energy therefore is limited to only the first term of the integral in (27a). Similarly, using (26), we note that the second term in the bending energy integral (27b) is identically zero for clamped boundary condition where $\nabla\text{w}=\boldsymbol{0}$ on $\partial\omega$. In fact, the total energy for the clamped problem is independent of $\nu$. Indeed, $\nu$ does not enter either the boundary conditions or the energy expression for a clamped boundary value problem. On the contrary, there is a $\nu$ dependence in the free boundary and the simply supported boundary problems through the boundary conditions as well as the second term in the bending energy integral (27b). For the flat solution, irrespective of the boundary condition, $U=U_{s}=({1}/{2E})\int_{\omega}(\Delta\Phi)^{2}\text{d}A$ with $\Phi$ determined from solving the system of Equations (9). For a circular plate of radius $R$, the total energy for the flat solution is ${Es^{2}R^{2}}/{32\pi}$ (which increases unboundedly with the size of the plate). The flat solution remains the stable solution to our problem prior to the buckling transition to the non-flat (buckled, $\text{w}\neq 0$) solutions [19, 6].

According to Seung and Nelson [6], for a plate with free boundary condition, the buckling transition is given in terms of a dimensionless number $y_{c}=R\sqrt{{Es}/{D}}$, where $y_{c}$ depends only on $\nu$. For a fixed plate size ($R$), disclination strength ($s$), and bending modulus ($D$) this formula can also be used to calculate the critical value of the stretching modulus $E$ and its variations with respect to $\nu$. For our present discussion, we fix $L=2$, $s=\pi/3$, $D=0.01$, and a mesh size of $64\times 64$ elements. The effective radius is calculated as $R=L/\sqrt{\pi}$. The calculated values for the critical $E$ for a range of $\nu$ values (corresponding to stable isotropic elastic materials) and for various choices of boundary conditions are given in Figure 11. The variation of critical $E$ with $\nu$ for the free boundary condition is consistent with the prediction of Seung and Nelson [6]. On the other hand, as expected, the critical $E$ for the clamped boundary does not vary with $\nu$. The trend in the variation of critical $E$ with $\nu$ for the simply supported boundary condition is similar to that for the free boundary, however with higher magnitudes. For any given $\nu$, the critical $E$ is always highest for the clamped boundary and lowest for the free boundary. More importantly, it is clear that the buckling transitions are significantly dependent on the nature of boundary conditions.