Section: New Results

Numerical Methods

Numerical methods for computing cyclic steady states

Participants : Ustim Khristenko, Patrick Le Tallec [correspondant] .

This work is focused on two techniques for fast computing of the steady cyclic states of evolution problems in non-linear mechanics with space-time periodicity conditions. This kind of problems can be faced in various applications, for instance in the rolling of a tyre with periodic sculptures as well as in a beating heart. Direct solvers for such problems are not very convenient, since they require the inversion of very large matrices. In order to avoid this, a cyclic solution is usually computed as an asymptotic limit of the associated initial value problem with prescribed initial data. However, when the relaxation time is high, convergence to the limit cycle can be very slow. The first technique considered is the Newton-Krylov method, looking for the unknown initial state that provides the space-time periodic solution. This initial state is defined by the space-time periodicity condition, solved with the Newton-Raphson technique. Since the associated Jacobian cannot be expressed explicitly, the method uses one of the matrix-free Krylov iterative solvers. Using information stored while computing the residual to solve the linear system makes its calculation time negligible with respect to the residual calculation time. The second method is the delayed feedback control: an observer-controller type modification of the standard evolution to the limit cycle by introducing a feedback control term, based on the periodicity error. The main result is the optimal form of the control term for a very general class of linear evolution problems, providing the fastest convergence to the cyclic solution. This control has also been adapted and tested for nonlinear problems. The methods discussed have been assessed using academic applications and they have also been implemented into the Michelin industrial code – applied to the rolling tyre model – as well as into the M3DISIM code for the cardiac contraction problem.

Solving isotropic elastodynamics using potentials

Participants : Sébastien Imperiale [correspondant] , Jorge Albella.

This work has the potential to provide an original efficient method for the computations of elastic waves propagation in soft media (such as biological tissues), based on the property that pressure and shear waves decouple in isotropic media. Towards this direction, we considered the numerical solution of 2D elastodynamics isotropic equations using the decomposition of the displacement fields into potentials. This appears as a challenge for finite element methods, and we have addressed here the particular question of free boundary conditions. A stable (mixed) variational formulation of the evolution problem is proposed.

The Arlequin method for transient wave scattering by small obstacles

Participants : Sébastien Imperiale [correspondant] , Jorge Albella.

In this work we extend the Arlequin method to overlapping domain decomposition technique for transient wave equation scattering by small obstacles. The main contribution of this work is to construct and analyze some variants of the Arlequin method from the continuous level to the fully discrete level. The constructed discretizations allow to solve wave propagation problems while using non-conforming and overlapping meshes for the background propagating medium and the surrounding of the obstacle, respectively. Hence we obtain a flexible and stable method in terms of the space discretization – an inf-sup condition is proven – while the stability of the time discretization is ensured by energy identities.

Construction of a fourth-order time scheme for dissipative wave equations

Participants : Sébastien Imperiale [correspondant] , Juliette Chabassier [Magique-3d] , Julien Diaz [Magique-3d] .

This works deals with the construction of a fourth-order, energy preserving, explicit time discretization for dissipative linear wave equations. This discretization is obtained by replacing the inversion of a matrix – that comes naturally after using the technique of the Modified Equation on the second order Leap Frog scheme applied to dissipative linear wave equations – by an explicit approximation of its inverse. The stability of the scheme is studied first using an energy analysis, then an eigenvalue analysis. Numerical results in 1D illustrate the good behavior regarding space/time convergence and the efficiency of the newly derived scheme compared to more classical time discretizations. A loss of accuracy is observed for non-smooth profiles of dissipation, and we propose an extension of the method that fixes this issue. Finally, we assess the good performance of the scheme for a realistic dissipation phenomenon in Lorentz materials.

Coupled variational formulations of linear elasticity and the DPG methodology

Participant : Patrick Le Tallec [correspondant] .

In this work, we develop a general approach akin to domain-decomposition methods to solve a single linear PDE, but where each subdomain of a partitioned domain is associated with a distinct variational formulation coming from a mutually well-posed family of so-called broken variational formulations of the original PDE. It can be exploited to solve challenging problems in a variety of physical scenarios where stability or a particular mode of convergence is desired in some part of the domain. The linear elasticity equations are solved in this work, but the approach can be applied to other equations as well. The broken variational formulations, which are essentially extensions of more standard formulations, are characterized by the presence of mesh-dependent broken test spaces and interface trial variables at the boundaries of the elements of the mesh. This allows necessary information to be naturally transmitted between adjacent subdomains, resulting in coupled variational formulations which are then proved to be globally well-posed. They are solved numerically using the DPG methodology, which is especially crafted to produce stable discretizations of broken formulations. Finally, expected convergence rates are verified in two different illustrative examples. This work has resulted in the publication [19].

A discontinuous Galerkin approach for cardiac electrophysiology

Participant : Radomir Chabiniok [correspondant] .

Cardiac electrophysiology simulations are numerically challenging due to the propagation of a steep electrochemical wave front, and thus require discretizations with small mesh sizes to obtain accurate results. In this work – in collaboration with the Institute for Computational Mechanics, Technical University Munich and published in [21] – we present an approach based on the Hybridizable Discontinuous Galerkin method (HDG), which allows an efficient implementation of high-order discretizations into a computational framework. In particular using the advantage of the discontinuous function space, we present an efficient p-adaptive strategy for accurately tracking the wave front. HDG allows to reduce the overall degrees of freedom in the final linear system to those only on the element interfaces. Additionally, we propose a rule for a suitable integration accuracy for the ionic current term depending on the polynomial order and the cell model to handle high-order polynomials. Our results show that for the same number of degrees of freedom coarse high-order elements provide more accurate results than fine low-order elements. Introducing p-adaptivity further reduces computational costs while maintaining accuracy by restricting the use of high-order elements to resolve the wave front. For a patient-specific simulation of a cardiac cycle, p-adaptivity reduces the average number of degrees of freedom by 95% compared to the non-adaptive model. In addition to reducing computational costs, using coarse meshes with our p-adaptive high-order HDG method also simplifies practical aspects of mesh generation and postprocessing.