Assessing the Structural Performance of Biodegradable Capsules

Abstract

Biodegradable materials pose challenges over all aspects of computational mechanics. In this study, the focus is on the resulting domain uncertainty. Model structures or devices are shells of revolution subject to random variation of the outer surface. The novelty of the proposed computational approach is the possibility to restrict the variation to specific parts of the structure using a posteriori filtering, which is applied to the random process whose realisations are the profiles of the shells. The dimensionally reduced stochastic elasticity problems are solved using a collocation method where every realisation is discretised separately. The collocation scheme is validated against standard Monte Carlo. The reliability of the simulations is further confirmed via a posteriori error estimates that are computed using the same collocation scheme. The quantities of interest on the nominal domain are the expected displacement fields and their variances.

1. Introduction

Biodegradable materials pose great challenges for computational mechanics. Due to the wide spectrum of different options in both materials and manufacturing, many of the existing studies available are focused on some specific question; see, for instance [1,2]. For an overview of the issues related to modelling problems related to such materials, see [3] and references therein. In this study, the focus is on the effect of loss of material, which is assumed to be measurably greater than what would normally be assumed to result from normal wear over time. The modelling approach adopted here is to take the degradation process as a random process that results in deformation of the object of interest.

One of the promising new developments in long-term drug delivery devices is the introduction of implantable capsules [4], which could be patient-specific and constructed using 3D printing with biodegradable materials. Assuming current manufacturing technologies, these capsules belong into the category of thick shells, a class of structures that falls in between full 3D elasticity and dimensionally reduced thin-shell models. Images of representative samples are shown in Figure 1. Two details of irregular boundaries are included within the set. It should be emphasised that the degradation visible in the images is due to the manufacturing process only. However, it is to be expected that with biodegradable materials, similar features will develop over time.

A full simulation model would have to include the effects of time-dependent chemical processes on the material properties, as well as the deformations that naturally would also evolve in time. The most suitable form for current stochastic finite element methods is to have the random input in the form of Karhunen–Loéve expansion, which should also be time-dependent. Parabolic or time-dependent problems in the context of model problems have been discussed by Gittelson et al. [5]. Problems with multiple, yet independent, spatial random variables in elasticity were recently discussed in [6]. In order to keep the discussion focused on uncertain domains, in this paper, it is assumed that the process is studied at some given time $t_0$ and the material properties are considered isotropic, and linear elasticity remains valid. Naturally, this also reduced the computational complexity significantly.

For standard (made of metals) shells, there have been many studies considering various imperfections, with the Delft Imperfection Data Bank as the prime example. A description of the data bank and a comprehensive review of measurement procedures is given in a thesis by de Vries [7]. To measure imperfections and build related models, first of all, one needs shells. Testing and manufacturing of the kind of capsules that are the focus of this study has only started, and there simply are not enough data to build even rudimentary material models. However, understanding simulation processes is likely to lead to more efficient testing and procedures, and thus speed up the development of realistic material models.

1.1. Novelty of This Study

Formally, the degradation process leads to a stochastic problem with an uncertain domain. This problem is solved using stochastic finite element techniques (SFEM). Within the stochastic collocation process, each realisation of the process is computed separately with exact geometry description. What makes the current study unique is that the geometric deformation process is coupled with a Butterworth-type filter, with which the effect can be localised to act only at given sections of the capsule. This leads to increased flexibility in considering devices with regions with different coatings, for example. Notice that this does not mean that the resulting systems would be partly deterministic and partly stochastic; as always, any uncertainty in the model results in a fully stochastic solution process.

Another novelty of the current study is the application of a posteriori error estimates that are represented on the same discretisations as the actual simulations. This makes it possible to analyse the error estimates in exactly the same way as the solutions, that is, one can compute both expectation and the variance of the error over the parameter space. This approach is shown to result in meaningful and highly understandable descriptions of the simulation errors that further increase the reliability of the interpretations of the simulated solutions.

1.2. Approaches to Domain Uncertainty

Besides material degradation, manufacturing imperfections also lead to problems where the shape of the object may not be perfectly known. Modelling domain uncertainty has been studied by many authors, mostly within model problems. If the domain perturbations are small, application of the perturbation technique as in, e.g., [8,9,10], is straightforward. In the perturbation approach, which is based on Eulerian coordinates, the shape derivative is used to linearise the problem under consideration around a nominal reference geometry. Statistical moments of the problem can be derived using simple first-order shape expansions.

Another approach is the domain-mapping approach. The shape uncertainty is transferred onto a fixed reference domain, and hence reflects the Lagrangian setting. In particular, it allows us to deal with large deformations; see [11,12], for example. Yet another approach, the one adopted in this study, is to use direct mapping, that is, every realisation is discretised separately. This is feasible in the context of high-order finite element methods or the p-version of FEM [13]. For a detailed analysis of this approach in the context of quasi-Monte Carlo methods (QMC), see [14]. In its basic variant, the deformed domain is pulled back to the nominal domain simply by mesh remapping. In complicated domains, conformal mappings in 2D provide the best alternative [15].

1.3. Outline of Paper

The rest of the paper is structured as follows: First, the simple model for domain degradation is introduced (Section 2). Next, the shells of revolution are briefly introduced within the context of this study (Section 3). Section 4 introduced the stochastic model and the finite element method used to perform the simulations. An extensive set of numerical simulations is analysed before the conclusions are drawn. The linear elasticity model is outlined in the Appendix A.

2. Modelling Degradation

Degradation is not simple to model as a random process. If one takes degradation as loss—which is quite reasonable—it has to be modelled as something that cannot change its sign. In other words, degradation is monotone, and in the context of this study, it cannot be reversed.

There are two immediate choices available. The first one is to model the process as a lognormal one. Second, a random field such as Karhunen–Loéve (KL) expansion can be squared. This latter approach is the one adopted here. Some realisations of KL-expansions are shown in Figure 2. Every initial configuration represents a nominal domain. The realisations by definition will have a smaller volume (or area, depending on the model), and the difference is the loss.

2.1. Full Material Model

The time-dependent formulation of a random variable in generalised or polynomial KL-expansion is

$$\begin{array}{c}\hfill a(t,\mathbf{x})=a_0(t,\mathbf{x})+{\left(\sum _{m=1}^{\infty }{a}_m(t,\mathbf{x})y_m\right)}^k,(t,\mathbf{x})\in I\times D,\mathbf{y}={\left(y_m\right)}_{m=1}^{\infty }\in \Gamma ,k\ge 1\in \mathbb{N},\end{array}$$

(1)

where $a_m$ are the stochastic fluctuations, I is the time interval, D is the spatial computational domain, and $\Gamma$ is the space of i.i.d. random variables $y_m$. For the standard affine case, $k=1$. If one wanted to include multiple random variables, say, $a(t,\mathbf{x})$ and $b(t,\mathbf{x})$, then the random information $c(t,\mathbf{x})$ could simply be given as a product, for instance, $c(t,\mathbf{x})=a(t,\mathbf{x})\times b(t,\mathbf{x})$. For details, see [5,6].

Here, the focus is on uncertain domains, and therefore, only the shape of the capsule is taken to be random. Moreover, this effect is considered at a fixed point in time.

2.2. Filter

Even though time has been removed from the process, spatial localisation makes it possible to consider situations where different parts of the body are either exposed to their environment at different times or intensities. Butterworth filters are well suited to this task [16].

Suppose that $\Gamma$ is a positively oriented boundary of some domain $G\subset \mathbb{C}$. By the Cauchy integral formula

$$r\left(\xi \right)=\frac{1}{2\pi i}{\oint }_{\Gamma }{(z-\xi )}^{-1}dz=\left\{\begin{array}{c}1,\xi \in G,\hfill \\ 0,\xi \in G^c.\hfill \end{array}\right.$$

(2)

The idea is to approximate $r\left(\xi \right)$ with $r_N\left(\xi \right)={\sum }_{k=0}^{N-1}{w}_k{(z_k-\xi )}^{-1}$, for some $w_k,z_k\in \mathbb{C}$.

Definition 1

(Butterworth filter). Selecting the points and weights for $k=0,1,\dots ,N-1,$

$$z_k=\gamma e^{i({\theta }_k+\varphi )}+y,w_k=\gamma e^{i({\theta }_k+\varphi )}/N,$$

(3)

with ${\theta }_k=2\pi k/N$ and $\varphi =\pm \pi /N$. These weights and points are obtained using the N-point uniform trapezoid rule approximation of the contour integral in (2) when Γ is set to the circle $\Gamma =\{\gamma e^{i(\theta +\varphi )}+y:\theta \in [0,2\pi )\}$ enclosing the interval of interest.

In Figure 2, the effect of filtering on the degradation process is illustrated. Two filters are applied with variation in N and $\gamma$.

3. Model Problem: Shell of Revolution

Every shell of revolution is defined by its profile function, denoted here simply as $f\left(x\right)$, where x is the axial coordinate, $x\in [x_0,x_1]$ (see Figure 3). Full 3D analysis, especially with uncertain domain, is very expensive. Since it is assumed that the shells are thick, in other words, the thin structure assumptions are not valid, the reduced model has to be chosen with care. Under the assumption that the loading can be expanded as a Fourier series in the angular direction, it is possible to transfer the problem into cylindrical coordinates, where the computational domain is a 2D section that revolves around the axis of revolution. The elasticity formulation used in simulations is outlined in Appendix A.

The computational domain $\widehat{\Omega }$ is defined as the region bounded by the uncertain top profile function $f^{\omega }\left(x\right)$ and the deterministic bottom profile function $f_b\left(x\right)$. Formally, one has (see Figure 4)

$$\widehat{\Omega }=\left\{(x,r)\right|a\le x\le b,f_b\left(x\right)<r<f^{\omega }\left(x\right)\}.$$

(4)

In the following, the $f_b\left(x\right)$ is always constant, representing the inner radius of the capsule.

In deterministic problems, the thickness of the shell, d, is assumed to be constant. If the ratio of the thickness and the radius of the shell, R, is less than 1/10, it is often advantageous to use dimensionally reduced models where the thickness is simply a (dimensionless) parameter $t=d/R$. This division between thin and thick shells is not arbitrary, but based on the asymptotic accuracy of the reduced models. For instance, in Figure 4, the local dimensionless thickness is within 1/10 and 1/5.

Boundary Layers

Another option would be to use so-called hierarchic models to partly reduce the computational cost of the deterministic part [17]. This has not been considered here, since the connection of the local curvature changes and the induced boundary layers is not clear. The full 3D model (2D section in the Fourier setting) will for certain be able to capture these effects, provided the underlying deterministic finite element method is sufficiently accurate.

In every shell, there is a layer of characteristic length scale equal to the thickness of the shell (first layer). As the thickness approaches the limit for thick shells, the second layer∼square root of the thickness is likely to emerge.

4. Stochastic Finite Element Methods

Standard stochastic finite element methods can be divided into intrusive and nonintrusive ones. Nonintrusive ones have the benefit of preserving the existing software infrastructure and investment. Intrusive approaches at the moment appear to have some benefits if adaptivity is available. To date, adaptivity in domain uncertainty problems has not been studied. There are two important issues affecting the choice of approach: the cost of integration of the parameter-dependent part, and the properties of the parametrisation, that is, how does the domain change as a function of the random input.

The stochastic collocation method considered here is based on the assumption that the dependence on the abstract variable $\omega$ can be parameterised using a countable set of mutually independent random variables. Such a parametrisation can sometimes be known a priori or it can be estimated statistically.

4.1. Deterministic Part: p-Version

In order to accurately model and simulate the effects of changes in the domain, a method capable of accurate geometry description is needed. The p- and $hp$-finite element methods are ideally suited to this requirement. The paper by Babuška and Suri [18] gives an accessible overview of the method. For a more detailed exposition, see Schwab [19], and for those familiar with engineering approach, the book by Szabo and Babuška [13] is recommended.

The selection of shape functions is not unique. The Szabo and Babuška-style hierarchic shape functions are based on the integrated Legendre polynomials: For $x\in [-1,1]$

$${\varphi }_n\left(\xi \right)=\sqrt{\frac{2n-1}{2}}{\int }_{-1}^{\xi }{P}_{n-1}\left(t\right)dt,n=2,3,\dots$$

(5)

The normalising coefficients are chosen so that ${\int }_{-1}^1\frac{d{\varphi }_i\left(\xi \right)}{d\xi }\frac{d{\varphi }_j\left(\xi \right)}{d\xi }d\xi ={\delta }_{ij},i,j\ge 2.$ For a quadrilateral reference element over the domain $[-1,1]\times [-1,1]$, the shape functions are divided into three categories: nodal shape functions, side modes, and internal modes. There are four nodal shape functions, $4(p-1)$ side modes, and $(p-1)(p-1)$ interior modes per quadrilateral with order $(p\ge 2)$. Taken alone, the nodal shapes define the standard four-node quadrilateral finite element. For the internal modes have the tensor product structure

$$N_{i,j}^0(\xi ,\eta )={\varphi }_i\left(\xi \right){\varphi }_j\left(\eta \right),i=2,\dots ,p,j=2,\dots ,p.$$

(6)

The internal shape functions are often referred to as bubble-functions.

The Legendre polynomials have the property $P_n(-x)={(-1)}^n{P}_n\left(x\right)$, which is inherited by the integrated Legendre polynomials. It follows that since in 2D all internal edges of the mesh are shared by two different elements, additional bookkeeping is required to ensure that each edge has the same global parameterisation in both elements.

The option of increasing discretisation accuracy with the polynomial order makes it possible to have as large elements as possible, provided the possibly curved boundary segments are represented accurately. The linear blending function method of Gordon and Hall [20] is the standard choice for this purpose.

In the general case, all sides of an element can be curved; however, in the examples considered in this study, only one edge is curved, as in Figure 4. If every side is parameterised:

$$x=x_i\left(t\right),y=y_i\left(t\right),-1\le t\le 1,i=1,\dots ,\mathrm{number} \mathrm{of} \mathrm{sides},$$

(7)

then using capital letters as coordinates of the corner points, $(X_i,Y_i)$, the mapping for the global x-coordinates of a quadrilateral can be written as

$$\begin{array}{cc}\hfill x& =\frac{1}{2}(1-\eta )x_1\left(\xi \right)+\frac{1}{2}(1+\xi )x_2\left(\eta \right)+\frac{1}{2}(1+\eta )x_3\left(\xi \right)+\frac{1}{2}(1-\xi )x_4\left(\eta \right)\hfill \\ \hfill & -\frac{1}{4}(1-\xi )(1-\eta )X_1-\frac{1}{4}(1+\xi )(1+\eta )X_2-\frac{1}{4}(1+\xi )(1+\eta )X_3\hfill \\ \hfill & -\frac{1}{4}(1-\xi )(1+\eta )X_4,\hfill \end{array}$$

(8)

and symmetrically for the y-coordinate. If the side parametrisations represent straight edges, the mapping simplifies to the standard bilinear mapping of quadrilaterals.

Consider the abstract problem setting with u defined on the standard piecewise polynomial finite element space on some discretisation $\mathcal{T}$ of the computational domain $\Omega$. Assuming that the exact solution $u\in H_0^1\left(D\right)$ has finite energy, we arrive at the approximation problem: Find $\widehat{u}\in V$ such that

$$a(\widehat{u},v)=l\left(v\right)(=a(u,v))(\forall v\in V),$$

(9)

where $a(·,·)$ and $l\left(·\right)$, are the bilinear form and the load potential, respectively. Additional degrees of freedom can be introduced by enriching the space V. This is accomplished via introduction of an auxiliary subspace or “error space” $W\subset H_0^1\left(D\right)$ such that $V\cap W=\left\{0\right\}$. We can then define the error problem: Find $\epsilon \in W$ such that

$$a(\epsilon ,v)=l\left(v\right)-a(\widehat{u},v)(=a(u-\widehat{u},v))(\forall v\in W).$$

(10)

This can be interpreted as a projection of the residual to the auxiliary space. Here, the solution $ϵ$ is the error function, which is defined on the same discretisation as the solution $\widehat{u}$. It follows that any postprocessing that can be performed for the solution is also available for the error function.

4.2. Stochastic Collocation

For the uncertain profile function $f^{\omega }$, a representation of the form

$$\begin{array}{c}\hfill f^{\omega }\left(\mathbf{x}\right)=f_0\left(\mathbf{x}\right)-{\left(\sum _{m=1}^{\infty }{f}_m\left(\mathbf{x}\right)Y_m\left(\omega \right)\right)}^2,\end{array}$$

(11)

is assumed. Here, ${\left\{f_m\right\}}_{m\ge 0}$ are spatial stochastic fluctuations decaying in $L^{\infty }$-norm as $m\to \infty$, and ${\left\{Y_m\right\}}_{m\ge 1}$ is a family of random variables assumed to be i.i.d. The series (11) can be viewed as a special form of a Karhunen–Loevé expansion of the random underlying random field, in which case the functions $f_m\left(\mathbf{x}\right)$ are the eigenfunctions of some covariance function $C(x_1,x_2)$. For a Galerkin SFEM perspective and further references, see [21]. Squaring of the sum of the random input ensures that the perturbation does not change its sign. In some simple geometric settings, these eigenfunctions are known explicitly, but they can also be resolved numerically as eigenvectors of the sampled covariance matrix; see, e.g., ([22], Chapter 8), [23,24]. For numerical computations, the series is truncated after M terms. For the model to be applicable as $m\to \infty$, the decay of of the stochastic fluctuations must be greater than quadratic, which is the critical decay rate (see [14]). Notice that in all truncated formulations, the problem is well posed, even with decay rate $=2$.

In the following, an anisotropic Smolyak-type collocation operator defined with respect to a finite multi-index set (see, e.g., [25,26,27,28]) is considered. In principle, a simpler formulation would also be sufficient for the purpose of our numerical experiments. However, the general form of the operator is kept here, since it allows for efficient computation even in high-dimensional parameter spaces.

For the sake of accuracy, it is standard practice to choose the collocation points to be the abscissae of orthogonal polynomials; see [29]. Here, the collocation method is formulated using Legendre polynomials, which are the optimal choice when the input random variables are uniform. If a lognormal random field model had been adopted instead of the squared one, one should use Hermite polynomials instead, but otherwise, the collocation method remains the same.

Let $L_p$ be the univariate Legendre polynomial of degree p. Denote by ${\left\{{\chi }_k^{\left(p\right)}\right\}}_{k=0}^p$ the zeros of $L_{p+1}$ and by ${\left\{w_k^{\left(p\right)}\right\}}_{k=0}^p$ the associated Gauss–Legendre quadrature weights. The one-dimensional Lagrange interpolation operators ${\mathcal{I}}_p^{\left(m\right)}$ are defined via $\left({\mathcal{I}}_p^{\left(m\right)}v\right)\left(Y_m\right)={\sum }_{k=0}^{p}v\left({\chi }_k^{\left(p\right)}\right){\ell }_k^{\left(p\right)}\left(Y_m\right),$ where ${\left\{{\ell }_k^{\left(p\right)}\right\}}_{k=0}^p$ are the related Lagrange basis polynomials of degree p. Observe that ${\int }_{-1}^1\left({\mathcal{I}}_p^{\left(m\right)}v\right)\left(Y_m\right)\frac{d{Y}_m}{2}={\sum }_{k=0}^{p}v\left({\chi }_k^{\left(p\right)}\right)w_k^{\left(p\right)}.$

Now, let $\mathcal{A}\subset {\mathbb{N}}_0^M$ be a finite set of multi-indices. For $\alpha ,\beta \in \mathcal{A}$, one can write $\alpha \le \beta$ if ${\alpha }_m\le {\beta }_m$ for all $m=1,\dots ,M$. It is implicitly assumed that $\mathcal{A}$ is monotone in the following sense: $\exists \alpha \in \mathcal{A}\mathrm{such}\mathrm{that}\beta \le \alpha \Rightarrow \beta \in \mathcal{A}.$ The sparse collocation operator is defined as ${\mathcal{I}}_{\mathcal{A}}:={\sum }_{\alpha \in \mathcal{A}}{⨂}_{m=1}^M\left({\mathcal{I}}_{{\alpha }_m}^{\left(m\right)}-{\mathcal{I}}_{{\alpha }_m-1}^{\left(m\right)}\right),$ where ${\mathcal{I}}_{-1}^{\left(m\right)}:=0.$ This may be rewritten in a computationally more convenient form:

$${\mathcal{I}}_{\mathcal{A}}=\sum _{\alpha \in \mathcal{A}}\sum _{\beta \in {\mathcal{G}}_{\alpha }}{(-1)}^{{\left|\right|\alpha -\beta \left|\right|}_1}\underset{m=1}{\overset{M}{⨂}}{\mathcal{I}}_{{\beta }_m}^{\left(m\right)},$$

(12)

where ${\mathcal{G}}_{\alpha }:=\{\beta \in {\mathbb{N}}_0^M|\alpha -1\le \beta \le \alpha \}$. The collocation points are of the form

$$\begin{array}{c}\hfill {\chi }_{\gamma }^{\left(\alpha \right)}=\left({\chi }_{{\gamma }_1}^{\left({\alpha }_1\right)},\dots ,{\chi }_{{\gamma }_M}^{\left({\alpha }_M\right)}\right)\in \Gamma \end{array}$$

(13)

for some $\gamma \in {\mathbb{N}}_0^M$ such that $\gamma \le \alpha \in \mathcal{A}$. Similarly, the tensorised quadrature weights are defined as $w_{\gamma }^{\left(\alpha \right)}=w_{{\gamma }_1}^{\left({\alpha }_1\right)}\cdots w_{{\gamma }_M}^{\left({\alpha }_M\right)}$. Higher moments, such as the expected value and variance, for the collocated solution, may now be computed by applying the quadrature rule

$$\mathcal{E}\left[\underset{m=1}{\overset{M}{⨂}}{\mathcal{I}}_{{\alpha }_m}^{\left(m\right)}v\right]=\sum _{\gamma \le \alpha }v\left({\chi }_{\gamma }^{\left(\alpha \right)}\right)w_{\gamma }^{\left(\alpha \right)}$$

(14)

to the terms in (12).

The accuracy of the collocated approximation is ultimately determined by the smoothness of the solution as well as the choice of the multi-index set $\mathcal{A}\subset {\mathbb{N}}_0^M$; see [26,28] for a detailed analysis. In this paper, the tensor product grids defined by

$$\mathcal{A}=\{\alpha \in {\mathbb{N}}_0^M|{\alpha }_m\le p,m=1,\dots ,M\},p\in {\mathbb{N}}_0$$

(15)

are used in the simulations below.

4.3. Validation: Monte Carlo

Monte Carlo methods are robust in the statistical sense, yet often inefficient in terms of computational complexity. The proposed collocation method is validated against the Monte Carlo (see Figure 5). The collocation scheme is used to compute the expected value of the $L^2$-norm of the transverse deflection. This is taken to be the reference value, and the observed error within the Monte Carlo is taken to be the $L^2$-error of the transverse deflection. As can be seen, the Monte Carlo error decreases, as predicted by the theory.

5. Numerical Simulations

In the set of numerical simulations, the stochastic fluctuations are set to be $f_m\left(\xi \right)=sin(log\left(100\right)m\xi )/{(m+1)}^{\sigma }$, and the KL-expansions are truncated at $M=5$. The external loading is assumed to constant, $l\left(x\right)=1$, and acting on the w-component. The simulations are divided into four groups with initial working hypothesis:

• Effect of decay of parameters—as $\sigma$ increases, the observed amplitudes should decrease.
• Effect of local thickness—As $f_b\left(\mathbf{x}\right)\to f^{\omega }\left(\mathbf{x}\right)$, the local changes in curvature start dominating. In other words, aggregate quantities of interest such as norms of displacement fields should start increasing.
• Localised variation—If the variation is removed from the fixed boundaries, the layer effects are significantly reduced.
• Higher Fourier modes—This should primarily affect the deterministic part if the stochastic setup is kept constant.

The fifth and final section after the simulations concerns the statistics of the error. Since the error estimate has the same representation as the actual finite element solution, exactly the same procedures can be applied to it. Indeed, it is possible to consider the reliability of the simulations by considering the expectation of the error and even its variance. This is of interest particularly when there is an a priori understanding of the underlying results that is contradicted by the simulations.

The nominal domain is $\Omega =[-1,1]\times [f_b,1]$ in all cases. Loading is a transverse unit traction load acting on the top surface with the angular wave number $K=2$. The Poisson number $\nu =1/3$. Since the Young’s modulus is assumed to be constant, it has been set to a unit value $E=1$, and hence, all displacements are given in scaled units that should be interpreted only in relative terms. Simulation data are tabulated in

Table 1. Simulation setup.

FEM	Simulations: $p=4$, auxiliary space: $p=6$ (interior), $p=5$ (edges).
Mesh	Nodes: 63, edges: 102, interiors: 40.
DOF	Approximation: 2187, auxiliary space: 2226.
Loading	Transverse unit traction load acting on the top surface ($K=2$).
Stochastic Dimension	$M=5$, full tensor basis, 243 collocation points.
Filters	(boundary) $N=4$, $\gamma =1/4$, $\varphi =\pi /4$, and $y=\pm 1$,
	(interior) $N=10$, $\gamma =2/3$, $\varphi =\pi /4$, and $y=0$.

, and the specific parametrisations of various simulation sets are given in

Table 2. Outline of the simulation parametrisation.

Simulation	Parameter	Range
Decay	$\sigma$	$\sigma \in [2,4]$
Thickness	d	(upper bound) $d\le [0.1,0.2]$
Location	$x_0,R$	$x_0=0,R=2/3$
Fourier Mode	K	$K=4$

.

Finally, in all simulation sets, the parameter runs from low (blue) to high (purple).

5.1. Effect of Decay of Parameters

The simulation results are summarised in Figure 6, Figure 7 and Figure 8. In all cases, the perturbed domain $\Omega =[-1,1]\times [0.8,f^{\omega }\left(x\right)]$. As expected, as the decay rate $\sigma \to \infty$, the variance decreases. Since the overall variation is small, the expectations are not strongly affected by the decay rate.

The quantity of interest is the $L^2$-norm of the transverse displacement. The results are computed using the true perturbed domain and the pulled-back field on the nominal domain. The smaller perturbations are reflected on the norms as well. The convergence graph of the modelling error (see Figure 7b) shows that the difference between the true and the mapped fields is not significant. This suggests that this admittedly crude approach can be used as the first approximation.

Of particular interest is the deformation profile detail of Figure 7d. The shortest boundary layer proportional to thickness is clearly visible. This shortest length layer is often omitted from the considerations in the reduced models, but in this “thick” shell regime, it does play a role. Notice that the next axial layer proportional to square root of thickness extends over the whole domain, if both ends are considered.

Variance results are in line with the a priori knowledge. Variance results are computed on the nominal domain only, and hence, numerical instability at the boundaries is observed. Interestingly, in Figure 8c, two things are evident: first, the increase in variance is logarithmic; second, the logarithmic scaling also applies to the distance from the boundary. Geometrically, this means that with smaller values of $\sigma$, the variance in the boundary layer profiles is higher.

5.2. Effect of Local Thickness

In this set of simulations, everything else is kept constant except the bottom radius. As the bottom radius increases, the local effective thickness decreases. Since the profile function is the same ($\sigma =2.5$ constant), only thickness variation effects are observed.

As the local thickness tends to zero, the longer-term boundary layers start to play a more important role. This cannot be detected in the expectation, but in Figure 9d, the variance corresponding to the thinnest case displays more curvature (notice the loglog-scale!). For a thinner shell, the second layer is shorter, but the moment it induces causes relatively stronger displacements closer to the boundary. This is an excellent example of usefulness of uncertainty quantification.

5.3. Localised Variation

Following with the setup of the previous set of simulations, now the filtering is changed to allow for variation only away from the boundaries. Comparison of Figure 9d and Figure 10d reveals the effect of the shortest layer. Since in Figure 10d there is very little variance at the boundary, only the second layer effect is visible. In Figure 9d, the stronger overall effect is due to the linear combination of the first two layers.

5.4. Higher Fourier Mode Loading

One further parameter to consider is the wave number K of the loading. Again, the simulation setup is the one used in the previous set. Since the change in loading effectively affects only the deterministic part, it is not surprising that at $K=4$, the result changes significantly only in the $L^2$-norm of the expectation of the transverse displacement. The maximal amplitudes are smaller, since the higher angular oscillations induces membrane effects, that is, the structure is slightly more resistant to bending (see Figure 11).

5.5. Error Estimation

Finally, the auxiliary space error estimation provides a unique opportunity to study the reliability of the simulation via the statistics of the error function itself. Here, the setup is the same as with the decay rate simulations. Notice that the number of degrees of freedom for error estimation is greater than that of the actual solutions. This is not unusual in the p-version with a small number of elements. The surface plots of Figure 12 indicate that the errors concentrate on the boundaries as one would expect, and indeed, the crude mapping of the domain leads to high variance of the error along the top surface.

The decay rate affects the size of the error in a predictable way. In Figure 13 the specific feature of the auxiliary space error estimation is visible; namely, since the model degrees of freedom are not included in the approximation space, the error representation always has an oscillating nature.

Also, the logarithmic dependence of the variance to the decay rate is naturally also present in the error function. Here, the numerical instabilities prevent full resolution of the behaviour close to the boundary (Figure 14).

6. Conclusions

Structures made of biodegradable materials can be analysed using current simulation techniques, yet with significantly increased complexity. In this study, it was naively assumed that the materials can be modelled as isotropic–possibly after homogenisation—and that reliable data exist for derivation of the necessary degradation model. Both of these assumptions are difficult to realise in practice.

For capsules, the degradation leads to domain uncertainty. It is clear that even in the so-called thick shell class of structures, the boundary layers do play a role, and their presence is seen in higher statistical moments such as variance. In such multiparametric problems, having reliable error estimation techniques available is of great importance. Here, the novelty lies in the application of the error function, which is represented on the same discretisation as the individual finite element solutions. This leads to a possibility to study the error using the same tools as the actual solution. In particular, in the examples above it have been identified not only where the error is likely to be largest, but also where the variance of the error is largest.

Many of the proposed designs for long-term drug delivery capsules have geometric properties that simplify their analysis. It is clear that full 3D analysis under realistic degradation models remains very expensive in all phases of simulation spanning, from modelling to numerical simulations.