Which paper dose the Chapelle 2018 prescribe?

[nabw]

Hi, it is the one found here: https://doi.org/10.1007/s10255-019-0799-5

[nabw]

This code has been accumulating rust for a while. What are you interested in doing? Perhaps there is a better way.

Let me know.

[shaoyaoqian]

We aim to simulate the fluid-structure interaction of the human ventricle, modeling the ventricular tissue using poroelasticity. Our approach is based on the numerical method presented in the paper titled 'A poroelastic immersed finite element framework for modeling cardiac perfusion and fluid–structure interaction' (link: here). (It should be noted that there are many errors in this paper. )

We have worked on it for months, but didn't get the correct numerical results. We tried an example in paper 'Effective and energy-preserving time discretization for a general nonlinear poromechanical formulation'.

The $m$ calculated by our program is:

https://github.com/nabw/poroelasticity-linear-solvers/assets/115222128/c19c7bc6-6c19-4a75-8b1b-29b66d1c7f5f

[shaoyaoqian]

Here is a shot, if you can not open the mp4 video. It is obviously wrong.

[shaoyaoqian]

By the way, we implemented the immersed boundary method with FEniCS as the backend finite element solver.

[nabw]

That quite ambitious! It is not clear by what you say what the error could be. In addition, the numerical approximation of such nonlinear models is still quite open, so many things can go wrong.

I have some doubts regarding their consitutive choices, but the numerical method seems quite transparent. What is it that is not working in your formulation?

[nabw]

We have worked a bit on robust numerical methods for poroelasticity, both linear and nonlinear:

https://doi.org/10.1016/j.cma.2021.114183 https://doi.org/10.4995/YIC2021.2021.12324

Maybe that helps.

[shaoyaoqian]

I noticed that you have authored another paper, accessible at http://arxiv.org/abs/2111.04206. Additionally, I came across a corresponding GitHub repository based on the paper, located at https://github.com/ruizbaier/PoroelasticModelForAcuteMyocarditis.

While reviewing the examples in the paper, I encountered some issues:

1. In the first demonstration, I observed that the boundary condition for pressure at the top of the box is specified as 0.2 MPa. Numerically, it should be expressed as 200,000 Pa, using the same unit as other parameters. We encountered difficulties running the program with this boundary condition. Besides, the pressure appears to close to zero in the figure of the paper.

   

2. Regarding the second example, our program yields an absolute pressure value approximately 25% less than the result presented in the paper.

   

   

I will attach our programs.

[shaoyaoqian]

example 2.1

from fenics import *

parameters["form_compiler"]["representation"] = "uflacs"
parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["quadrature_degree"] = 2

# parameter-dependent functions

normalize_vec = lambda u : u/sqrt(dot(u,u))
kappa  = lambda J,phi: kappa0*(1 + (1-phi0)**2/phi0**3*J*phi**3*(J-phi)**2)  # isotropic Kozeny-Carman

# time constants
t = 0.0
dt = 0.1
Tfinal = 1
cont = 0

# ********* model constants ******* #

E = Constant(1.e4)
nu = Constant(0.33)
ell = Constant(0.0)

rhos = Constant(2.e-3)
rhof = Constant(1.e-3)
kappa0 = Constant(1.e-5)
alpha = Constant(0.25)
muf = Constant(1.e-3)
bb = Constant((0, 0))
I = Identity(2)

# neo-Hookean
mu_s = Constant(E/(2.*(1. + nu)))
lmbdas = Constant(E*nu/((1. + nu)*(1. - 2.*nu)))

# ********* Mesh and I/O ********* #
mesh = UnitSquareMesh(100, 100)
bdry = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)

class Boundary_1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

class Boundary_2(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1.0)

class Boundary_3(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0) 

class Boundary_4(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 1.0)

boundary_1 = Boundary_1()
boundary_2 = Boundary_2()
boundary_3 = Boundary_3()
boundary_4 = Boundary_4()

boundary_1.mark(bdry, 1)
boundary_2.mark(bdry, 2)
boundary_3.mark(bdry, 3)
boundary_4.mark(bdry, 4)

fileO = File("outputs/boundary.pvd")
fileO << bdry

ds = Measure("ds", subdomain_data= bdry)  
nn = FacetNormal(mesh)

fileO_p = File("outputs/drainage/p.pvd")
fileO_u = File("outputs/drainage/u.pvd")
fileO_phi = File("outputs/drainage/phi.pvd")

# ********* Finite dimensional spaces ********* #

P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
Bub = FiniteElement("Bubble", mesh.ufl_cell(), 4)
P1b = VectorElement(P1 + Bub)
P2vec = VectorElement("CG", mesh.ufl_cell(), 1)
Hh = FunctionSpace(mesh, MixedElement([P2vec, P1, P1]))

print("**************** Total Dofs = ", Hh.dim())

Sol = Function(Hh)
dSol = TrialFunction(Hh)
u, phi, p = split(Sol)
v, psi, q = TestFunctions(Hh)

# ******* Boundary conditions ********** #
p_in = Constant(200)
t_1 = Expression('2000*sin(pi*t)', t = t, degree=0)
u_D = project(Constant((0, 0)), Hh.sub(0).collapse())
bcU = DirichletBC(Hh.sub(0), u_D, bdry, 3)
bcP1 = DirichletBC(Hh.sub(2), p_in, bdry, 3)
bcP2 = DirichletBC(Hh.sub(2), 0, bdry, 1)
bcP3 = DirichletBC(Hh.sub(2), 0, bdry, 2)
bcH = [bcU, bcP1, bcP2, bcP3]

# ******* Initial conditions ********** #

phi0 = Expression('0.6 + 0.1*sin(x[0]*x[1])', degree=2)
uold = project(Constant((0, 0)), Hh.sub(0).collapse())
phiold = project(phi0, Hh.sub(1).collapse())
pold = interpolate(Constant(0), Hh.sub(2).collapse())

# ********  Weak form ********** #

ndim = u.geometric_dimension()

F = I + grad(u)
J = det(F)
invF = inv(F)

Peff = mu_s*(F - invF.T) + lmbdas*ln(J)*invF.T

# Poromechanics
F1 = inner(Peff - alpha*p*J*invF.T, grad(v)) * dx \
    + psi*(J-1-phi+phi0)*dx \
    + rhof*J/dt*(phi-phiold)*q * dx \
    + rhof*dot(phi*J*invF*kappa(J, phi)/muf*invF.T*grad(p), grad(q)) * dx \
    - rhos*dot(bb, v) * dx \
    - dot(-J*invF.T*t_1*nn, v) * ds(4) \
    - rhof*J*ell*q*dx

Tang = derivative(F1, Sol, dSol)

# ********* Time loop ************* #

while (t < Tfinal):

    t += dt
    t_1.t = t

    print("t=%.2f" % t)

    # ********* Solving ************* #

    solve(F1 == 0, Sol, J=Tang, bcs=bcH,
          solver_parameters={'newton_solver': {'linear_solver': 'mumps',
                                               'absolute_tolerance': 1.0e-7,
                                               'relative_tolerance': 1.0e-7}})
    uh, phih, ph = Sol.split()

    # ********* Writing ************* #

    fileO_p << ph
    fileO_u << uh
    fileO_phi << phih

    cont += 1

    # ********* Updating ************* #
    assign(uold, uh)
    assign(phiold, phih)
    assign(pold, ph)

# ************* End **************** #

example 2.2

from fenics import *

parameters["form_compiler"]["representation"] = "uflacs"
parameters["allow_extrapolation"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["quadrature_degree"] = 2

# parameter-dependent functions

normalize_vec = lambda u : u/sqrt(dot(u,u))
kappa  = lambda J,phi: kappa0*(1 + (1-phi0)**2/phi0**3*J*phi**3*(J-phi)**2)  # isotropic Kozeny-Carman

# time constants
t = Expression("t", t=0,degree=1)
dt = 1.0
Tfinal = 100
cont = 0

# ********* model constants ******* #

# E = Constant(3.e4)
E = Constant(6000)
nu = Constant(0.2)
ell = Constant(0.0)

rhos = Constant(2.e-3)
rhof = Constant(1.e-3)
kappa0 = Constant(1.e-8)
alpha = Constant(1)
muf = Constant(1.e-3)
bb = Constant((0, 0))
I = Identity(2)

# neo-Hookean
mu_s = Constant(E/(2.*(1. + nu)))
lmbdas = Constant(E*nu/((1. + nu)*(1. - 2.*nu)))

# ********* Mesh and I/O ********* #
mesh = RectangleMesh(Point(0.0, 0.0), Point(8.0, 5.0), 100, 100)
bdry = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)

class Boundary_1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

class Boundary_2(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 8.0)

class Boundary_3(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0) 

class Boundary_4(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 5.0)

class Boundary_5(SubDomain):
    def inside(self, x, on_boundary):
       return on_boundary and near(x[1], 5.0) and (0.0 <= x[0] <= 1.0)  

boundary_1 = Boundary_1()
boundary_2 = Boundary_2()
boundary_3 = Boundary_3()
boundary_4 = Boundary_4()
boundary_5 = Boundary_5()

boundary_1.mark(bdry, 1)
boundary_2.mark(bdry, 2)
boundary_3.mark(bdry, 3)
boundary_4.mark(bdry, 4)
boundary_5.mark(bdry, 5)

fileO = File("outputs/boundary.pvd")
fileO << bdry

ds = Measure("ds", subdomain_data= bdry)  
nn = FacetNormal(mesh)

fileO_p = File("outputs/drainage2/p.pvd")
fileO_u = File("outputs/drainage2/u.pvd")
fileO_phi = File("outputs/drainage2/phi.pvd")

# ********* Finite dimensional spaces ********* #

P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
Bub = FiniteElement("Bubble", mesh.ufl_cell(), 4)
P1b = VectorElement(P1 + Bub)
P2vec = VectorElement("CG", mesh.ufl_cell(), 1)
Hh = FunctionSpace(mesh, MixedElement([P2vec, P1, P1]))

print("**************** Total Dofs = ", Hh.dim())

Sol = Function(Hh)
dSol = TrialFunction(Hh)
u, phi, p = split(Sol)
v, psi, q = TestFunctions(Hh)

# ******* Boundary conditions ********** #
u_D = project(Constant((0, 0)), Hh.sub(0).collapse())
bcU = DirichletBC(Hh.sub(0), u_D, bdry, 3)
bcU1 = DirichletBC(Hh.sub(0).sub(0), Constant(0), bdry, 1)
bcU2 = DirichletBC(Hh.sub(0).sub(0), Constant(0), bdry, 2)
bcP1 = DirichletBC(Hh.sub(2), 0, bdry, 4)
bcP2 = DirichletBC(Hh.sub(2), 0, bdry, 5)
# bcP2 = DirichletBC(Hh.sub(2), 0, bdry, 1)
# bcP3 = DirichletBC(Hh.sub(2), 0, bdry, 2)
# bcP4 = DirichletBC(Hh.sub(2), 0, bdry, 3)
bcH = [bcU, bcP1,bcP2, bcU1, bcU2]

# ******* Initial conditions ********** #
phi0 = Constant(0.33)
uold = project(Constant((0, 0)), Hh.sub(0).collapse())
phiold = project(phi0, Hh.sub(1).collapse())
pold = interpolate(Constant(0), Hh.sub(2).collapse())

# ********  Weak form ********** #

ndim = u.geometric_dimension()

F = I + grad(u)
J = det(F)
invF = inv(F)

Peff = mu_s*(F - invF.T) + lmbdas*ln(J)*invF.T
t_1 = 0.2*(lmbdas+2.0*mu_s)*sin(2*pi*t/15.0)

# Poromechanics
F1 = inner(Peff - alpha*p*J*invF.T, grad(v)) * dx \
    + psi*(J-1-phi+phi0)*dx \
    + rhof*J/dt*(phi-phiold)*q * dx \
    + rhof*dot(phi*J*invF*kappa(J, phi)/muf*invF.T*grad(p), grad(q)) * dx \
    - rhos*dot(bb, v) * dx \
    - dot(-J*t_1*invF.T*nn, v) * ds(5) \
    - rhof*J*ell*q*dx

Tang = derivative(F1, Sol, dSol)

# ********* Time loop ************* #

while (t.t < Tfinal):

    t.t += dt

    print("t=%.2f" % t.t)

    # ********* Solving ************* #

    solve(F1 == 0, Sol, J=Tang, bcs=bcH,
          solver_parameters={'newton_solver': {'linear_solver': 'mumps',
                                               'absolute_tolerance': 1.0e-7,
                                               'relative_tolerance': 1.0e-7}})
    uh, phih, ph = Sol.split()

    # ********* Writing ************* #

    fileO_p << ph
    fileO_u << uh
    fileO_phi << phih

    cont += 1

    # ********* Updating ************* #
    assign(uold, uh)
    assign(phiold, phih)
    assign(pold, ph)

# ************* End **************** #

[nabw]

Hi, I'm very happy that you are looking at our work!

Did you try running exactly the code at the repo? This should be a first test, as I see for example that you are using a very low quadrature degree for assembly (2 instead of 3).

You also mention certain difficulties when running the code. What do you mean exactly by this?

[shaoyaoqian]

Hi, I am back.

I think some parameters and boundary conditions detailed in the paper didn't match the pictures.

When going through the paper more carefully, we are unable to formulate the second equation.

Here is the equation we formulated, without $\frac{1}{J}$ and $\phi_f$

We formulated it from

$$ \frac{\partial}{\partial t}\left(\rho^f \phi\right)+\nabla \cdot\left(\rho^f \phi \underline{u}^f\right)=\rho^f s $$

$$ \underline{w}=-\underline{\underline{K}} \cdot \nabla p $$

$$ \underline{w}=\phi\left(\underline{u}^f-\underline{u}^s\right), $$

This equation is from a paper of Chapelle, named A poroelastic model valid in large strains with applications to perfusion in cardiac modeling.

Could you tell me how to formulate the equation in your paper.

[nabw]

The $\phi$ in the permeability depends on the model used, so (1) it is well justified and (2) it should not make a big difference in the fluid response beyond slightly decreasing the permeability. Additionally, the 1/J comes from the pullback of the forms, which depends on 1) the definition of $\phi$, 2) the definition of the source terms and 3) the type of time derivative under consideration.

Note that the phi in Chapelle's work does not mean the same as it does in ours (Lagrangian vs Eulerian).

It is difficult to help you without more precise questions regarding "non-matching pictures".

Best!

[nabw]

Btw this model is based on MacMinn et al "Large Deformations of a Soft Porous Material", not on Chapelle's formulation, which they later improved in https://doi.org/10.1016/j.euromechflu.2014.02.009 .