Nella maggior parte dei fenomeni è spesso necessario portare in conto, oltre alle caratteristiche meccaniche, diverse altre grandezza. Ad esempio è importante analizzare la variazione nel tempo di grandezze come la temperatura, flussi termici, intensità di corrente, flusso di massa, gradienti di concentrazione, ecc.
La maggior tessuti biologici possono essere visti come mezzi porosi, permeabili e deformabili al cui interno si infiltrano fluidi come sangue e liquido interstiziale. Portare in conto questi comportamenti in un modello computazione, come con la poroelasticità lineare, permette di simulare e descrivere in modo più preciso il loro comportamento. Teorie della poroelasticità sono largamente utilizzare nell’ingegneria civile, petrolifera e biomedica. Esistono diverse sfaccettature di questa teoria tra le quali i modelli basati sui lavori di Biot. Le equazioni di governo e le condizioni al contorno permettono poi di adattare il modello computazioni ai diversi ambienti che possono essere analizzati.
Indice
Background
Teorie dei mezzi porosi, ovvero teorie che portano in conto il comportamento meccanico e termodinamico di miscele binarie di sostanze immiscibili, svolgono ruoli importanti in molti campi dell’ingegneria. I primi riferimenti a mezzi porosi iniziano intorno al 1800 per vedere la prima formulazione di una teoria vera e propria solamente nei primi anni del ventesimo secolo.
I recenti sviluppi nella meccanica computazionale hanno permesso di portare in conto fenomeni più o meno complessi anche nelle simulazioni. La teoria della poroelasticità trova origine nei lavori di Paul Fillunger, Karl von Terzaghi e Maurice Biot basati sull’idea di ripartire le tensioni totali tra quelle che ricadono sullo scheletro poroso e quelle che invece ricadono sui fluido che riempie i pori. Portando in conto il comportamento costitutivo, sia dello scheletro poroso che del fluido, si assicura la compatibilità delle deformazioni.
Processi termo chemo elastici
Nei materiali poroelastici quindi viene considerata la presenza di un’altra sostanze, oltre al materiale principale, che influenza il comportamento meccanico. Questa viene tipicamente misurata con la concentrazione, variabile di stato che va ad aggiungersi alla deformazione e alla temperatura per descrivere correttamente il comportamento costitutivo. Si parla quindi di trasformazioni termo-chemo-elastiche.
La concentrazione è chiaramente controllabile dall’esterno per cui soggetta ad un bilancio di massa che tiene conto sia della frazione di molecole prodotte o consumate nelle reazioni sia di quelle che fluiscono attraverso la superficie.
\frac{\partial c}{\partial t}=-\operatorname{div}(\mathbf{j})+r_{c}Questo bilancio di massa va ad aggiungersi al bilancio di energia espresso dal primo principio della termodinamica. È quindi necessario considerare anche la presenza di un’energia chimica che va ad aggiungersi al bilancio che vede protagonisti energia cinetica, energia potenziale, potenza delle forze esterne e calore scambiato. Dunque l’energia chimica segue dalla localizzazione dell’integrale della potenza, ovvero di:
\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 vQuindi la prima legge della termodinamica può essere espressa come:
\rho \dot{e}=\boldsymbol{\sigma}: \dot{\varepsilon}+r-\operatorname{div}(\mathbf{q})+\mu \dot{c}-\operatorname{grad}(\mu) \cdot \mathbf{j}Allora si porta in conto anche il secondo principio della termodinamica tramite la disuguaglianza di Clausius-Duhem e l’introduzione dell’energia libera:
\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 0Risposta disaccoppiata
Se si considerano i due fenomeni, chimico ed elastico, entrambi presenti ma disaccoppiati, l’energia può essere espressa come la somma dei due contributi.
\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)
Sebbene il contributo elastico sia sempre quello dovuto all’energia di deformazione, si aggiunge un nuovo contributo dato da una forza di diffusione che risulta proporzionale al gradiente di concentrazione e inversamente proporzionale alla concentrazione stessa.
\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)Accoppiamento chemomeccanico
Allora è possibile considerare la variazione di concentrazione legata ad una variazione di porosità. Ovvero in condizioni di poroelasticità, considerando i pori saturi, questo corrisponde a misurare la variazione del contenuto fluido. Quindi è possibile legare il campo di porosità direttamente alla concentrazione tramite il volume molare ed introdurre una nuova energia libera che risulta giustificata sia termodinamicamente che dal comportamento costitutivo ottenibile.
\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}Pertanto la deformazione viene divisa in parte volumetrica e deviatorica, dove il contributo volumetrico va ad accoppiare i due fenomeni proprio a causa del legame tra la variazione di porosità e quella di 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)Ovvero la pressione nei pori dipende dalla porosità per il tramite del modulo di Biot (M). È legata anche alla variazione volumetrica sempre per il tramite del modulo di Biot ma anche per un altro coefficiente, ovvero il coefficiente di Biot (α).
Il comportamento materiale presenta in modo evidente l’accoppiamento meccanico dato dalla pressione del fluido. Considerando infatti un carico idrostatico applicato è possibile considerare i due casi limite:
- Condizioni undrained, il fluido non può uscire dai pori. Quindi non varierà la porosità. La pressione sarà in equilibrio con il carico idrostatico interno e la variazione di volume sarà associata al campo idrostatico per il tramite del modulo di Bulk.
- Condizioni drained, il fluido è libero di uscire. Ci sarà un cambio di porosità e la pressione nei pori non cambia. Il cambio di porosità, proporzionale alla deformazione volumetrica, si rispecchia in un modulo di Bulk minore che si traduce in una risposta meccanica in cui il materiale sarà meno rigido.
Modelli computazionali multifisici
Dovendo risolvere contemporaneamente sia il bilancio di forze che di massa si riparte dall’energia libera e dalla forma debole del bilancio di forza. Si arriva quindi ad una forma analoga per il bilancio di massa scritto in forma debole:
\frac{\partial c}{\partial t}-\operatorname{Div}\left(k_{c} \nabla \mu\right)=0Questo dovrà essere implementato nella formulazione FEM.
La forma debole dipende da alcune variabili di storia e questo verrà portato in conto considerandone la derivata incrementale. Inoltre, per portare in conto del bilancio di massa, si parte da una nuova pseudoenergia con una struttura leggermente diversa. Questa pseudenergia, sebbene sia priva di un senso fisico, permette di utilizzare la formulazione variazionale ovvero arrivare al calcolo di un residuo che porta in conto correttamente i bilanci di forze e di massa.
\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}
Implementando questa formulazione si ottiene una risposta multicampo ma con i due campi disaccoppiate. Questo lo si vede anche ricavando la matrice tangente che sarà di fatto simmetrica e costituita da due blocchi non nulli che rappresentano l’equazioni di governo del bilancio di forze e del bilancio di massa.
Poroelasticità lineare
Unendo il bilancio di massa ad una legge costitutiva che racchiude la variazione del contenuto fluido legandolo alla pressione nei pori si ottiene un’equazione che esprime la variazione della pressione nei pori nello spazio e nel tempo. Questa variazione è a sua volta legata alla deformazione volumetrica e si parla di regime di poroelasticità.
Allora considerando la presenza di un legame tra la pressione e la deformazione questo viene espresso dal legame costitutivo.
\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.Quindi la formulazione variazione è esprimibile a partire da una pseudoenergia, definita come:
\Psi_{p s e u d o}=\Psi_{e l}+\Psi_{p m}+\Psi_{p f}Ovvero alla quale si è aggiunto un nuovo termine:
\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}}Discretizzando le funzioni incognite sia in termini di spostamenti che di pressione è possibile scrivere il residuo come:
\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}}Applicazioni chemo-elastiche
Considerando la presenza di un campo di deformazioni elastiche e la presenza di un campo di forze di diffusione ma considerandoli come campi indipendenti è possibile analizzare un problema chemoelastico. Consideriamo 1/4 di cilindro e opportune condizioni al bordo di simmetria che rendono la simulazione equivalente alla simulazione di un cilindro cavo con al suo interno una pressione definita.
Questo esempio è tipico di gran parte delle applicazioni biomedicali come tessuti o pareti dei vasi dove strutture cilindriche cave sono sottoposte al passaggio di un fluido, a diffusione dei soluti e a gradienti pressori che ne influenzano il comportamento, anche nel regime elastico.
Consideriamo quindi come volume di riferimento 1/4 di corona circolare con una pressione sul bordo interno.
La pressione viene interna viene applicata tramite un moltiplicatore di carico che ne permette un iniziale aumento graduale come rappresentato in figura. La risoluzione affronta quindi passi di carico incrementali.
La pressione interna ha due importanti conseguenze: allarga il cilindro e tende a spingere il fluido verso la parte esterna. Essendoci un flusso nullo al bordo esterno vedremo un aumento della concentrazione al bordo esterno a causa della spinta del campo pressorio interno. È presente una prima fase in cui c’è un aumento di concentrazione legato alla condizione di carico. Successivamente la pressione al bordo interno si porta al valore massimo e questo si rispecchia in un aumento progressivo della concentrazione, punto per punto, nel passare del tempo.
Chiaramente la concentrazione, che cresce nel tempo, varierà anche nello spessore della parete e si vede come scende spostandosi lungo il raggio per finire con un punto a tangente orizzontale. Questo rispecchia proprio il fatto che il flusso (sua derivata) al bordo è nullo.
Disaccoppiamento
L’allargamento della sezione circolare si stabilizza al secondo passo di carico dove la soluzione del problema di equilibrio elastico porta al minimo energetico. Ovvero si arriva alla condizione di equilibrio con il carico pressorio applicato.
Sebbene questo equilibrio soddisfi il bilancio di forze non soddisfa il bilancio di massa che si traduce in un progressivo aumento della concentrazioni al bordo. Ovvero c’è un continuo trasferimento di massa tra il bordo interno, dove è applicata la condizione di carico, e il bordo esterno che si trova a pressioni inferiori.
È quindi fortemente evidente, non essendo in regime di poroelasticità, il disaccoppiamento tra il bilancio di forze e il bilancio di massa. Inizialmente la pressione porta ad una deformazione che si rispecchia nel soddisfacimento del bilancio di forze. Rimane però insoddisfatto il bilancio di massa che si rispecchia in un progressivo aumento della concentrazione al bordo.
Poroelasticità e accoppiamento della risposta
Considerando invece un regime di poroelasticità lineare, dove la concentrazione è associata ad una porosità di riferimento, le deformazioni volumetriche si riflettono in variazioni della concentrazione e viceversa. Ovvero l’accoppiamento della risposta si rispecchia proprio nel fatto che le variazioni di porosità possono essere dovute sia alle deformazioni volumetriche che ad un carico pressorio nei pori.
\Omega_{\mathrm{f}}\left(c-c_{0}\right)=\zeta=\alpha \varepsilon_{v o l}+\frac{p}{M}Questo accoppiamento si rispecchia un comportamento elastico influenzato dalla presenza e dal comportamento del fluido. Ovvero il comportamento elastico risente della possibilità del fluido di scorrere e di subire gli effetti del campo pressorio.
Implementazione della poroelasticità
Implementando numericamente la poroelasticità in una formulazione FEM i primi risultati dell’accoppiamento appaiono dalla matrice dei coefficienti K. Tale matrice tangente perde la sua simmetria a causa dei contributi misti che legano spostamenti e campo pressorio.
\left[\begin{array}{c|c }
\bold K_{uu} &\bold K_{up}\\
\hline
\bold K_{pu} &\bold K_{pp}
\end{array}\right]Punch test
Per analizzare questo comportamento è possibile considerare un punch test. Ovvero consideriamo un volume di riferimento rettangolare sul quale viene applicata una certa tensione.
Il carico applicato tende a comprimere il volume di riferimento che risponderà con una deformazione in campo elastico. È interessante osservare come questa deformazione risente del campo pressorio.
Senza applicare particolari condizioni al bordo si avrebbe un flusso nullo su tutto il perimetro. Ovvero si è in condizioni non drenanti. Allora la pressione sarà in equilibrio con il carico idrostatico interno e la variazione di volume sarà associata al campo idrostatico per il tramite del modulo di Bulk.
La pressione nel punto centrale, sotto il carico applicato aumenta inizialmente fino ad un valore massimo dove si arriva ad un primo equilibrio tra il carico interno e la deformazione elastica. A questo effetto però si accoppia anche la pressione nei pori che, essendo aumentata, tenderà a ridistribuire il fluido. Il fluido, anche se non può effettivamente uscire dai bordi, può ridistribuirsi all’interno del volume di riferimento.
L’aumento di pressione dovuto al carico tenderà a far svuotare la parte immediatamente sotto il carico (dove è aumentata notevolmente la pressione) e a far distribuire parte del fluido nella zone laterali. Questo si rispecchia anche in un ulteriore abbassamento della zona di carico fino al raggiungimento dell’equilibrio idrostatico.
Condizioni drenanti
In regime di poroelasticità è possibile analizzare anche condizioni di tipo drained dove al fluido è concesso anche di passare tramite le superfici di bordo. Chiaramente a seconda della pressione applicata si avranno differenti condizioni di equilibrio. Partendo dal bordo inferiore è possibile analizzare cosa succede al variare della pressione. È possibile quindi implementare una campagna di simulazione che analizza 9 differenti casi, tra -250 e 250, passando anche per le condizioni undrained.
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"]}]];Differenti condizioni pressorie portando a differenti condizioni di equilibrio. Questo è facilmente osservabile analizzando gli spostamenti della zona di carico o l’equilibrio idrostatico. Al crescere della pressione al bordo si osserva un abbassamento della zona di carico sempre minore.
Questo lo si vede ancora meglio nel tempo. Al variare del gradiente pressorio al bordo si sposta sempre più il punti di equilibrio e infatti si raggiungono pressioni massime sempre più alte. L’andamento nel tempo rimane simile ma viene traslato ad abbassamenti sempre minori aumentando il gradiente pressorio al bordo. Tale considerazione è vera fino ad un certo limite in cui le condizioni si invertono e si arriva ad un equilibrio idrostatico a pressioni maggiori di quanto andrebbe a generare il solo carico esterno. Ovvero si rispecchia addirittura in un relativo innalzamento finale della zona di carico. Ovvero a pressioni al bordo elevate il volume di riferimento tende a gonfiarsi a causa dell’elevato introito di flusso che non viene più respinto dal carico esterno.
Questo fenomeno è ancora più evidente osservando le deformate e il campo di pressione con differenti condizioni di pressione al bordo.
Variazioni della simmetria
In questo caso è evidente, oltre al forte accoppiamento di poroelasticità, una perfetta condizione di simmetria orizzontale. Tale simmetria viene meno andando a scegliere condizioni di carico non simmetriche, ad esempio andando ad introdurre una condizione drained soltanto sul bordo destro.
L’andamento dell’abbassamento e dell’equilibrio pressorio segue una logica analoga al caso precedente sebbene le differenze diventano molto minori. Cioè è dovuto allo spostamento dell’equilibrio idrostatico sul bordo destro andando a ridurre l’effetto del carico esterno. Ovvero nella posizione centrale il comportamento materiale risente principalmente del carico di compressione applicato. Allo stesso modo il bordo sinistro presenta delle condizioni molto simili ad un caso undrained. Al bordo destro invece è evidente l’effetto del gradiente pressorio che tende a gonfiare o sgonfiare il volume di riferimento a seconda della sua direzione.
È evidente anche la perdita di simmetria orizzontale osservando le deformate e il campo di pressione con differenti condizioni di pressione al bordo.
Si crea quindi un gradiente pressorio complessivo principalmente orizzontale che però risente del carico, più o meno fortemente, nella zona centrale dove quest’ultimo è applicato.
Chiaramente questa situazione cambia se le condizioni vengono applicati ad entrambi i bordi. Applicando condizioni discordi si andrà semplicemente ad unire la condizioni di pressione negativa e positiva ma su due lati opposti.
Condizioni concordi invece andranno ad amplificare gli effetti delle condizioni drenanti sul singolo lato mantenendo però condizioni di simmetria.
Gradiente pressorio e carico esterno
È interessante notare come anche applicando pressioni positive molto alte, ovvero condizione per la quale si potrebbe pensare di ri-sollevare il carico, lo spostamento verticale della zona di carico rimane quasi costante. Questo è dovuto al fatto che c’è elasticità lineare su tutta la struttura mentre la zona di carico, a causa della tensione esterna, diventa più rigida. ovvero, una volta bilanciato il carico applicato la sovrapressione va a spingere sui bordi gonfiando, o sgonfiando, il resto della struttura ma lasciando la zona di carico quasi inalterata.
Disaccoppiamento
Infine è possibile tornare a studiare un caso in cui c’è disaccoppiamento tra elasticità e campo di pressioni andando a vincolare la pressione ad essere nulla in ogni nodo.
If[
OptionValue["uncoupled"] == True,
SMTAddEssentialBoundary["p", 1 -> 0]
];Questo porta ad un problema elastico che non risente minimamente del gradiente pressorio ma soltando del carico applicato e dell’elasticità lineare.
Codice
Il codice è disponibile alla repository GitHub: https://github.com/mastroalex/poroelasticity
Fonti
- 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
