#include <deal.II/base/work_stream.h>\r
#include <deal.II/base/mpi.h>\r
#include <deal.II/base/quadrature_point_data.h>\r
-#include <deal.II/base/std_cxx11/shared_ptr.h>\r
\r
#include <deal.II/differentiation/ad.h>\r
\r
+#include <deal.II/distributed/shared_tria.h>\r
+\r
#include <deal.II/dofs/dof_renumbering.h>\r
#include <deal.II/dofs/dof_tools.h>\r
#include <deal.II/dofs/dof_accessor.h>\r
#include <deal.II/grid/grid_in.h>\r
#include <deal.II/grid/grid_out.h>\r
#include <deal.II/grid/manifold_lib.h>\r
-#include <deal.II/grid/tria.h>\r
#include <deal.II/grid/tria_accessor.h>\r
#include <deal.II/grid/tria_boundary_lib.h>\r
#include <deal.II/grid/tria_iterator.h>\r
gravity_value = prm.get_double("gravity value");\r
\r
if ( (fluid_type == "Markert") && ((init_intrinsic_perm == 0.0) || (viscosity_FR == 0.0)) )\r
- throw std::runtime_error ("Markert seepage velocity formulation requires the definition of "\r
- "'initial intrinsic permeability' and 'fluid viscosity' greater than 0.0.");\r
+ AssertThrow(false, ExcMessage("Markert seepage velocity formulation requires the definition of "\r
+ "'initial intrinsic permeability' and 'fluid viscosity' greater than 0.0."));\r
\r
if ( (fluid_type == "Ehlers") && ((init_darcy_coef == 0.0) || (weight_FR == 0.0)) )\r
- throw std::runtime_error ("Ehler seepage velocity formulation requires the definition of "\r
- "'initial Darcy coefficient' and 'fluid weight' greater than 0.0.");\r
+ AssertThrow(false, ExcMessage("Ehler seepage velocity formulation requires the definition of "\r
+ "'initial Darcy coefficient' and 'fluid weight' greater than 0.0."));\r
\r
const std::string eigen_solver_type = prm.get("eigen solver");\r
if (eigen_solver_type == "QL Implicit Shifts")\r
}\r
iteration += 1;\r
if (iteration > 15 )\r
- throw std::runtime_error ("No convergence in local Newton iteration for the "\r
- "viscoelastic exponential time integration algorithm.");\r
+ AssertThrow(false, ExcMessage("No convergence in local Newton iteration for the "\r
+ "viscoelastic exponential time integration algorithm."));\r
}\r
\r
NumberType aux_J_e_1 = 1.0;\r
permeability_term = get_darcy_flow_current(F) / weight_FR;\r
\r
else\r
- throw std::runtime_error ("Material_Darcy_Fluid --> Only Markert "\r
- "and Ehlers formulations have been implemented.");\r
+ AssertThrow(false, ExcMessage(\r
+ "Material_Darcy_Fluid --> Only Markert "\r
+ "and Ehlers formulations have been implemented."));\r
\r
return ( -1.0 * permeability_term * det_F\r
* (grad_p_fluid - get_body_force_FR_current()) );\r
seepage_velocity = get_seepage_velocity_current(F,grad_p_fluid);\r
}\r
else\r
- throw std::runtime_error ("Material_Darcy_Fluid --> Only Markert and Ehlers formulations have been implemented.");\r
+ AssertThrow(false, ExcMessage(\r
+ "Material_Darcy_Fluid --> Only Markert and Ehlers "\r
+ "formulations have been implemented."));\r
\r
dissipation_term = ( invert(permeability_term) * seepage_velocity ) * seepage_velocity;\r
dissipation_term *= 1.0/(det_F*det_F);\r
return body_force;\r
}\r
private:\r
- std_cxx11::shared_ptr< Material_Hyperelastic<dim, NumberType> > solid_material;\r
- std_cxx11::shared_ptr< Material_Darcy_Fluid<dim, NumberType> > fluid_material;\r
+ std::shared_ptr< Material_Hyperelastic<dim, NumberType> > solid_material;\r
+ std::shared_ptr< Material_Darcy_Fluid<dim, NumberType> > fluid_material;\r
};\r
\r
// @sect3{Nonlinear poro-viscoelastic solid}\r
void run();\r
\r
protected:\r
- typedef Sacado::Fad::DFad<double> ADNumberType;\r
+ using ADNumberType = Sacado::Fad::DFad<double>;\r
\r
std::ofstream outfile;\r
std::ofstream pointfile;\r
const unsigned int this_mpi_process;\r
mutable ConditionalOStream pcout;\r
\r
- //A collection of the parameters used to describe the problem setup\r
+ // A collection of the parameters used to describe the problem setup\r
const Parameters::AllParameters ¶meters;\r
\r
- //Declare an instance of dealii Triangulation class (mesh)\r
- Triangulation<dim> triangulation;\r
+ // Declare an instance of dealii Triangulation class (mesh)\r
+ parallel::shared::Triangulation<dim> triangulation;\r
\r
// Keep track of the current time and the time spent evaluating certain functions\r
Time time;\r
this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)),\r
pcout(std::cout, this_mpi_process == 0),\r
parameters(parameters),\r
- triangulation(Triangulation<dim>::maximum_smoothing),\r
+ triangulation(mpi_communicator,Triangulation<dim>::maximum_smoothing),\r
time(parameters.end_time, parameters.delta_t),\r
timerconsole( mpi_communicator,\r
pcout,\r
timerconsole.enter_subsection("Setup system");\r
timerfile.enter_subsection("Setup system");\r
\r
- // Partition triangulation\r
- GridTools::partition_triangulation (n_mpi_processes, triangulation);\r
-\r
-\r
//Determine number of components per block\r
std::vector<unsigned int> block_component(n_components, u_block);\r
block_component[p_fluid_component] = p_fluid_block;\r
\r
//Setup the initial quadrature point data using the info stored in parameters\r
FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>\r
- cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+ cell (IteratorFilters::LocallyOwnedCell(),\r
dof_handler_ref.begin_active()),\r
- endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+ endc (IteratorFilters::LocallyOwnedCell(),\r
dof_handler_ref.end());\r
for (; cell!=endc; ++cell)\r
{\r
- Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());\r
+ Assert(cell->is_locally_owned(), ExcInternalError());\r
+ Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());\r
+\r
const std::vector<std::shared_ptr<PointHistory<dim, ADNumberType> > >\r
lqph = quadrature_point_history.get_data(cell);\r
Assert(lqph.size() == n_q_points, ExcInternalError());\r
solution_total);\r
\r
FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>\r
- cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+ cell (IteratorFilters::LocallyOwnedCell(),\r
dof_handler_ref.begin_active()),\r
- endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+ endc (IteratorFilters::LocallyOwnedCell(),\r
dof_handler_ref.end());\r
for (; cell != endc; ++cell)\r
{\r
- Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());\r
+ Assert(cell->is_locally_owned(), ExcInternalError());\r
+ Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());\r
+\r
assemble_system_one_cell(cell, scratch_data, per_task_data);\r
copy_local_to_global_system(per_task_data);\r
}\r
ScratchData_ASM<ADNumberType> &scratch,\r
PerTaskData_ASM &data) const\r
{\r
- Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());\r
+ Assert(cell->is_locally_owned(), ExcInternalError());\r
\r
data.reset();\r
scratch.reset();\r
void Solid<dim>::update_end_timestep()\r
{\r
FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>\r
- cell (IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+ cell (IteratorFilters::LocallyOwnedCell(),\r
dof_handler_ref.begin_active()),\r
- endc (IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+ endc (IteratorFilters::LocallyOwnedCell(),\r
dof_handler_ref.end());\r
for (; cell!=endc; ++cell)\r
{\r
- Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError());\r
+ Assert(cell->is_locally_owned(), ExcInternalError());\r
+ Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());\r
+\r
const std::vector<std::shared_ptr<PointHistory<dim, ADNumberType> > >\r
lqph = quadrature_point_history.get_data(cell);\r
Assert(lqph.size() == n_q_points, ExcInternalError());\r
class FilteredDataOut : public DataOut<dim, DH>\r
{\r
public:\r
- FilteredDataOut (const unsigned int subdomain_id)\r
- :\r
- subdomain_id (subdomain_id)\r
+ FilteredDataOut ()\r
{}\r
\r
virtual ~FilteredDataOut() {}\r
typename DataOut<dim, DH>::active_cell_iterator\r
cell = this->dofs->begin_active();\r
while ((cell != this->dofs->end()) &&\r
- (cell->subdomain_id() != subdomain_id))\r
+ (!cell->is_locally_owned()))\r
++cell;\r
return cell;\r
}\r
{\r
if (old_cell != this->dofs->end())\r
{\r
- const IteratorFilters::SubdomainEqualTo predicate(subdomain_id);\r
+ const IteratorFilters::LocallyOwnedCell predicate{};\r
return\r
++(FilteredIterator<typename DataOut<dim, DH>::active_cell_iterator>\r
(predicate,old_cell));\r
else\r
return old_cell;\r
}\r
-\r
- private:\r
- const unsigned int subdomain_id;\r
};\r
\r
template<int dim, class DH=DoFHandler<dim> >\r
class FilteredDataOutFaces : public DataOutFaces<dim,DH>\r
{\r
public:\r
- FilteredDataOutFaces (const unsigned int subdomain_id)\r
- :\r
- subdomain_id (subdomain_id)\r
+ FilteredDataOutFaces ()\r
{}\r
\r
virtual ~FilteredDataOutFaces() {}\r
{\r
typename DataOutFaces<dim,DH>::active_cell_iterator\r
cell = this->dofs->begin_active();\r
- while ((cell!=this->dofs->end()) && (cell->subdomain_id()!=subdomain_id))\r
+ while ((cell!=this->dofs->end()) && (!cell->is_locally_owned()))\r
++cell;\r
return cell;\r
}\r
{\r
if (old_cell!=this->dofs->end())\r
{\r
- const IteratorFilters::SubdomainEqualTo predicate(subdomain_id);\r
+ const IteratorFilters::LocallyOwnedCell predicate{};\r
return\r
++(FilteredIterator<typename DataOutFaces<dim,DH>::active_cell_iterator>\r
(predicate,old_cell));\r
else\r
return old_cell;\r
}\r
-\r
- private:\r
- const unsigned int subdomain_id;\r
};\r
\r
//Class to compute gradient of the pressure\r
\r
//Iterate through elements (cells) and Gauss Points\r
FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>\r
- cell(IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+ cell(IteratorFilters::LocallyOwnedCell(),\r
dof_handler_ref.begin_active()),\r
- endc(IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+ endc(IteratorFilters::LocallyOwnedCell(),\r
dof_handler_ref.end()),\r
- cell_v(IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+ cell_v(IteratorFilters::LocallyOwnedCell(),\r
vertex_handler_ref.begin_active()),\r
- cell_v_vec(IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+ cell_v_vec(IteratorFilters::LocallyOwnedCell(),\r
vertex_vec_handler_ref.begin_active());\r
//start cell loop\r
for (; cell!=endc; ++cell, ++cell_v, ++cell_v_vec)\r
{\r
- if (cell->subdomain_id() != this_mpi_process) continue;\r
- material_id(cell->active_cell_index())=\r
- static_cast<int>(cell->material_id());\r
+ Assert(cell->is_locally_owned(), ExcInternalError());\r
+ Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());\r
+\r
+ material_id(cell->active_cell_index())=\r
+ static_cast<int>(cell->material_id());\r
\r
fe_values_ref.reinit(cell);\r
\r
}\r
\r
// Add the results to the solution to create the output file for Paraview\r
- FilteredDataOut<dim> data_out(this_mpi_process);\r
+ FilteredDataOut<dim> data_out;\r
std::vector<DataComponentInterpretation::DataComponentInterpretation>\r
comp_type(dim,\r
DataComponentInterpretation::component_is_part_of_vector);\r
\r
//Iterate through elements (cells) and Gauss Points\r
FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>\r
- cell(IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+ cell(IteratorFilters::LocallyOwnedCell(),\r
dof_handler_ref.begin_active()),\r
- endc(IteratorFilters::SubdomainEqualTo(this_mpi_process),\r
+ endc(IteratorFilters::LocallyOwnedCell(),\r
dof_handler_ref.end());\r
//start cell loop\r
for (; cell!=endc; ++cell)\r
{\r
-\r
- if (cell->subdomain_id() != this_mpi_process) continue;\r
+ Assert(cell->is_locally_owned(), ExcInternalError());\r
+ Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError());\r
\r
fe_values_ref.reinit(cell);\r
\r
std::vector<double> displ_incr(dim, 0.0);\r
(void)boundary_id;\r
(void)direction;\r
- throw std::runtime_error ("Displacement loading not implemented for Ehlers verification examples.");\r
+ AssertThrow(false, ExcMessage("Displacement loading not implemented for Ehlers verification examples."));\r
\r
return displ_incr;\r
}\r
std::vector<double> displ_incr(dim, 0.0);\r
(void)boundary_id;\r
(void)direction;\r
- throw std::runtime_error ("Displacement loading not implemented for Ehlers verification examples.");\r
+ AssertThrow(false, ExcMessage("Displacement loading not implemented for Ehlers verification examples."));\r
\r
return displ_incr;\r
}\r
std::vector<double> displ_incr(dim, 0.0);\r
(void)boundary_id;\r
(void)direction;\r
- throw std::runtime_error ("Displacement loading not implemented for Franceschini examples.");\r
+ AssertThrow(false, ExcMessage("Displacement loading not implemented for Franceschini examples."));\r
\r
return displ_incr;\r
}\r
We implemented a nonlinear poro-viscoelastic formulation with the aim of characterising brain tissue response to cyclic loading. Our model captures both experimentally observed fluid flow and conditioning aspects of brain tissue behavior in addition to its well-established nonlinear, preconditioning, hysteretic, and tension-compression asymmetric characteristics.
-The tissue is modelled as a biphasic material consisting in an immiscible aggregate of a nonlinear viscoelastic solid skeleton saturated with pore fluid. The governing equations are linearised using automatic differentiation and solved monolithically for the unknown solid displacements and fluid pore pressure values.
+The tissue is modelled as a biphasic material consisting of an immiscible aggregate of a nonlinear viscoelastic solid skeleton saturated with pore fluid. The governing equations are linearised using automatic differentiation and solved monolithically for the unknown solid displacements and fluid pore pressure values.
-A detailed description of the formulation, it's verification, and the results obtained can be found in:
-* E. Comellas, S. Budday, J.-P. Pelteret, G. A. Holzapfel and P. Steinamnn (2020), Modeling the porous and viscous responses of human brain tissue behavior, Computer Methods in Applied Mechanics and Engineering (accepted for publication).
+A detailed description of the formulation, its verification, and the results obtained can be found in:
+* E. Comellas, S. Budday, J.-P. Pelteret, G. A. Holzapfel and P. Steinmann (2020), Modeling the porous and viscous responses of human brain tissue behavior, Computer Methods in Applied Mechanics and Engineering, 113128. DOI: [10.1016/j.cma.2020.113128](https://doi.org/10.1016/j.cma.2020.113128);
In this paper we show that nonlinear poroelasticity alone can reproduce consolidation experiments, yet it is insufficient to capture stress conditioning due to cyclic loading. We also discuss how the poroelastic response exhibits preconditioning and hysteresis in the fluid flow space, with porous and viscous effects being highly interrelated.
-Quick facts on the code
------------------------
+Quick facts about the code
+--------------------------
* Biphasic material following the Theory of Porous Media
* Nonlinear finite viscoelasticity built on Ogden hyperelasticity
* Darcy-like fluid flow
-* Spatial discretisation with continuous Q2P1 Lagrangian finite elements
+* Spatial discretisation with continuous Q2-P1 Lagrangian finite elements
* Temporal discretisation with a stable implicit one-step backward differentiation method
* Newton-Raphson scheme to solve the nonlinear system of governing equations
-* Forward mode automatic differentiation with the number of derivative components chosen at run-time (Sacado library within Trilinos package) to linearise the equations
+* Forward mode automatic differentiation with the number of derivative components chosen at run-time (Sacado library within Trilinos package) to linearise the governing equations (and, implicitly, the constitutive laws)
* Trilinos direct solver for the (non-symmetric) linear system of equations using a monolithic scheme
-* Parallelization through Threading Building Blocks and across nodes via MPI (Trilinos package)
-* Based on step-44 and code gallery contributions 'Quasi-Static Finite-Strain Compressible Elasticity' and 'Quasi-Static Finite-Strain Quasi-incompressible Visco-elasticity'
+* Parallelization through Threading Building Blocks and across nodes via MPI (using Trilinos linear algebra)
+* Based on step-44 and the code gallery contributions 'Quasi-Static Finite-Strain Compressible Elasticity' and 'Quasi-Static Finite-Strain Quasi-incompressible Visco-elasticity'
* Only works in 3D
Running the code
----------------
### Requirements
-* Version 9.1 or greater of deal.II
-* C++11, MPI, Metis and Trilinos with Sacado must be enabled
+* MPI and Trilinos (built with the Sacado library) must be enabled
### Compiling and running
Similar to the example programs, run
The 'run-multi-calc.py' and 'runPoro.sh' files provided must both be in the main directory. This will automatically generate the required input files and run them in sequence.
+
Reference for this work
-----------------------
-If you use this program as a basis for your own work, please consider citing the paper referenced in the introduction. The initial version of this work was contributed to the deal.II project by E. Comellas and J.-P. Pelteret.
+If you use this program as a basis for your own work, please consider citing the paper referenced in the introduction. The initial version of this work was contributed to the deal.II project by E. Comellas and J-P. Pelteret.
* G.A. Holzapfel (2001), Nonlinear Solid Mechanics. A Continuum Approach for Engineering, John Wiley & Sons. ISBN: [978-0-471-82319-3](http://eu.wiley.com/WileyCDA/WileyTitle/productCd-0471823198.html);
+Results
+----------------------
+The results shown here are a selection of those presented and discussed in the paper referenced in the introduction.
+
+### Consolidation experiments
+We reproduce the uniaxial consolidation experiments of brain tissue from the seminal paper by Franceschini et al. (2006) using a reduction of our formulation to nonlinear poroelasticity without the viscous component.
+
+The following geometry (320 cells and 9544 degrees of freedom), boundary conditions, and loading are used:
+
+
+
+The material properties are:
+* First Lamé parameter: 334kPa
+* Ogden parameters: μ<sub>1</sub>=1.044kPa, μ<sub>2</sub>=1.183kPa, α<sub>1</sub>=4.309, and α<sub>2</sub>=7.736
+* Initial solid volume fraction: 0.80
+* Isotropic initial intrinsic permeability: 8.0e-11mm<sup>2</sup>
+* Deformation-dependency control parameter for the permeability: 40
+* Fluid viscosity: 0.89mPa·s
+
+We consider the effect of gravity in our formulation, with the simplifying assumption that the fluid and solid density are both 0.997mg/mm<sup>3</sup>.
+
+Using these parameters and simulation conditions, we obtain a reasonable fit to the experimental curve:
+
+
+
+What is interesting about these results is that they show that nonlinear poroelasticity alone, i.e. with a purely hyperelastic solid component, can capture the consolidation behavior of brain tissue in response to an [oedometric test](https://en.wikipedia.org/wiki/Oedometer_test).
+
+However, as we will see in the next section, the viscous component is necessary to also capture stress conditioning due to cyclic loading.
+
+
+### Cyclic experiments under multiple loading modes
+Budday et al. (2017) performed an exhaustive set of cyclic loading experiments on brain tissue. They showed in this and subsequent publications that a finite viscoelastic Ogden model is able to predict essential features like nonlinearity, preconditioning, hysteresis, and tension-compression asymmetry. However, such a monophasic viscoelastic formulation can only implicitly capture the effects due to the extracellular fluid flow within the tissue. This was the motivation behind the formulation presented here: to provide a comprehensive model for brain tissue response that will enable us to explore the porous and viscous effects in its mechanical response to a wide range of loading scenarios.
+
+As an example, we show here the results for cyclic compressive loading. The following geometry (512 cells and 15468 degrees of freedom), boundary conditions, and loading are used:
+
+
+
+The material properties are:
+* First Lamé parameter: 24.5kPa
+* Ogden hyperelastic parameters: μ<sub>∞,1</sub>=-83.9Pa, and α<sub>∞,1</sub>=-11.92
+* Ogden viscous parameters: μ<sub>1</sub>=-2100Pa, and α<sub>1</sub>=-2.2
+* Deformation-independent viscosity of the solid: 14kPa·s
+* Initial solid volume fraction: 0.75
+* Isotropic initial intrinsic permeability: 1.0e-8mm<sup>2</sup>
+* Deformation-dependency control parameter for the permeability: 50
+* Fluid viscosity: 0.89mPa·s
+
+To simplify the problem, we neglect the effect of gravity for this example. Here is an animation of our results, visualised with Paraview:
+
+[](https://youtu.be/1a4LbMXrnrs)
+
+To compare with the experimental results, our code computes the nominal stress on the loading surface and provides the values for each timestep in the "data-for-gnuplot.sol" output file. We can then plot compressive nominal stress versus the averaged vertical stretch:
+
+
+
+
+
+We see that viscosity is required in the solid component to reproduce the preconditioning and hysteretic response seen in the experiments.
+
+We were somewhat surprised by these results, because we were expecting to see preconditioning and hysteresis with poroelasticity alone. We discussed this in a conference talk:
+
+* E. Comellas, J.-P. Pelteret, S. Budday and P. Steinmann (2018). Unsolved issues in the numerical modelling of experimentally-observed porous effects in brain tissue behaviour, 6th European Conference on Computational Mechanics and 7th European Conference on Computational Fluid Dynamics (ECCM-ECFD 2018), Glasgow (UK), 11th–15th June 2018. DOI: [10.13140/RG.2.2.18553.29283](https://doi.org/10.13140/RG.2.2.18553.29283);
+
+So, we set out to explore why this is not the case. For that, we studied the porous and viscous dissipation, which are computed in the code and also provided in the "data-for-gnuplot.sol" output file. We determined that there is dissipation occurring in the poroelastic case, but we were barely seeing it in the stress-stretch plot. Even if we played around with the material parameters to increase the amount of dissipation, the slight hysteresis in our plot barely changed. *Where is all that energy going?*
+
+We found the answer by looking closely into the thermodynamic basis of our constitutive equations. The viscous dissipation is a function of the viscous part of the stress tensor and the viscosity of the solid component. However, the porous dissipation depends on the porous material parameters and the seepage velocity, i.e. the relative velocity of the fluid with respect to the deforming solid component, which is proportional to the gradient of the pressure. Now, the plots we were studying are in the displacement-related space. We were looking for the porous dissipation in the wrong place! Our unknowns, the displacements and pressure, are slightly coupled. This is why we see a bit of porous dissipation in the stress-stretch plot for the poroelastic example.
+
+Indeed, when we computed an equivalent to the nominal stress but in the pressure-related space, we found the "missing" dissipation. The "accumulated outgoing fluid volume" in the pressure-space is the equivalent to the nominal stress in the displacement-space. We plot it versus the fluid reaction pressure on the loading surface, which can be loosely seen as an equivalent to the averaged vertical stretch in the previous plots.
+
+
+
+
+### Our takeaways
+In developing this formulation and running the examples provided in this code gallery, we have realised the importance of the loading and boundary conditions selected in the problems. Inhomogeneous displacements and pressure distributions require a fine enough discretization as well as a careful set-up of representative loading and boundary conditions. Without inhomogeneous stretch and pressure maps, the predicted viscous and porous behaviours of the material will not be realistic.
+
+We have started to explore the relation between viscous and porous effects in response to different types of loading. Interestingly, the fluid accumulated outside the sample increases for each loading cycle in the poroelastic model, while it decreases for the poro-viscoelastic one.
+
+
+
+
+The evolution of pressure in the centre of the loading surface, where we observe a pressure concentration due to loading, has a noticeably different pattern for the two models. Note also the delay between the peak pressure values and peak displacement ones.
+
+
+
+
+We are not sure yet how much of these differences is due to the incorporation of the viscous component, and how much can be explained by the choice of material parameters. (Note that the solid part of the poroelastic example has been adjusted with the Franceschini experiment, but the poro-viscoelastic one used the Budday experiments.) We will keep exploring with more numerical examples, but additional experiments are needed to adjust our material parameters in order to better distinguish the porous and viscous parts of brain tissue response to loading.
+
+We hope the community finds this code useful, whether it is for brain mechanics or other applications, and look forward to learning and discussing about new insights gained about the interrelations of the viscous, porous, and elastic parts of the nonlinear poro-viscoelastic formulation.
projects / code-gallery.git / commitdiff
