`Law of Laplace`

Module: tutorials.04_EM_tissue.02_laplace.run

Section author: Matthias Gsell <matthias.gsell@medunigraz.at>

This example verifies the computed stress and the computed work using a simple spherical shell geometry and time-varying Neumann boundary condition applied on the inner boundary. To verify the computed stress, we use the law of Laplace, see [2], to estimate the wall stress and compare it against the mean computed stress. For the work varification, we compute the external work analytically and compare it against the computed internal work.

Stress Verification

To verify the stress, we use the law of Laplace to estimate the wall stress (wall tension).The law of Laplace can be used to describe the wall stress of a spherical shell in terms of the surface pressure applied on the inner boundary, see Figure Fig. 132. To estimate the wall tension, we consider the geometry at time $t = T$ with pressure $p(t=T) = p$ and radii $r(t=T) = r, R(t=T) = R$ .

Fig. 132 Deformation at time step $t=0$ and $t=T$ .

Derivation Of Laplace’s Law

To derive Laplace’s law we want to use the balance of forces. We assume that the center of the spherical shell is the origin and we denote the inner radius of the deformed geometry by $r$ and the outer radius by $R$ with $R > r$ , see figure Fig. 133.

Fig. 133 Geometrical setting.

By $h$ we denote the wall thickness which is given by $h = R - r$ . Next, we clip the geometry using the x-y-plane as a clipping plane and obtain two equally sized shells, see figure Fig. 134.

Fig. 134 Acting forces.

The force vector acting on the inner surface is given in spherical coordinates by

$\mathbf{F}(\theta, \varphi) = p \left( \begin{array}{c} \sin(\theta) \, \cos(\varphi) \\ \sin(\theta) \, \sin(\varphi) \\ \cos(\theta) \end{array} \right)$

with $\theta \in [0, \pi]$ and $\varphi \in [0, 2 \pi]$ . The components of the total force vector $\mathbf{F}^{tot}$ are

(41) $F_x^{tot} & = \int \limits_{\theta=0}^{\frac{\pi}{2}} \int \limits_{\varphi=0}^{2 \pi} \big( \mathbf{F} \cdot \mathbf{e}_x \big) \, r^2 \, \sin(\theta) \mathrm{d} \varphi \, \mathrm{d} \theta = p \, r^2 \int \limits_{\theta=0}^{\frac{\pi}{2}} \int \limits_{\varphi=0}^{2 \pi} \sin^2(\theta) \, \cos(\varphi) \mathrm{d} \varphi \, \mathrm{d} \theta = 0, \\ F_y^{tot} & = \int \limits_{\theta=0}^{\frac{\pi}{2}} \int \limits_{\varphi=0}^{2 \pi} \big( \mathbf{F} \cdot \mathbf{e}_y \big) \, r^2 \, \sin(\theta) \mathrm{d} \varphi \, \mathrm{d} \theta = p \, r^2 \int \limits_{\theta=0}^{\frac{\pi}{2}} \int \limits_{\varphi=0}^{2 \pi} \sin^2(\theta) \, \sin(\varphi) \mathrm{d} \varphi \, \mathrm{d} \theta = 0, \\ F_z^{tot} & = \int \limits_{\theta=0}^{\frac{\pi}{2}} \int \limits_{\varphi=0}^{2 \pi} \big( \mathbf{F} \cdot \mathbf{e}_z \big) \, r^2 \, \sin(\theta) \mathrm{d} \varphi \, \mathrm{d} \theta = p \, r^2 \int \limits_{\theta=0}^{\frac{\pi}{2}} \int \limits_{\varphi=0}^{2 \pi} \sin(\theta) \, \cos(\theta) \mathrm{d} \varphi \, \mathrm{d} \theta = p \, r^2 \, \pi \\$

with $\theta \in (0, \frac{\pi}{2})$ since we only integrate over one half of the shell.

For the derivation of Laplace’s law we approximate the wall stress $\sigma$ by its mean value $\bar{\sigma}$ , see figure Fig. 134. Since the geometry and the acting forces are symmetric, the tangential stress in any direction must be the same and there will be zero shear stress. Due to the assumption, the mean wall stress at the cut face just acts just in z-direction, thus we have $\bar{\sigma} \cdot \mathbf{e}_x = 0$ and $\bar{\sigma} \cdot \mathbf{e}_y = 0$ . The components of the total wall stress vector $\bar{\sigma}^{tot}$ are

(42) $\bar{\sigma}_x^{tot} & = \int \limits_{r=r}^{R} \int \limits_{\varphi=0}^{2 \pi} \big( \bar{\sigma} \cdot \mathbf{e}_x \big) \, r^2 \, \sin(\theta) \mathrm{d} \varphi \, \mathrm{d} r = 0, \\ \bar{\sigma}_y^{tot} & = \int \limits_{r=r}^{R} \int \limits_{\varphi=0}^{2 \pi} \big( \bar{\sigma} \cdot \mathbf{e}_y \big) \, r^2 \, \sin(\theta) \mathrm{d} \varphi \, \mathrm{d} r = 0, \\ \bar{\sigma}_z^{tot} & = \int \limits_{r=r}^{R} \int \limits_{\varphi=0}^{2 \pi} \big( \bar{\sigma} \cdot \mathbf{e}_z \big) \, r^2 \, \sin(\theta) \mathrm{d} \varphi \, \mathrm{d} r = \pi \, \big( R^2 - r^2 \big) \, \sigma_W$

with $\sigma_W := \bar{\sigma} \cdot \mathbf{e}_z$ . From equations (41) and (42) we have

(43) $\mathbf{F}^{tot} = \bar{\sigma}^{tot} \quad \Leftrightarrow \quad p \, r^2 = (R^2 - r^2) \, \sigma_W.$

From (43) we conclude that

(44) $\sigma_W = \frac{p \, r^2}{R^2 - r^2} = \frac{p \, r}{2 \, h \left( 1 + \frac{h}{ 2 \, r } \right) }$

and for $h \ll r$ we obtain the law of Laplace

(45) $\sigma_W \approx \sigma_{L} = \frac{p \, r}{2 \, h}$

which can be used to estimate the wall stress for thin walled spherical shells.

Modeling Error

The law of Laplace, see equation (45), is just an approximation to the more accurate formula (44). The approximation is justified if the wall thickness is mush smaller than the inner radius. To figure out how strong the relative modeling error depends on $R$ and $r$ , we rewrite the error as

(46) $err_{\text{rel}} := \left| \frac{\sigma_W - \sigma_L}{\sigma_W} \right| = \frac{R - r}{2 \, r} = \frac{h}{2 \, r}$

which yields to the condition

$\frac{h}{r} \stackrel{!}{<} 2 \, \varepsilon$

if we want to ensure that the relative modeling error is less than $\varepsilon \in \mathbb{R}_+$ .

Example

Assume, that we have a spherical shell with $r = 30 mm$ and we want to use the law of Laplace, equation (45), to estimate the surface tension. To ensure that the error is less than one present, that is $\varepsilon = 0.01$ , we have to ensure that the wall thickness satisfies $h < 0.6 mm$ . For $\varepsilon = 0.1$ we have $h < 6 mm$ and for $\varepsilon = 0.25$ we have $h < 15 mm$ .

Verification

To verify the stress, we first compute the stress tensor $\sigma_C$ elementwise and project the tensor from the cartesian coordinate system to the spherical coordinate system, i.e.

(47) $\sigma_S = \left( \begin{array}{ccc} \sigma_{rr} & \sigma_{r \varphi} & \sigma_{r \theta } \\ \sigma_{\varphi r} & \sigma_{\varphi \varphi} & \sigma_{\varphi \theta } \\ \sigma_{\theta r} & \sigma_{\theta \varphi} & \sigma_{\theta \theta} \end{array} \right) = \mathbf{P}^T \, \underbrace{ \left( \begin{array}{ccc} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{zx} & \sigma_{zy} & \sigma_{zz} \end{array} \right)}_{= \mathbb{\sigma}_C} \mathbf{P}$

with the projection matrix $\mathbf{P} = \big( \mathbf{e}_r \; \mathbf{e}_{\varphi} \; \mathbf{e}_{\theta} \big)$ . The mean stress tensor then is

$\bar{\sigma}_S := \frac{1}{| \Omega(T) |} \int \limits_{\Omega(T)} \sigma_S \, \mathrm{d} \, \mathbf{x}$

where the integration is applied componentwise.

Experiments

For the experiments we choose the following setups with two different geometries. The first geometry is a thin walled spherical shell and the second geometry is a rather thick walled spherical shell. The material model is the Demiray model with $\kappa = 1000.0$ , $a=34.0$ and $b=10.0$ .

	$r(t=0)$ [mm]	$R(t=0)$ [mm]	$h(t=0)$ [mm]	$p = p(t=T)$ [kPa]	$T$ [ms]
Experiment 1	15.0	15.5	0.5	2.0	100.0
Experiment 2	15.0	15.5	0.5	4.0	100.0
Experiment 3	15.0	30.0	15.0	2.0	100.0
Experiment 4	15.0	30.0	15.0	4.0	100.0

Using the method described in this section, we obtain the following results.

	$\sigma_W$ [kPa]	$\sigma_L$ [kPa]	$\bar{\sigma}_{rr}$ [kPa]	$\bar{\sigma}_{\theta \theta}$ [kPa]	$\bar{\sigma}_{\varphi \varphi}$ [kPa]
Experiment 1	39.2682	38.7745	-0.9303	38.8500	38.8516
Experiment 2	83.8585	84.8469	-1.8753	83.7231	83.7220
Experiment 3	0.6937	1.0304	-0.3601	0.6260	0.6271
Experiment 4	1.4428	2.1226	-0.7198	1.2951	1.2973

We see, that for the thin walled geometry the analytically computed wall stress $\sigma_W$ and the Laplace approximation $\sigma_L$ are almost equal to the computed mean tangential stresses $\bar{\sigma}_{\theta \theta}$ and $\bar{\sigma}_{\varphi \varphi}$ . The radial stress $\bar{\sigma}_{rr}$ is compared to the tangential stresses negligible. The relative modeling error for the first two experiments are $err_{\text{rel}} = 0.0061$ and $err_{\text{rel}} = 0.0056$ , see equation (46). For the experiments with a thick walled geometry we have a rather high discrepancy in the analytical and computed stresses. The relative modeling errors for the last two experiments are $err_{\text{rel}} = 0.4853$ and $err_{\text{rel}} = 0.4714$ .

Work Verification

To verify the computed work, we compute the external work analytically and compare it with the computed internal work. To compute the external work, we use the well known principle of pressure-volume work which is $\delta W = p(V) \, \mathrm{d} V$ where $V$ denotes the volume of the inner sphere ( $V(t) = \frac{4 \, \pi \, r(t)^3}{3}$ ). The integral representation then reads

(48) $W_{\text{ext}} = \int \limits_{V = V(t=0)}^{V(t=T)} p(V) \, \mathrm{d} V.$

Next, we assume that the pressure is a piecewise linear function of the radius, that is

(49) $p(V) = p(r) = \frac{p_{i}-p_{i-1}}{r_{i}-r_{i-1}} \, (r - r_{i-1}) + p_{i-1}$

for $r \in [r_{i-1}, r_{i}]$ and $i = 0, \ldots, n$ . We get the pressure and radius values $p_i, r_i$ from the simulation.

Using spherical coordinates and equation (49), we can write equation (48) as

(50) $W_{\text{ext}} & = \int \limits_{r = r(0)}^{r(T)} \, \int \limits_{\theta = 0}^{\pi} \, \int \limits_{\varphi = 0}^{2 \pi} r^2 \, \sin(\theta) \, p(r) \, \mathrm{d} \varphi \, \mathrm{d} \theta \, \mathrm{d} r = 4 \, \pi \, \sum \limits_{i=1}^{n} \int \limits_{r = r_{i-1}}^{r_i} r^2 \, p(r) \, \mathrm{d} r \\ & = 4 \, \pi \, \sum \limits_{i=1}^{n} \, \left( \frac{p_{i} - p_{i-1}}{12\, (r_{i} - r_{i-1})} \big( 3 \, r_{i}^4 - 4 \, r_{i}^3\, r_{i-1} + r_{i-1}^4 \big) + \frac{p_{i-1}}{3} \big( r_{i}^3 - r_{i-1}^3 \big) \right)$

which provides an analytical approximation of the external work in [kPa mm^3] = 10^(-6) [J].

To compute the internal work, we use the internal mechanical power $\mathcal{P}_{\text{int}}$ which is

(51) $\mathcal{P}_{\text{int}}(t) = \int \limits_{\Omega_{\text{ref}}} \mathbf{S} : \dot{\mathbf{E}} \, \mathrm{d} \mathbf{x},$

see [1] (Section 4.4) for further details. The internal mechanical work $W_{\text{int}}$ is then given by

(52) $W_{\text{int}} = \int \limits_{t=0}^{T} \mathcal{P}_{\text{int}}(t) \, \mathrm{d} t.$

Note

The formula in [3] for the internal heart power and its modification yield

(53) $\mathcal{P}_{\text{int},L}(t) = \sigma_L(t) \, V(t) \, 2 \, \dot{\bar{\varepsilon}}(t)$

where $V(t)$ denotes the volume of the geometry at time $t$ and $\dot{\bar{\varepsilon}}(t)$ denotes the mean strain rate with the mean strain $\bar{\varepsilon}(t) = \frac{1}{2} \, \left( \frac{r(t)}{r(0)} + \frac{R(t)}{R(0)} \right)$ . First, we assume that the stress, the volume and the radii depend piecewise linearly on time.

For the mean strain we obtain

(54) $\bar{\varepsilon}(t) \approx \frac{1}{2} \, \left( \frac{r_{i}-r_{i-1}}{r(0)} + \frac{R_{i}-R_{i-1}}{R(0)} \right) \, \frac{t-t_{i-1}}{t_{i}-t_{i-1}} + \frac{1}{2} \, \left( \frac{r_{i-1}}{r(0)} + \frac{R_{i-1}}{R(0)} \right) \qquad \text{ and } \qquad \dot{\bar{\varepsilon}}(t) \approx \frac{1}{2} \, \left( \frac{r_{i}-r_{i-1}}{r(0)} + \frac{R_{i}-R_{i-1}}{R(0)} \right) \, \frac{1}{t_{i}-t_{i-1}}$

with $t \in (t_{i-1}, t_i)$ and $r_i := r(t_i), R_i := R(t_i)$ . Thus, for $t \in (t_{i-1}, t_i)$ , we have

$\mathcal{P}_{\text{int},L}(t) \approx \left( \big( \sigma_{L,i}-\sigma_{L,i-1} \big) \, \frac{t-t_{i-1}}{t_i - t_{i-1}} + \sigma_{L,i-1} \right) \, \left( \big( V_i-V_{i-1} \big) \, \frac{t-t_{i-1}}{t_i - t_{i-1}} + V_{i-1} \right) \, \left( \frac{r_{i}-r_{i-1}}{r(0)} + \frac{R_{i}-R_{i-1}}{R(0)} \right) \, \frac{1}{t_{i}-t_{i-1}}$

with $\sigma_{L,i} := \sigma_L(t_i)$ and $V_i := V(t_i)$ .

By integrating over $t \in (t_{i-1}, t_i)$ we get

$\int \limits_{t=t_{i-1}}^{t_i} \mathcal{P}_{\text{int},L}(t) \, \mathrm{d} t \approx \frac{1}{3} \, \left( \frac{r_{i}-r_{i-1}}{r(0)} + \frac{R_{i}-R_{i-1}}{R(0)} \right) \, \left( \big( \sigma_{L,i} + \frac{1}{2} \, \sigma_{L,i-1} \big) \big( V_i + \frac{1}{2} \, V_{i-1} \big) + \frac{3}{4} \sigma_{L,{i-1}} V_{i-1} \right)$

and finally for the internal work $W_{\text{int},L}$ we have

$W_{\text{int},L} & = \int \limits_{t=0}^{T} \mathcal{P}_{\text{int},L}(t) \, \mathrm{d} t = \sum \limits_{i=1}^{n} \int \limits_{t=t_{i-1}}^{t_i} \mathcal{P}_{\text{int},L}(t) \, \mathrm{d} t \\ & \approx \sum \limits_{i=1}^{n} \frac{1}{3} \, \left( \frac{r_{i}-r_{i-1}}{r(0)} + \frac{R_{i}-R_{i-1}}{R(0)} \right) \, \left( \big( \sigma_{L,i} + \frac{1}{2} \, \sigma_{L,i-1} \big) \big( V_i + \frac{1}{2} \, V_{i-1} \big) + \frac{3}{4} \sigma_{L,{i-1}} V_{i-1} \right).$

To obtain $W_{\text{int},W}$ we proceed in the same way by replacing $\sigma_L$ with $\sigma_W$ .

Experiments

Using the same geometrical settings as in Experiments we obtain the following results.

	$W_{\text{ext}}$ [J]	$W_{\text{int}}$ [J]	$W_{\text{int},L}$ [J]	$W_{\text{int},W}$ [J]
Experiment 1	0.00407	0.00420	0.00434	0.00428
Experiment 2	0.01072	0.01112	0.01181	0.01167
Experiment 3	0.00074	0.00070	0.00097	0.00065
Experiment 4	0.00303	0.00285	0.00396	0.00268

Meshes

Download the following Python-script and run it!

References

[1]	Holzapfel, G.A., Nonlinear Solid Mechanics: A Continuum Approach for Engineering (2000), West Sussex, England: John Wiley & Sons, Ltd

[2]	Grossman, W. and Jones, D. and McLaurin, L.P., Wall stress and patterns of hypertrophy in the human left ventricle (1975), American Society for Clinical Investigation

[3]

Fernandes, J.F. and Goubergrits, L. and Brüning, J. and Hellmeier, F. and Nordmeyer, S. and da Silva, T.F. and Schubert, St. and Berger, F. and Kuehne, T. and Kelm, M. and others, Beyond Pressure Gradients: The Effects of Intervention on Heart Power in Aortic Coarctation (2017), Public Library of Science