From: Jean-Paul Pelteret Date: Sat, 23 May 2020 19:18:24 +0000 (+0200) Subject: Documentation and code improvements X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F72%2Fhead;p=code-gallery.git Documentation and code improvements --- diff --git a/Nonlinear_PoroViscoelasticity/CMakeLists.txt b/Nonlinear_PoroViscoelasticity/CMakeLists.txt index a0a416d..6c48e74 100644 --- a/Nonlinear_PoroViscoelasticity/CMakeLists.txt +++ b/Nonlinear_PoroViscoelasticity/CMakeLists.txt @@ -26,16 +26,12 @@ ENDIF() # Are all dependencies fullfilled? -IF(NOT DEAL_II_WITH_CXX11 OR - NOT DEAL_II_WITH_MPI OR - NOT DEAL_II_WITH_METIS OR +IF(NOT DEAL_II_WITH_MPI OR NOT DEAL_II_WITH_TRILINOS OR NOT DEAL_II_TRILINOS_WITH_SACADO) MESSAGE(FATAL_ERROR " Error! The deal.II library found at ${DEAL_II_PATH} was not configured with - DEAL_II_WITH_CXX11 = ON DEAL_II_WITH_MPI = ON - DEAL_II_WITH_METIS = ON DEAL_II_WITH_TRILINOS = ON DEAL_II_TRILINOS_WITH_SACADO = ON One or all of these are OFF in your installation but are required for this code gallery example." diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/Franceschini-consolidation.png b/Nonlinear_PoroViscoelasticity/doc/images-results/Franceschini-consolidation.png new file mode 100644 index 0000000..18ce85e Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/Franceschini-consolidation.png differ diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-reac_p-poro-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-reac_p-poro-COMP.png new file mode 100755 index 0000000..0efb2ef Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-reac_p-poro-COMP.png differ diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-COMP.png new file mode 100755 index 0000000..61f472e Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-COMP.png differ diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-visco-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-visco-COMP.png new file mode 100755 index 0000000..4576b68 Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_fluid-time-poro-visco-COMP.png differ diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-COMP.png new file mode 100755 index 0000000..c446723 Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-COMP.png differ diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-visco-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-visco-COMP.png new file mode 100755 index 0000000..aed8542 Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_pres-displ-vs-time-poro-visco-COMP.png differ diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-experiments-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-experiments-COMP.png new file mode 100644 index 0000000..75cceee Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-experiments-COMP.png differ diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-COMP.png new file mode 100755 index 0000000..1868863 Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-COMP.png differ diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-visco-COMP.png b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-visco-COMP.png new file mode 100755 index 0000000..901553c Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/brain_cube_stress-stretch-poro-visco-COMP.png differ diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/consolidation-mesh.png b/Nonlinear_PoroViscoelasticity/doc/images-results/consolidation-mesh.png new file mode 100644 index 0000000..6902ba6 Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/consolidation-mesh.png differ diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/cyclic-loading-mesh.png b/Nonlinear_PoroViscoelasticity/doc/images-results/cyclic-loading-mesh.png new file mode 100644 index 0000000..afae378 Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/cyclic-loading-mesh.png differ diff --git a/Nonlinear_PoroViscoelasticity/doc/images-results/poro-visco-compression.gif b/Nonlinear_PoroViscoelasticity/doc/images-results/poro-visco-compression.gif new file mode 100644 index 0000000..5b0c2cb Binary files /dev/null and b/Nonlinear_PoroViscoelasticity/doc/images-results/poro-visco-compression.gif differ diff --git a/Nonlinear_PoroViscoelasticity/nonlinear-poro-viscoelasticity.cc b/Nonlinear_PoroViscoelasticity/nonlinear-poro-viscoelasticity.cc index 74a431f..c96f2b1 100644 --- a/Nonlinear_PoroViscoelasticity/nonlinear-poro-viscoelasticity.cc +++ b/Nonlinear_PoroViscoelasticity/nonlinear-poro-viscoelasticity.cc @@ -33,10 +33,11 @@ #include #include #include -#include #include +#include + #include #include #include @@ -47,7 +48,6 @@ #include #include #include -#include #include #include #include @@ -448,12 +448,12 @@ namespace NonLinearPoroViscoElasticity gravity_value = prm.get_double("gravity value"); if ( (fluid_type == "Markert") && ((init_intrinsic_perm == 0.0) || (viscosity_FR == 0.0)) ) - throw std::runtime_error ("Markert seepage velocity formulation requires the definition of " - "'initial intrinsic permeability' and 'fluid viscosity' greater than 0.0."); + AssertThrow(false, ExcMessage("Markert seepage velocity formulation requires the definition of " + "'initial intrinsic permeability' and 'fluid viscosity' greater than 0.0.")); if ( (fluid_type == "Ehlers") && ((init_darcy_coef == 0.0) || (weight_FR == 0.0)) ) - throw std::runtime_error ("Ehler seepage velocity formulation requires the definition of " - "'initial Darcy coefficient' and 'fluid weight' greater than 0.0."); + AssertThrow(false, ExcMessage("Ehler seepage velocity formulation requires the definition of " + "'initial Darcy coefficient' and 'fluid weight' greater than 0.0.")); const std::string eigen_solver_type = prm.get("eigen solver"); if (eigen_solver_type == "QL Implicit Shifts") @@ -1009,8 +1009,8 @@ namespace NonLinearPoroViscoElasticity } iteration += 1; if (iteration > 15 ) - throw std::runtime_error ("No convergence in local Newton iteration for the " - "viscoelastic exponential time integration algorithm."); + AssertThrow(false, ExcMessage("No convergence in local Newton iteration for the " + "viscoelastic exponential time integration algorithm.")); } NumberType aux_J_e_1 = 1.0; @@ -1222,8 +1222,9 @@ namespace NonLinearPoroViscoElasticity permeability_term = get_darcy_flow_current(F) / weight_FR; else - throw std::runtime_error ("Material_Darcy_Fluid --> Only Markert " - "and Ehlers formulations have been implemented."); + AssertThrow(false, ExcMessage( + "Material_Darcy_Fluid --> Only Markert " + "and Ehlers formulations have been implemented.")); return ( -1.0 * permeability_term * det_F * (grad_p_fluid - get_body_force_FR_current()) ); @@ -1250,7 +1251,9 @@ namespace NonLinearPoroViscoElasticity seepage_velocity = get_seepage_velocity_current(F,grad_p_fluid); } else - throw std::runtime_error ("Material_Darcy_Fluid --> Only Markert and Ehlers formulations have been implemented."); + AssertThrow(false, ExcMessage( + "Material_Darcy_Fluid --> Only Markert and Ehlers " + "formulations have been implemented.")); dissipation_term = ( invert(permeability_term) * seepage_velocity ) * seepage_velocity; dissipation_term *= 1.0/(det_F*det_F); @@ -1423,8 +1426,8 @@ namespace NonLinearPoroViscoElasticity return body_force; } private: - std_cxx11::shared_ptr< Material_Hyperelastic > solid_material; - std_cxx11::shared_ptr< Material_Darcy_Fluid > fluid_material; + std::shared_ptr< Material_Hyperelastic > solid_material; + std::shared_ptr< Material_Darcy_Fluid > fluid_material; }; // @sect3{Nonlinear poro-viscoelastic solid} @@ -1439,7 +1442,7 @@ namespace NonLinearPoroViscoElasticity void run(); protected: - typedef Sacado::Fad::DFad ADNumberType; + using ADNumberType = Sacado::Fad::DFad; std::ofstream outfile; std::ofstream pointfile; @@ -1525,11 +1528,11 @@ namespace NonLinearPoroViscoElasticity const unsigned int this_mpi_process; mutable ConditionalOStream pcout; - //A collection of the parameters used to describe the problem setup + // A collection of the parameters used to describe the problem setup const Parameters::AllParameters ¶meters; - //Declare an instance of dealii Triangulation class (mesh) - Triangulation triangulation; + // Declare an instance of dealii Triangulation class (mesh) + parallel::shared::Triangulation triangulation; // Keep track of the current time and the time spent evaluating certain functions Time time; @@ -1674,7 +1677,7 @@ namespace NonLinearPoroViscoElasticity this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)), pcout(std::cout, this_mpi_process == 0), parameters(parameters), - triangulation(Triangulation::maximum_smoothing), + triangulation(mpi_communicator,Triangulation::maximum_smoothing), time(parameters.end_time, parameters.delta_t), timerconsole( mpi_communicator, pcout, @@ -1993,10 +1996,6 @@ namespace NonLinearPoroViscoElasticity timerconsole.enter_subsection("Setup system"); timerfile.enter_subsection("Setup system"); - // Partition triangulation - GridTools::partition_triangulation (n_mpi_processes, triangulation); - - //Determine number of components per block std::vector block_component(n_components, u_block); block_component[p_fluid_component] = p_fluid_block; @@ -2174,13 +2173,15 @@ namespace NonLinearPoroViscoElasticity //Setup the initial quadrature point data using the info stored in parameters FilteredIterator::active_cell_iterator> - cell (IteratorFilters::SubdomainEqualTo(this_mpi_process), + cell (IteratorFilters::LocallyOwnedCell(), dof_handler_ref.begin_active()), - endc (IteratorFilters::SubdomainEqualTo(this_mpi_process), + endc (IteratorFilters::LocallyOwnedCell(), dof_handler_ref.end()); for (; cell!=endc; ++cell) { - Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError()); + Assert(cell->is_locally_owned(), ExcInternalError()); + Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError()); + const std::vector > > lqph = quadrature_point_history.get_data(cell); Assert(lqph.size() == n_q_points, ExcInternalError()); @@ -2481,13 +2482,15 @@ namespace NonLinearPoroViscoElasticity solution_total); FilteredIterator::active_cell_iterator> - cell (IteratorFilters::SubdomainEqualTo(this_mpi_process), + cell (IteratorFilters::LocallyOwnedCell(), dof_handler_ref.begin_active()), - endc (IteratorFilters::SubdomainEqualTo(this_mpi_process), + endc (IteratorFilters::LocallyOwnedCell(), dof_handler_ref.end()); for (; cell != endc; ++cell) { - Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError()); + Assert(cell->is_locally_owned(), ExcInternalError()); + Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError()); + assemble_system_one_cell(cell, scratch_data, per_task_data); copy_local_to_global_system(per_task_data); } @@ -2526,7 +2529,7 @@ namespace NonLinearPoroViscoElasticity ScratchData_ASM &scratch, PerTaskData_ASM &data) const { - Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError()); + Assert(cell->is_locally_owned(), ExcInternalError()); data.reset(); scratch.reset(); @@ -2747,13 +2750,15 @@ namespace NonLinearPoroViscoElasticity void Solid::update_end_timestep() { FilteredIterator::active_cell_iterator> - cell (IteratorFilters::SubdomainEqualTo(this_mpi_process), + cell (IteratorFilters::LocallyOwnedCell(), dof_handler_ref.begin_active()), - endc (IteratorFilters::SubdomainEqualTo(this_mpi_process), + endc (IteratorFilters::LocallyOwnedCell(), dof_handler_ref.end()); for (; cell!=endc; ++cell) { - Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError()); + Assert(cell->is_locally_owned(), ExcInternalError()); + Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError()); + const std::vector > > lqph = quadrature_point_history.get_data(cell); Assert(lqph.size() == n_q_points, ExcInternalError()); @@ -2799,9 +2804,7 @@ namespace NonLinearPoroViscoElasticity class FilteredDataOut : public DataOut { public: - FilteredDataOut (const unsigned int subdomain_id) - : - subdomain_id (subdomain_id) + FilteredDataOut () {} virtual ~FilteredDataOut() {} @@ -2812,7 +2815,7 @@ namespace NonLinearPoroViscoElasticity typename DataOut::active_cell_iterator cell = this->dofs->begin_active(); while ((cell != this->dofs->end()) && - (cell->subdomain_id() != subdomain_id)) + (!cell->is_locally_owned())) ++cell; return cell; } @@ -2822,7 +2825,7 @@ namespace NonLinearPoroViscoElasticity { if (old_cell != this->dofs->end()) { - const IteratorFilters::SubdomainEqualTo predicate(subdomain_id); + const IteratorFilters::LocallyOwnedCell predicate{}; return ++(FilteredIterator::active_cell_iterator> (predicate,old_cell)); @@ -2830,18 +2833,13 @@ namespace NonLinearPoroViscoElasticity else return old_cell; } - - private: - const unsigned int subdomain_id; }; template > class FilteredDataOutFaces : public DataOutFaces { public: - FilteredDataOutFaces (const unsigned int subdomain_id) - : - subdomain_id (subdomain_id) + FilteredDataOutFaces () {} virtual ~FilteredDataOutFaces() {} @@ -2851,7 +2849,7 @@ namespace NonLinearPoroViscoElasticity { typename DataOutFaces::active_cell_iterator cell = this->dofs->begin_active(); - while ((cell!=this->dofs->end()) && (cell->subdomain_id()!=subdomain_id)) + while ((cell!=this->dofs->end()) && (!cell->is_locally_owned())) ++cell; return cell; } @@ -2861,7 +2859,7 @@ namespace NonLinearPoroViscoElasticity { if (old_cell!=this->dofs->end()) { - const IteratorFilters::SubdomainEqualTo predicate(subdomain_id); + const IteratorFilters::LocallyOwnedCell predicate{}; return ++(FilteredIterator::active_cell_iterator> (predicate,old_cell)); @@ -2869,9 +2867,6 @@ namespace NonLinearPoroViscoElasticity else return old_cell; } - - private: - const unsigned int subdomain_id; }; //Class to compute gradient of the pressure @@ -3028,20 +3023,22 @@ namespace NonLinearPoroViscoElasticity //Iterate through elements (cells) and Gauss Points FilteredIterator::active_cell_iterator> - cell(IteratorFilters::SubdomainEqualTo(this_mpi_process), + cell(IteratorFilters::LocallyOwnedCell(), dof_handler_ref.begin_active()), - endc(IteratorFilters::SubdomainEqualTo(this_mpi_process), + endc(IteratorFilters::LocallyOwnedCell(), dof_handler_ref.end()), - cell_v(IteratorFilters::SubdomainEqualTo(this_mpi_process), + cell_v(IteratorFilters::LocallyOwnedCell(), vertex_handler_ref.begin_active()), - cell_v_vec(IteratorFilters::SubdomainEqualTo(this_mpi_process), + cell_v_vec(IteratorFilters::LocallyOwnedCell(), vertex_vec_handler_ref.begin_active()); //start cell loop for (; cell!=endc; ++cell, ++cell_v, ++cell_v_vec) { - if (cell->subdomain_id() != this_mpi_process) continue; - material_id(cell->active_cell_index())= - static_cast(cell->material_id()); + Assert(cell->is_locally_owned(), ExcInternalError()); + Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError()); + + material_id(cell->active_cell_index())= + static_cast(cell->material_id()); fe_values_ref.reinit(cell); @@ -3274,7 +3271,7 @@ namespace NonLinearPoroViscoElasticity } // Add the results to the solution to create the output file for Paraview - FilteredDataOut data_out(this_mpi_process); + FilteredDataOut data_out; std::vector comp_type(dim, DataComponentInterpretation::component_is_part_of_vector); @@ -3532,15 +3529,15 @@ namespace NonLinearPoroViscoElasticity //Iterate through elements (cells) and Gauss Points FilteredIterator::active_cell_iterator> - cell(IteratorFilters::SubdomainEqualTo(this_mpi_process), + cell(IteratorFilters::LocallyOwnedCell(), dof_handler_ref.begin_active()), - endc(IteratorFilters::SubdomainEqualTo(this_mpi_process), + endc(IteratorFilters::LocallyOwnedCell(), dof_handler_ref.end()); //start cell loop for (; cell!=endc; ++cell) { - - if (cell->subdomain_id() != this_mpi_process) continue; + Assert(cell->is_locally_owned(), ExcInternalError()); + Assert(cell->subdomain_id() == this_mpi_process, ExcInternalError()); fe_values_ref.reinit(cell); @@ -4092,7 +4089,7 @@ namespace NonLinearPoroViscoElasticity std::vector displ_incr(dim, 0.0); (void)boundary_id; (void)direction; - throw std::runtime_error ("Displacement loading not implemented for Ehlers verification examples."); + AssertThrow(false, ExcMessage("Displacement loading not implemented for Ehlers verification examples.")); return displ_incr; } @@ -4320,7 +4317,7 @@ namespace NonLinearPoroViscoElasticity std::vector displ_incr(dim, 0.0); (void)boundary_id; (void)direction; - throw std::runtime_error ("Displacement loading not implemented for Ehlers verification examples."); + AssertThrow(false, ExcMessage("Displacement loading not implemented for Ehlers verification examples.")); return displ_incr; } @@ -4480,7 +4477,7 @@ namespace NonLinearPoroViscoElasticity std::vector displ_incr(dim, 0.0); (void)boundary_id; (void)direction; - throw std::runtime_error ("Displacement loading not implemented for Franceschini examples."); + AssertThrow(false, ExcMessage("Displacement loading not implemented for Franceschini examples.")); return displ_incr; } diff --git a/Nonlinear_PoroViscoelasticity/readme.md b/Nonlinear_PoroViscoelasticity/readme.md index 4981df8..372c96d 100644 --- a/Nonlinear_PoroViscoelasticity/readme.md +++ b/Nonlinear_PoroViscoelasticity/readme.md @@ -6,35 +6,34 @@ Overview 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 @@ -70,10 +69,11 @@ python run-multi-calc.py 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. @@ -87,3 +87,91 @@ Recommended literature * 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: + +![Consolidation geometry](./doc/images-results/consolidation-mesh.png) + +The material properties are: +* First Lamé parameter: 334kPa +* Ogden parameters: μ1=1.044kPa, μ2=1.183kPa, α1=4.309, and α2=7.736 +* Initial solid volume fraction: 0.80 +* Isotropic initial intrinsic permeability: 8.0e-11mm2 +* 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/mm3. + +Using these parameters and simulation conditions, we obtain a reasonable fit to the experimental curve: + +![Consolidation results](./doc/images-results/Franceschini-consolidation.png) + +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: + +![Cyclic loading geometry](./doc/images-results/cyclic-loading-mesh.png) + +The material properties are: +* First Lamé parameter: 24.5kPa +* Ogden hyperelastic parameters: μ∞,1=-83.9Pa, and α∞,1=-11.92 +* Ogden viscous parameters: μ1=-2100Pa, and α1=-2.2 +* Deformation-independent viscosity of the solid: 14kPa·s +* Initial solid volume fraction: 0.75 +* Isotropic initial intrinsic permeability: 1.0e-8mm2 +* 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: + +[![Cyclic loading video results](./doc/images-results/poro-visco-compression.gif)](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: + +![Cyclic loading experimental results](./doc/images-results/brain_cube_stress-stretch-experiments-COMP.png) +![Cyclic loading poroelastic results](./doc/images-results/brain_cube_stress-stretch-poro-COMP.png) +![Cyclic loading poro-viscoelastic results](./doc/images-results/brain_cube_stress-stretch-poro-visco-COMP.png) + +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. + +![Cyclic loading poroelastic hysteresis](./doc/images-results/brain_cube_fluid-reac_p-poro-COMP.png) + + +### 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. + +![Cyclic loading poroelastic accum fluid](./doc/images-results/brain_cube_fluid-time-poro-COMP.png) +![Cyclic loading poro-viscoelastic accum fluid](./doc/images-results/brain_cube_fluid-time-poro-visco-COMP.png) + +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. + +![Cyclic loading poroelastic accum fluid](./doc/images-results/brain_cube_pres-displ-vs-time-poro-COMP.png) +![Cyclic loading poro-viscoelastic accum fluid](./doc/images-results/brain_cube_pres-displ-vs-time-poro-visco-COMP.png) + +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.