Printing non-Euclidean solids

TL;DR

This paper explores how incompatible surface growth can be used to control the behavior of non-Euclidean solids, with applications in 3D printing arteries and explosive plants, by tuning geometrical frustration during deposition.

Contribution

It introduces a mechanics framework for incompatible surface growth to achieve targeted configurations in non-Euclidean solids, including explicit protocols for biological and technological applications.

Findings

Stress uniformity in arteries under inhomogeneous loads

Complete energy release in explosive plants with a single cut

Topological incompatibility is essential for physiological targets

Abstract

Geometrically frustrated solids with non-Euclidean reference metric are ubiquitous in biology and are becoming increasingly relevant in technological applications. Often they acquire a targeted con- figuration of incompatibility through surface accretion of mass as in tree growth or dam construction. We use the mechanics of incompatible surface growth to show that geometrical frustration develop- ing during deposition can be fine-tuned to ensure a particular behavior of the system in physiological (or working) conditions. As an illustration, we obtain an explicit 3D printing protocol for arteries, which guarantees stress uniformity under inhomogeneous loading, and for explosive plants, allowing a complete release of residual elastic energy with a single cut. Interestingly, in both cases reaching the physiological target requires the incompatibility to have a topological (global) component.

Equations

$$\text{div}{\mathbf{\sigma }}_p=0\text{in}{\mathcal{B}}_p,{\mathbf{\sigma }}_p\text{n}={\mathbf{s}}_p\text{on}{\omega }_p$$

$$\left\{\begin{array}[]{llll}\text{div}\,\dot{\bm{\sigma}}=\textbf{0}&&\text{in}&\mathscr{B}(t)\\ \dot{\bm{\sigma}}\textbf{\emph{n}}=|\nabla\tau|^{-1}\text{div}(\mathring{\bm{\sigma}}_{a}+\textbf{\emph{p}})&&\text{on}&\omega(t).\end{array}\right.$$

$$\mathring{\mathbf{\eta }}-\nabla \tau \times {\left[\text{curl}\dot{\mathbf{ϵ}}\right]}_{\omega }^{\text{T}}-\text{curl}\left[\nabla \tau \times {\dot{\mathbf{ϵ}}}_{\omega }\right]={\mathbf{\eta }}_p\text{in}\mathcal{B}(t_f).$$

$$\dot{u}(r,R)=\frac{R^{2}f(R)}{\mu r_i^2+(\mu +\lambda )R^2}\frac{r^2-r_i^2}{2r}.$$

$${\mathring{\sigma }}_{a_{\theta }}(R)=\frac{4\mu (\mu +\lambda ){\int }_{r_i}^R\left(\mu r_i^2+(\mu +\lambda )r^2\right){\phi }_p(r)dr}{(2\mu +\lambda )(\mu r_i^2+(\mu +\lambda )R^2)}.$$

$$\frac{{\mathring{\sigma }}_{a_{\theta }}(R)}{p}=\frac{r_i{r}_f(R-r_i)(\mu r_i+(\mu +\lambda )R)}{(r_f-r_i)R(\mu r_i^2+(\mu +\lambda )R^2)}$$

$$\frac{{\mathring{\sigma }}_{a_{\theta }}(R)}{c}=-\frac{2\mu (\lambda +\mu )}{\lambda +2\mu }\frac{(\mu +\lambda )(r_i^2-R^2)+2r_i^2\mu \log (r_i/R)}{\mu r_i^2+(\mu +\lambda )R^2}.$$

$$\pi (R)=\frac{2(\mu +\lambda )}{r_i^2(2\mu +\lambda )}{\int }_{r_i}^R\left(\mu r_i^2+(\mu +\lambda )r^2\right){\phi }_p(r)dr.$$

Full text

Externally unloaded elastic solids can be still endogenously pre-stressed by distributed self-equilibrated force couples. In living organisms such pre-conditioning is a way of achieving specific targets in physiological regimes Gladman , for instance, vegetable leafs require residual stresses to open Oliver while arteries need a pre-load to ensure transmural uniformity of hoop stress ChuFung . Residual stresses are equally important in engineering applications, where pre-loading is used either for reinforcement or to delay the onset of failure ResStr2 ; ResStrConcr . Furthermore, challenging new applications, like the design of programmable bio-mimetic materials, depend crucially on our ability to create complex patterns of residual stresses Ge ; Kempaiah ; LindNature16 ; Geitmann ; Danescu .

In this Letter we address the question how a particular distribution of residual stresses can be produced in a solid as a result of surface accretion of mass, a central process in both natural growth and 3D printing.

The source of residual stresses in elastic solids is the incompatibility of the reference configuration preventing its isometric embeding into the Euclidean space. A reference (natural) state is characterised by a metric tensor and in Euclidean solids this metric is flat Sharon2007 ; Efrati2009 . In “non-Euclidean solids”, a term apparently coined by Poincaré Poincare , the reference metric is curved, and the associated geometrical frustration manifests itself through residual stresses Hillel2016 .

A reference curvature can be “embedded” into a solid by using rather well understood techniques of differential swelling, inhomogeneous thermal expansion, bulk growth and remodelling Gladman ; Santangelo2016 ; KupfermanPNAS ; YavariThermal ; C16 . Geometrical frustration can also emerge as a result of surface accretion, as in tree growth, roll winding and dam construction Skalak3 ; Archer ; Correa ; Gibson ; Arut . In the case of surface growth, the relation between the physics of deposition and the resulting incompatibility is implicit and it is not clear which accretion protocol leads to a desired distribution of residual stresses.

Surface growth is often modeled in a holonomic format of elastically coherent phase transitions, that take place without generation of incompatibility Ganghoffer2012 ; Ciarletta2013 . A more general, non-holonomic approach should allow for the incompatibility to be acquired at the moment of the creation of a new continuum particles, see for instance BG ; Trincher ; Gamba . In this Letter we consider a general inverse problem of this type and view the deposition stress as a tensorial control parameter. We find an explicit link between the implemented deposition strategy and the resulting incompatibility. The obtained relations not only reveal the mechanisms of biological adaptation associated with surface growth but can also guide additive manufacturing of programmable meta-materials.

To illustrate the general theory, we study in detail the process of artificial 3D printing of arteries. Presently, the circumferential wrapping of sheets of living cells is used to reproduce their natural layered structure NatureMain2015 ; Lheureux2009 ; Lheureux2012 . However, in physiological conditions the resulting transmural stress distribution is far from realistic Zaucha ; VarnerNelson2014 . Instead, our approach allows reaching the target physiological state precisely, and we show that the proposed strategy is compatible with available manufacturing technologies. As another illustration, we design a growth protocol ensuring that a single cut in a hollow cylinder results in a complete release of the residually stored elastic energy. This prototypical problem is relevant for the understanding of explosive seed dispersal and other functions of plant actuators Hofhuis2016 ; Shahaf2011 ; Xuxu2016 . Quite remarkably, we find that in both biological examples the crucial role is played by a global contribution to incompatibility usually associated with topological defects in crystals (disclinations) Anthony .

Consider a body $\mathscr{B}_{p}$ that in physiological (or working) conditions is subjected to tractions ${\bf s}_{p}$ applied on its boundary $\omega_{p}$. The symmetric physiological stress field $\bm{\sigma}_{p}(\textbf{\emph{x}})$ must satisfy the equilibrium conditions

[TABLE]

where n is the outward normal. In the absence of additional equations, the stress remains under-determined. If the deformation is Euclidean (compatible, neighbor preserving, elastic) the problem can be closed by supplementing (1) with relations expressing the stress in terms of the gradients of vector-valued displacements. This leads to a distribution of stresses that can be altered only by changing the shape of the body or by varying the elastic properties of the material. An alternative way to control the stress state of the body and its deformed shape, apparently favored by biological systems, is to give up the compatibility and employ inelastic deformations.

Suppose (for simplicity) that the Hooke’s law still holds for incremental deformations and define the (linear) elastic strain $\bm{\epsilon}_{p}=\mathbb{C}^{-1}\bm{\sigma}_{p}$, where $\mathbb{C}$ is the elasticity tensor. The signature of the non-Euclidean character of the stress is the nonzero incompatibility $\bm{\eta}_{p}=\text{curl}\,\text{curl}\bm{\epsilon}_{p}$, the linear counterpart of the reference Riemann curvature, which satisfies the Bianchi identities $\text{div}\,\bm{\eta}_{p}={\bf 0}$. If a target incompatibility $\bm{\eta}_{p}$ is prescribed, its three independent components remain unconstrained by (1) and can be used to “engineer” a particular physiological state of stress. Embedding a strain incompatibility into the solid can be viewed as a way of “programming” the material, whose ultimate performances will depend on the loads and the shape of the body, e.g. Danescu .

Note that $\bm{\eta}_{p}$ should be understood in the sense of distributions VGvolterra because the target incompatibility may contain both diffuse and singular contributions. Furthermore, singular defect lines may have a global effect if they carry topological charges characterized by the nonzero Burgers $\textbf{B}_{p}(\textbf{\emph{x}}_{o})=\lim_{h\rightarrow 0}\int_{D_{h}}\textbf{\emph{y}}\times\bm{\eta}_{p}^{\textsf{T}}\textbf{\emph{n}}\,da$ and the Frank $\bm{\Omega}_{p}(\textbf{\emph{x}}_{o})=\lim_{h\rightarrow 0}\int_{D_{h}}\bm{\eta}_{p}^{\textsf{T}}\textbf{\emph{n}}\,da$ vectors Volterra ; Weingarten ; KupfermanPNAS ; SM where y is the position vector of points of an asymptotically shrinking oriented disk $D_{h}$ of diameter $h$, enclosing the singular point $\textbf{\emph{x}}_{o}$ note . Since the associated residual stresses cannot be removed by cutting singular lines out of the body, in non-simply connected bodies such topological charges may be located outside the domain $\mathscr{B}_{p}$.

To model a non-holonomically growing body we introduce a sequence of (incremental) configurations $\mathscr{B}(t)$. The time-like parameter $t$ changes in the interval $(t_{i},t_{f})$, denoting the beginning and the end of the accretion process. In particular, $\mathscr{B}(t_{f})=\mathscr{B}_{p}$. We denote by $\tau(\textbf{\emph{x}})$ the instant when the accreting surface $\omega(t)$, whose evolution is assumed to be known, passes through a point x.

At $t\geq\tau(\textbf{\emph{x}})$ the stress tensor in the growing body is $\bm{\sigma}(\textbf{\emph{x}},t)=\mathring{\bm{\sigma}}(\textbf{\emph{x}})+\int_{\tau(\textbf{\emph{x}})}^{t}\dot{\bm{\sigma}}(\textbf{\emph{x}},s)\,ds$, where the deposition stress $\mathring{\bm{\sigma}}$ can be further decomposed into a sum of two (rank-two) contributions: $\mathring{\bm{\sigma}}(\textbf{\emph{x}})=\textbf{\emph{p}}(\textbf{\emph{x}})+\mathring{\bm{\sigma}}_{a}(\textbf{\emph{x}})$, see Fig.1. While the “passive” contribution $\textbf{\emph{p}}={\bf s}\otimes\textbf{\emph{n}}+\textbf{\emph{n}}\otimes{\bf s}-({\bf s}\cdot\textbf{\emph{n}})\textbf{\emph{n}}\otimes\textbf{\emph{n}}$ is fully defined by the tractions ${\bf s}=\bm{\sigma}\textbf{\emph{n}}$, the “active” surface stress $\mathring{\bm{\sigma}}_{a}$, satisfying $\mathring{\bm{\sigma}}_{a}\textbf{\emph{n}}=\mathring{\bm{\sigma}}_{a}^{\textsf{T}}\textbf{\emph{n}}=\textbf{0}$, carries three independent degrees of freedom that can be used to “implant defects” in the upcoming layers.

The breakdown of $\mathring{\bm{\sigma}}$ into active and passive contributions is somewhat arbitrary, because in both natural and technological conditions the three components of p may also play the role of active agents. In some other cases the implied freedom simply disappears: during coherent structural transformations the whole stress $\mathring{\bm{\sigma}}$ is determined by the evolution laws of the interface Trusk83 ; Ciarletta2013 ; for solidification, it is natural to assume that the active term adjusts to ensure that $\mathring{\bm{\sigma}}$ is hydrostatic KF ; for surface growth in plants and in some additive manufacturing processes, the adhering layer can be treated as a pre-stressed elastic membrane, whose state of deformation is controlled by the requirement of equilibrium under suitable anchoring and adhesion conditions Manz3 ; GorielyKuhl . Here we neglect the potential constraints imposed by the deposition mechanism and we treat the independent components of $\mathring{\bm{\sigma}}_{a}$ as free control parameters.

Since equilibrium must hold at each stage of growth, we must have $\text{div}\,\bm{\sigma}=\textbf{0}$ in $\mathscr{B}(t)$ and $\bm{\sigma}=\mathring{\bm{\sigma}}$ on $\omega(t)$; note that the “whole” stress tensor is prescribed on the growing surface. We assume that inelastic phenomena leading to the accumulation of incompatibilities take place only at the instant of the deposition. Away from the accreting surface the incremental behavior is assumed to be linearly elastic. If we neglect the effect of prestress on the incremental behavior Ho , we can write $\dot{\bm{\sigma}}=\mathbb{C}\dot{\bm{\epsilon}}$, where $\dot{\bm{\epsilon}}=(\nabla\dot{\textbf{\emph{u}}}+\nabla\dot{\textbf{\emph{u}}}^{\textsf{T}})/2$ and $\dot{\textbf{\emph{u}}}$ is the incremental displacement. At a given $t$ we can formulate an incremental problem of elasticity theory in the form Trincher ; SM

[TABLE]

By solving a sequence of such problems we obtain the total elastic strain $\bm{\epsilon}(\textbf{\emph{x}},t)=\mathring{\bm{\epsilon}}(\textbf{\emph{x}})+\int_{\tau(x)}^{t}\dot{\bm{\epsilon}}(\textbf{\emph{x}},s)\,ds$, where the accreting contribution $\mathring{\bm{\epsilon}}=\mathbb{C}^{-1}(\mathring{\bm{\sigma}}_{a}+\textbf{\emph{p}})$, which is also assumed to be small, is generally incompatible. If we now require that the final incompatibility $\bm{\eta}(\textbf{\emph{x}},t_{f})=\text{curl}\,\text{curl}\bm{\epsilon}(\textbf{\emph{x}},t_{f})$ equals its target physiological value $\bm{\eta}_{p}(\textbf{\emph{x}})$, we obtain a constraint on the instantaneous incompatibility of the arriving material $\mathring{\bm{\eta}}=\text{curl}\,\text{curl}\mathring{\bm{\epsilon}}$ in the form

[TABLE]

Here we used the notation $A_{\omega}(\textbf{\emph{x}}):=A(\textbf{\emph{x}},\tau(\textbf{\emph{x}}))$, see SM for details. Since in (3) , $\dot{\bm{\epsilon}}(\textbf{\emph{x}},t)$ implicitly depends on $\mathring{\bm{\sigma}}_{a}(\textbf{\emph{x}})$ through (2), we now have a nonlocal relation between the three independent controls of $\mathring{\bm{\sigma}}_{a}$ and the three independent targets in $\bm{\eta}_{p}$. In the cases when these relations should be understood in the sense of distributions, we implicitly require that in each singular point $\textbf{\emph{x}}_{0}$, $\textbf{B}(\textbf{\emph{x}}_{0},t_{f})=\textbf{B}_{p}(\textbf{\emph{x}}_{0})$ and $\bm{\Omega}(\textbf{\emph{x}}_{0},t_{f})=\bm{\Omega}_{p}(\textbf{\emph{x}}_{0})$. Equation (3), which is the main result of this Letter, defines the deposition strategy that ensures the attainment of a desired stress distribution in physiological conditions.

As an illustration, consider the process of layered manufacturing of an artery NatureMain2015 ; Lheureux2012 . For simplicity, the artery will be modeled as a hollow, infinitely long cylinder loaded in plane strain. This is equivalent to replacing a cylinder by a disk which makes the problem two-dimensional and fully explicit. We assume that the deposition starts on a rigid mandrel of radius $r_{i}$ and that the disk grows outwards till the final (physiological) radius $r_{f}$ is reached. It is convenient to use as a time-like parameter $R=R(t)$, representing the current radius of the accreting surface (line), so that $\tau(R)\equiv R$. Intermediate configurations of the artery are then represented by $r_{i}\leq r\leq R\leq r_{f}$. See Fig.2 for the “macroscopic” rendering of this process: our “microscopic” formulation corresponds to the limit when the thickness of the attached layers $h\to 0$.

If the elastic solid is isotropic and the deposition strategy respects polar symmetry, the incremental displacement reduces to its radial component $\dot{u}(r,R)$. In this case the incremental radial and hoop strains are $\dot{\varepsilon}_{r}=\partial_{r}\dot{u}$ and $\dot{\varepsilon}_{\theta}=\dot{u}/r$, respectively, where the superposed dot denotes $\partial/\partial R$. The incremental stress rates $\dot{\sigma}_{r/\theta}=2\mu\dot{\varepsilon}_{r/\theta}+\lambda(\dot{\varepsilon}_{r}+\dot{\varepsilon}_{\theta})$, where $\lambda$ and $\mu$ are the Lamé moduli, must satisfy the equilibrium equation $\partial_{r}\dot{\sigma}_{r}+\left(\dot{\sigma}_{r}-\dot{\sigma}_{\theta}\right)/r=0$.

Observe that the applied tractions have only radial component $s(R)$, and that the surface component of the deposition stress is fully characterized by its hoop component $\mathring{\sigma}_{a_{\theta}}(R)\sim f_{t}/h$, see Fig.2. Then (2)2 reduces to $\dot{\sigma}_{r}(R,R)=g(R)$ where $g(R)=s^{\prime}(R)+\left(s(R)-\mathring{\sigma}_{a_{\theta}}(R)\right)/R$. Assuming that displacements are fixed on the rigid mandrel, $\dot{u}(r_{i},R)=0$, we obtain an explicit solution of the incremental problem

[TABLE]

To specify the deposition protocol $\mathring{\sigma}_{a_{\theta}}(r)$ we need to satisfy (3) and match the target topological constraints imposed through $\textbf{B}_{p}$ and $\bm{\Omega}_{p}$. In view of our symmetry assumptions, the strain incompatibility tensor reduces to $\bm{\eta}=\eta(r)\textbf{\emph{k}}\otimes\textbf{\emph{k}}$, where the unit vector k is aligned with the cylinder axis, $\eta={\varepsilon_{\theta}}^{\prime\prime}+(2{\varepsilon_{\theta}}^{\prime}-{\varepsilon_{r}}^{\prime})/r$ and ${}^{\prime}=\partial/\partial r$. Eq. (3) reduces to $\mathring{\eta}(R)-(R\,\dot{\varepsilon}_{\theta}(R,R))^{\prime}/R=\eta_{p}(R)$ and conditions on a potential line singularity at $r=0$ take the form $\textbf{B}_{p}=\textbf{0}$ and $\bm{\Omega}_{p}=2\pi\varphi(r_{i})r_{i}\textbf{\emph{k}}$. Here we introduce the function $\varphi(r)={\varepsilon_{\theta}}^{\prime}+({\varepsilon_{\theta}}-\varepsilon_{r})/r$, see SM for additional details. Since $\eta={(\varphi\,r)}^{\prime}/r$, we can recast (3) in the form $\mathring{\varphi}(R)-\dot{\varepsilon}_{\theta}(R,R)=\varphi_{p}(R),$ where $\mathring{\varphi}(r)$ refers to the arriving material and $\varphi_{p}(r)$ to the physiological target state, while $\dot{\varepsilon}_{\theta}(r,R)$ is calculated from (4). Note that if $\varphi_{p}(r_{i})\neq 0$, the target incompatibility has a nonzero topological (global) component.

If we now assume for determinacy that $s(R)=0$ and $\mathring{\sigma}_{a_{\theta}}(r_{i})=0$, we can express the function $\mathring{\varphi}(R)$ in terms of $\varphi_{p}(r)$. This gives the desired deposition strategy securing the attainment of a generic incompatibility:

[TABLE]

For arteries, the physiological state is characterised by a finite internal pressure $p$ acting on $r=r_{i}$ and a much smaller external pressure acting on $r=r_{f}$, which we assume to be equal to zero. Under these conditions, stresses in a purely elastic tube would be transmurally inhomogeneous (Fig.3), which is incompatible with experiments Holzapfel2007 pointing towards homogeneity of the hoop stress ChuFung . To find the physiologically justified incompatibility which guarantees that ${\sigma^{p}_{\theta}}^{\prime}=0$, we combine this target condition with the equilibrium equation ${\sigma_{r}^{p}}^{\prime}+(\sigma_{r}^{p}-\sigma_{\theta}^{p})/r=0$ and obtain that $\sigma^{p}_{r}=-p(r_{f}-r)r_{i}/(r(r_{f}-r_{i}))$ and $\sigma_{\theta}^{p}=p\,r_{i}/(r_{f}-r_{i})$. We can now compute the function $\varphi_{p}(r)$ and substitute it into (5). The resulting deposition strategy

[TABLE]

is illustrated in Fig.3.

Note that the singular component of the incompatibility does not vanish since $\varphi_{p}(r_{i})=p(2\mu+\lambda)r_{f}/\left(4\mu(\mu+\lambda)(r_{f}-r_{i})r_{i}\right)$; the physiological conditions then require the presence of a “ghost” wedge disclination (or its diffuse analog) aligned with the axis of the artery. Since the non-singular part of the incompatibility is also different from zero, the residual stresses cannot be relaxed by a single longitudinal cut turning the cylinder into a simply connected domain. This is consistent with experiments on arteries Holzapfel2007 , showing that the internal layer (media) has a greater opening angle than the external layer (adventitia), see figures in SM . Such behavior is also reproduced by our FEM simulations (Fig.3) for a disk manufactured following the proposed strategy, see SM for details on numerics.

As a second illustration, consider a rather different physiological target that may be relevant for explosive plants Hofhuis2016 . Keeping the same geometry as in the case of arteries, we demand that the distribution of incompatibility is such that the stored elastic energy due to residual stresses is fully released with a single global cut. This requirement will be met if we grow a hollow tube with $\eta_{p}=0$ in the bulk and nonzero $\Omega_{p}$. To this end we must choose $\varphi_{p}(r)=c/r$, with $c$ a constant characterising the magnitude of the stored/released energy. The resulting singular incompatibility field can be interpreted as a Volterra’s wedge disclination with an opening angle $\Omega_{p}=2\pi c$ SM ; VGvolterra ; Anthony . If we now substitute this incompatibility into (5), we obtain that the deposition strategy

[TABLE]

It is illustrated in Fig.4, where we also show by FEM simulation that a single longitudinal slicing of a pre-stressed cylinder with this (purely singular) incompatibility indeed leads to a complete release of the residual stresses, and that subsequent orthogonal slicing does not produce additional relaxation.

To illustrate yet another type of protocols where both tensorial components of the deposition stress, $\textbf{\emph{p}}(\textbf{\emph{x}})$ and $\mathring{\bm{\sigma}}_{a}(\textbf{\emph{x}})$, play an active role, we assume that the newly arriving continuum particles are hydrostatically pre-stressed, with the control parameter $\pi$ representing negative pressure. In the same geometrical setting as above we get $\pi(R)=s(R)=\mathring{\sigma}_{a_{\theta}}(R)$ and, following an almost identical line of reasoning, we obtain that under such boundary/deposition conditions the strain distribution characterized by a generic function $\varphi_{p}(R)$ can be reached if we use the protocol

[TABLE]

Clearly both targets considered above, hoop stress uniformity and a complete release of energy with a single cut, can be achieved in this framework as well. An important example of such hydrostatic “printing” is the crystallisation in a closed container, where the inhomogeneity of the deposition pressure is ensured by the finite compressibility of the melt, e.g. KF .

In conclusion, we outlined a new theoretical framework for controlled incompatible surface growth and obtained explicit relations that can be used to guide additive manufacturing. Acquiring an ability to generate complex patterns of residual stresses is a crucial step in both biological evolution and the design of bio-mimetic meta-materials. The proposed surface deposition strategy promises to bring a combination of unprecedented level of control, together with the ability to handle arbitrarily complex geometries. Future studies are needed to tailor our general theory to specific deposition technologies Kong , to extend it to finite strains SozioYavari , and to develop an energetic framework coupling the velocity of the accretion front with the corresponding driving forces DiCarloSurfBulk .

Acknowledgments. The authors thank M.Destrade, P.Recho and B.Shoikhet for helpful discussions. G.Z. was supported by the ERC Marie Curie Fellowship, INdAM $\&$ GNFM; L.T. was supported by the PSL grant CRITBIO.