In most phenomena it is often necessary to take into account, in addition to the mechanical characteristics, various other quantities. For example, it is important to analyze the variation over time of quantities such as temperature, thermal flows, current intensity, mass flow, concentration gradients, etc.
Most biological tissues can be seen as porous, permeable and deformable media within which fluids such as blood and interstitial fluid infiltrate. Taking these behaviors into account in a computation model, as with linear poroelasticity, allows to simulate and describe their behavior more precisely. Poroelasticity theories are widely used in civil, petroleum and biomedical engineering. There are several facets of this theory, including models based on Biot’s works. The governing equations and the boundary conditions then allow to adapt the computations model to the different environments that can be analyzed.
Table of contents
Background
Porous media theories, or theories that take into account the mechanical and thermodynamic behavior of binary mixtures of immiscible substances, play important roles in many of engineering. The first references to porous media begin around 1800 to see the first formulation of a real theory only in the early twentieth century.
Recent developments in computational mechanics have allowed us to take into account more or less complex phenomena even in simulations. The theory of poroelasticity originates in the works of Paul Fillunger, Karl von Terzaghi and Maurice Biot on the idea of dividing the total tensions between those that fall on the porous skeleton and those that fall on the fluids that fill the pores. By taking into account the constitutive behavior, both of the porous skeleton and of the fluid, the compatibility of the deformations is ensured.
Thermo-chemo-elastic phenomena
In poroelastic materials, therefore, the presence of another substance is considered, in addition to the main material, which influences the mechanical behavior. This is typically measured with concentration, a state variable that adds to the strain and temperature to correctly describe the constitutive behavior. We therefore speak of thermo-chemo-elastic transformations.
The concentration is clearly controllable from the outside and therefore subject to a mass balance that takes into account both the fraction of molecules produced or consumed in the reactions and those flowing through the surface.
\frac{\partial c}{\partial t}=-\operatorname{div}(\mathbf{j})+r_{c}This mass balance adds to the energy balance expressed by the first law of thermodynamics. It is therefore also necessary to consider the presence of a chemical energy which is added to the balance that sees kinetic energy, potential energy, power of external forces and heat exchanged as protagonists. Therefore the chemical energy follows from the localization of the integral of the power, that is of:
\int_{\Omega}\left[\mu r_{c}-\operatorname{div}(\mu \mathbf{j})\right] d v=\int_{\Omega}\left[\mu r_{c}-{\mu \operatorname{div}(\mathbf{j})-\operatorname{grad}(\mu)}{} \cdot \mathbf{j}\right] d v=\int_{\Omega}[\mu \dot{c}-\operatorname{grad}(\mu) \cdot \mathbf{j}] d vSo the first law of thermodynamics can be expressed as:
\rho \dot{e}=\boldsymbol{\sigma}: \dot{\varepsilon}+r-\operatorname{div}(\mathbf{q})+\mu \dot{c}-\operatorname{grad}(\mu) \cdot \mathbf{j}Then we also take into account the second law of thermodynamics through the Clausius-Duhem inequality and the introduction of free energy:
\mathcal{D}=\left(\boldsymbol{\sigma}-\frac{\partial \Psi}{\partial \boldsymbol{\varepsilon}}\right): \dot{\varepsilon}+\left(\mu-\frac{\partial \Psi}{\partial c}\right) \dot{c}-\operatorname{grad}(\mu) \cdot \mathbf{j}-\frac{\mathbf{q} \cdot \operatorname{grad}(T)}{T} \geq 0Uncoupled response
If we consider the two phenomena, chemical and elastic, both present but decoupled, the energy can be expressed as the sum of the two contributions.
\Psi(\varepsilon, c)=\Psi_{e l}(\varepsilon)+\Psi_{c h}(c)\\\:\\
\Psi_{e l}(\varepsilon)=\frac{1}{2}\left(K-\frac{2}{3} G\right)[\operatorname{Tr}(\varepsilon)]^{2}+G \operatorname{Tr}\left(\varepsilon^{T} \varepsilon\right) \\\:\\
\Psi_{c h}(c)=\left(\mu_{0}-R T\right) c+R T c \ln \left(\frac{c}{c_{e q}}\right)
Although the elastic contribution is always that due to the deformation energy, a new contribution is added given by a diffusion force which is proportional to the concentration gradient and inversely proportional to the concentration itself.
\mu=\frac{\partial \Psi}{\partial c}=\mu_{0}+R T \ln \left(\frac{c}{c_{e q}}\right)\\\:\\ \Rightarrow \mathbf{j}=-k_{D} \operatorname{grad}(\mu)=-\frac{R T}{c} \operatorname{grad}(c)Chemomechanical coupling
Then it is possible to consider the concentration variation linked to a porosity variation. In other words, in conditions of poroelasticity, considering the saturated pores, this corresponds to measuring the variation of the fluid content. Therefore it is possible to link the porosity field directly to the concentration through the molar volume and to introduce a new free energy which is justified both thermodynamically and by the constitutive behavior obtainable.
\Psi(\varepsilon, c)=\mu_{0}\left(c-c_{0}\right)+G\left\|\varepsilon_{d e v}\right\|^{2}+\frac{1}{2} K \varepsilon_{v o l}-\left(\alpha M \varepsilon_{v o l}\right) \Omega_{\mathrm{f}}\left(c-c_{0}\right)+\frac{1}{2} M\left[\Omega_{\mathrm{f}}\left(c-c_{0}\right)\right]^{2}Therefore the deformation is divided into volumetric and deviatoric part, where the volumetric contribution couples the two phenomena precisely because of the link between the variation in porosity and that of volume.
\boldsymbol{\sigma}=\frac{\partial \Psi}{\partial \varepsilon}=2 G \varepsilon+\left(K-\frac{2}{3} G\right) \operatorname{Tr}(\varepsilon) \mathbf{I}-\alpha M \Omega_{\mathrm{f}}\left(c-c_{0}\right) \mathbf{I}
\\\:\\
\mu=\frac{\partial \Psi}{\partial c}=\mu_{0}+\Omega_{\mathrm{f}}\left(-\alpha M \varepsilon_{v o l}+M \Omega_{\mathrm{f}}\left(c-c_{0}\right)\right)In other words, the pressure in the pores depends on the porosity through the Biot module ( M ). It is also linked to the volumetric variation always through the Biot module but also for another coefficient, namely the Biot coefficient ( α )
The material behavior clearly shows the mechanical coupling given by the fluid pressure. In fact, considering an applied hydrostatic load, it is possible to consider the two limit cases:
- Conditions undrained , the fluid cannot get out of the pores. Therefore the porosity will not vary. The pressure will be in equilibrium with the internal hydrostatic load and the volume variation will be associated with the hydrostatic field through the Bulk module.
- Conditions drained , the fluid is free to come out. There will be a change in porosity and the pressure in the pores does not change. The change in porosity, proportional to the volumetric deformation, is reflected in a lower Bulk modulus which results in a mechanical response in which the material will be less rigid.
Multiphysics computational models
Having to solve both the balance of forces and mass at the same time, we start again from the free energy and the weak form of the force balance. We then arrive at an analogous form for the mass balance written in weak form:
\frac{\partial c}{\partial t}-\operatorname{Div}\left(k_{c} \nabla \mu\right)=0This will have to be implemented in the FEM formulation.
The weak form depends on some history variables and this will be taken into account considering the incremental derivative. Furthermore, to take into account the mass balance, we start from a new pseudo-energy with a slightly different structure. This pseudenergy, although devoid of a physical sense, allows to use the variational formulation or to calculate a residual that correctly takes into account the balance of forces and mass.
\mathbf{R}=\left.\int_{\Omega_{e}}\left(\frac{\partial \Psi_{\text {pseudo }}}{\partial \mathbf{u}}, \frac{\partial \Psi_{\text {pseudo }}}{\partial c}\right) d \Omega\right|_{b_{c}, \mathbf{b}_{\nabla c}=\mathrm{const}}
\\\:\\ \Psi_{\text {pseudo }}=\frac{1}{2}\left(K-\frac{2}{3} G\right)[\operatorname{Tr}(\varepsilon)]^{2}+G \operatorname{Tr}\left(\varepsilon^{T} \varepsilon\right)\\+
\left(\mu_{o}-R T\right) c+R T c \ln \left(\frac{c}{c_{e q}}\right)\\+
c \underbrace{\left(\frac{\partial c}{\partial t}\right)}_{b_{c}}+\nabla c \cdot \underbrace{\left(k_{c} \nabla \mu\right)}_{\mathbf{b} \nabla c}
By implementing this formulation, a multi-field response is obtained but with the two fields decoupled. This can also be seen by obtaining the tangent matrix which will in fact be symmetrical and made up of two non-zero blocks representing the governing equations of the force balance and the mass balance.
Linear poroelasticity
By combining the mass balance with a constitutive law that encloses the variation of the fluid content by binding it to the pressure in the pores, an equation is obtained that expresses the variation of the pressure in the pores in space and time. This variation is in turn linked to the volumetric deformation and is referred to as the poroelasticity regime.
Then considering the presence of a link between pressure and deformation this is expressed by the constitutive link.
\left\{\begin{array}{l}
\boldsymbol{\sigma}=\boldsymbol{\sigma}_{\text {eff }}+\sigma_{p} \mathbf{I}=\boldsymbol{\sigma}_{\text {eff }}-\alpha p \mathbf{I} \\
\mathbf{q}=-k_{p} \nabla p \\
p=M\left(\zeta_{f}-\alpha \varepsilon_{v o l}\right)
\end{array} \right.Therefore the variation formulation can be expressed starting from a pseudo-energy, defined as:
\Psi_{p s e u d o}=\Psi_{e l}+\Psi_{p m}+\Psi_{p f}That is, to which a new term has been added:
\Psi_{p f}(p, \varepsilon):=p \underbrace{\left(\alpha \frac{\partial \varepsilon_{v o l}}{\partial t}+\frac{1}{M} \frac{\partial p}{\partial t}\right)}_{b_{p}}+\nabla p \cdot \underbrace{\left(k_{p} \nabla p\right)}_{\mathbf{b}_{\nabla p}}By discretizing the unknown functions both in terms of displacements and pressure, it is possible to write the residue as:
\mathbf{R}=\left.\int_{\Omega_{e}}\left(\frac{\partial \Psi_{p s e u d o}}{\partial \mathbf{u}}, \frac{\partial \Psi_{p s e u d o}}{\partial p}\right) d \Omega\right|_{b_{p}, \mathbf{b}_{\nabla p}, b_{v o l}=\mathrm{const}}Chemo-elastic applications
Considering the presence of an elastic deformation field and the presence of a diffusion force field but considering them as independent fields, it is possible to analyze a chemoelastic problem. Let us consider 1/4 of a cylinder and appropriate conditions at the symmetry edge that make the simulation equivalent to the simulation of a hollow cylinder with a defined pressure inside it.
This example is typical of most biomedical applications such as tissues or vessel walls where hollow cylindrical structures are subjected to the passage of a fluid, diffusion of solutes and pressure gradients that influence their behavior, even in the elastic regime.
We therefore consider 1/4 of a circular crown as the reference volume with a pressure on the inner edge.
Pressure field
The internal pressure is applied through a load multiplier which allows an initial gradual increase as shown in the figure. The resolution then addresses incremental load steps.
The internal pressure has two important consequences: it widens the cylinder and tends to push the fluid towards the outside. Since there is a zero flow at the outer edge, we will see an increase in concentration at the outer edge due to the thrust of the internal pressure field. There is a first phase in which there is an increase in concentration linked to the load condition. Subsequently, the pressure at the inner edge reaches its maximum value and this is reflected in a progressive increase in concentration, point by point, over time.
Clearly the concentration, which increases over time, will also vary in the thickness of the wall and you can see how it descends moving along the ray to end up with a point with a horizontal tangent. This reflects precisely the fact that the flow (its derivative) at the boundary is zero.
Decoupling
The widening of the circular section stabilizes at the second load step where the solution of the elastic equilibrium problem leads to the energy minimum. In other words, the condition of equilibrium is reached with the pressure load applied.
Although this equilibrium satisfies the force balance it does not satisfy the mass balance which results in a progressive increase in concentrations at the edge. That is, there is a continuous mass transfer between the inner edge, where the load condition is applied, and the outer edge which is at lower pressures.
It is therefore strongly evident, not being in the poroelasticity regime, the decoupling between the force balance and the mass balance. Initially the pressure leads to a deformation which is reflected in the satisfaction of the force balance. However, the mass balance remains unsatisfied, which is reflected in a progressive increase in concentration at the edge.
Poroelasticity and coupling of the response
Considering instead a linear poroelasticity regime, where the concentration is associated with a reference porosity, the volumetric deformations are reflected in variations of the concentration and vice versa. In other words, the coupling of the response is reflected precisely in the fact that the variations in porosity can be due both to volumetric deformations and to a pressure load in the pores.
\Omega_{\mathrm{f}}\left(c-c_{0}\right)=\zeta=\alpha \varepsilon_{v o l}+\frac{p}{M}This coupling reflects an elastic behavior influenced by the presence and behavior of the fluid. That is, the elastic behavior is affected by the possibility of the fluid to flow and to undergo the effects of the pressure field.
Implementation of poroelasticity
By numerically implementing the poroelasticity in an FEM formulation, the first coupling results appear from the matrix of coefficients K. This tangent matrix loses its symmetry due to the mixed contributions linking displacements and pressure field.
\left[\begin{array}{c|c }
\bold K_{uu} &\bold K_{up}\\
\hline
\bold K_{pu} &\bold K_{pp}
\end{array}\right]Punch test
To analyze this behavior, you can consider a punch test. That is, we consider a rectangular reference volume on which a certain voltage is applied.
The applied load tends to compress the reference volume which will respond with a deformation in the elastic range. It is interesting to observe how this deformation is affected by the pressure field.
Without applying particular conditions to the edge, there would be a zero flow over the entire perimeter. That is, it is in non-draining conditions. Then the pressure will be in equilibrium with the internal hydrostatic load and the volume variation will be associated with the hydrostatic field through the Bulk module.
The pressure in the central point, under the applied load, initially increases up to a maximum value where a first equilibrium is reached between the internal load and the elastic deformation. However, this effect is also coupled with the pressure in the pores which, having increased, will tend to redistribute the fluid. The fluid, even if it cannot actually exit the edges, can redistribute itself within the reference volume.
The pressure increase due to the load will tend to empty the part immediately under the load (where the pressure has increased considerably) and to distribute part of the fluid in the lateral areas. This is also reflected in a further lowering of the loading area until hydrostatic equilibrium is reached.
Draining conditions
In the poroelasticity regime it is also possible to analyze drained conditions where the fluid is also allowed to pass through the edge surfaces. Clearly, depending on the applied pressure, there will be different equilibrium conditions. Starting from the lower edge it is possible to analyze what happens when the pressure changes. It is therefore possible to implement a simulation campaign that analyzes 9 different cases, between -250 and 250, also passing through the undrained conditions.
tempVector = {};
pressureBC = {-250, -200, -100, -50, 0, 50, 100, 200, 250};
myPressureTotal = Table[0, {i, 1, Length[pressureBC]}];
myDisplacementTotal = Table[0, {i, 1, Length[pressureBC]}];
displacementMax = Table[0, {i, 1, Length[pressureBC]}];
pressureMax = Table[0, {i, 1, Length[pressureBC]}];
Do[
SetupSimulation[pressureBC[[i]], "drainBorder" -> "down"];
Print["pressure BC: ", pressureBC[[i]]];
tempVector = SolveSimulation[];
myPressureTotal[[i]] = tempVector[[1]];
myDisplacementTotal[[i]] = tempVector[[2]];
displacementMax[[i]] = Min[myDisplacementTotal[[i, All, 2]]];
pressureMax[[i]] = Max[myPressureTotal[[i, All, 2]]];
Print[Show[SMTShowMesh["FillElements" -> False, ImageSize -> 400],
SMTShowMesh["Field" -> "p", "DeformedMesh" -> True,
ImageSize -> 400]],
Show[SMTShowMesh["FillElements" -> False, ImageSize -> 400],
SMTShowMesh["Field" -> "v", "DeformedMesh" -> True,
ImageSize -> 400]]]
(*If[i==(Length[pressureBC]+1)/2,SMTAnimationOfResponse["Export \
to","mp4","FirstFrame"->"LastFrame","FileName"->"punch_test_null_BC",\
FrameRate->2]];*)
, {i, 1, Length[pressureBC]}]
Print[GraphicsRow[{ListPlot[Transpose[{pressureBC, displacementMax}],
PlotMarkers -> {"OpenMarkers"},
PlotLabel -> "min displacement"],
ListPlot[Transpose[{pressureBC, pressureMax}],
PlotMarkers -> {"OpenMarkers"}, PlotLabel -> "max pressure"]}]];
Print[GraphicsRow[{ListLinePlot[myDisplacementTotal,
PlotLabels -> pressureBC, PlotLabel -> "displacement"],
ListLinePlot[myPressureTotal, PlotLabels -> pressureBC,
PlotLabel -> "pressure"]}]];Different pressure conditions leading to different equilibrium conditions. This is easily observed by analyzing the displacements of the loading area or the hydrostatic balance. As the pressure at the edge increases, there is an increasingly smaller lowering of the loading area.
This is seen even better over time. As the pressure gradient at the edge changes, the equilibrium points move more and more and in fact higher and higher maximum pressures are reached. The trend over time remains similar but is translated to ever smaller drops by increasing the pressure gradient at the edge. This consideration is true up to a certain limit in which the conditions are reversed and a hydrostatic equilibrium is reached at pressures higher than what the external load alone would generate. That is, it is even reflected in a relative final elevation of the loading area. That is, at high edge pressures the reference volume tends to swell due to the high flow rate which is no longer rejected by the external load.
This phenomenon is even more evident by observing the deformations and the pressure field with different pressure conditions at the edge.
Variations in symmetry
In this case, in addition to the strong poroelasticity coupling, a perfect condition of horizontal symmetry is evident. This symmetry disappears by choosing non-symmetrical load conditions, for example by introducing a drained condition only on the right edge.
The trend of the lowering and pressure balance follows a logic similar to the previous case, although the differences become much smaller. That is, it is due to the displacement of the hydrostatic equilibrium on the right edge, reducing the effect of the external load. In other words, in the central position, the material behavior is mainly affected by the compressive load applied. Similarly, the left edge presents conditions very similar to an undrained case. On the right edge, on the other hand, the effect of the pressure gradient is evident, which tends to inflate or deflate the reference volume depending on its direction.
The loss of horizontal symmetry is also evident by observing the deformations and the pressure field with different pressure conditions at the edge.
An overall pressure gradient is therefore created mainly horizontal which however is affected by the load, more or less strongly, in the central area where the latter is applied.
Clearly this situation changes if conditions are applied to both edges. Applying discordant conditions will simply combine the negative and positive pressure conditions but on two opposite sides.
Agreeing conditions, on the other hand, will amplify the effects of the draining conditions on the single side while maintaining symmetry conditions.
Pressure gradient and external load
It is interesting to note that even applying very high positive pressures, that is a condition for which one could think of re-lifting the load, the vertical displacement of the loading area remains almost constant. This is due to the fact that there is linear elasticity on the whole structure while the loading area, due to the external tension, becomes more rigid. that is, once the applied load has been balanced, the overpressure pushes on the edges, inflating or deflating the rest of the structure but leaving the loading area almost unchanged.
Decoupling
Finally, it is possible to go back to studying a case in which there is a decoupling between elasticity and pressure field, binding the pressure to be zero in each node.
If[
OptionValue["uncoupled"] == True,
SMTAddEssentialBoundary["p", 1 -> 0]
];This leads to an elastic problem that is not affected in the least by the pressure gradient but only by the applied load and linear elasticity.
Code
The code is available at the GitHub repository: https://github.com/mastroalex/poroelasticity
References
- Meccanica Computazionale dei Tessuti e Biomateriali – Ingegneria medica – Università degli studi di Roma Tor Vergata
- A multiple-network poroelastic model for biological systems and application to subject-specific modelling of cerebral fluid transport – Guo, L. et al
- Development of porous media theories — A brief historical review – de Boer
- Poroelasticity – Cheng, A
