Section: New Results
Numerical methods for time domain wave propagation
Higher order time discretization of second order hyperbolic problems
Participants : Sébastien Impériale, Patrick Joly, Jerónimo Rodríguez.
Following the construction of optimal (namely maximizing the allowed CFL number) higher order schemes for second order hyperbolic problems of the form
we have extended in the case of first order system of the form
where, here, is an antisymmetric operator. The corresponding results are very similar.
Energy preserving schemes for hamiltonian systems
Participants : Patrick Joly, Juliette Chabassier.
A classical way of building numerical schemes for PDEs is to preserve, on a discrete level, a continuous conserved quantity, like an energy for instance. In 2008 we were interested in the so called “nonlinear hamiltonian systems of wave equations”, which can be written under the form , where the unknown is valued in , and the function U maps to itself. These systems preserve the continuous energy . This property leads, if the function U satisfies a sort of coercivity property, to a stability result for the solution of the PDE, hence the will to keep an analogous quantity conserved step to step in a discrete way.
We developed in 2008 an implicit, energy conservative, unconditionally stable, second order accurate scheme which can be written for any value of m . We tried to see our numerical scheme as a generalization of the classical -schemes for linear PDE, and we introduced a class of “preserving schemes” that fulfill a certain condition leading to a discrete energy preservation. We showed that explicit schemes do not fit this class except for linear models. We also showed that a “partially decoupled scheme”, based on a more intuitive gradient discretization that consists in approximating kU at the time tn with its linearization (directional finite difference of U along its kth variable between times n + 1 and n-1 and explicit along all other variables) does neither fit this class except for trivial models. These remarks point out our numerical scheme as a natural extension of -schemes for the nonlinear systems' case. This work has lead to international communications in 2009 (SANUM'09, WAVES'09). A technical INRIA report and an article submission are in progress.
We also extended the numerical scheme to more general energy preserving nonlinear systems : if the function maps now to , we call 1U its gradient along the first variable and 2U its gradient along the second variable . The new systems are written under the form . An efficient computation has been lead in C++, using a modified Newton's technique for each time step resolution as well as advanced algorithmic optimizations.
Hybrid Meshes for High-Order Discontinuous Galerkin Methods
Participants : Morgane Bergot, Gary Cohen, Marc Duruflé.
In the frame of her thesis directed by G. Cohen and M. Duruflé, M. Bergot studied a way to mix high-order hexahedral elements with tetrahedra, wedges and pyramids for discontinuous Galerkin methods. Hexahedra and tetrahedra have been well studied , and wedge elements are classically obtained by the tensorial product of a triangular element and a 1D element. The main effort is then devoted to the study of pyramidal elements, considering the following elements in the same mesh:
classical hexahedral elements with Gauss-Lobatto points;
Hesthaven tetrahedra, with Gauss-Lobatto points on the edges;
wedges resulting from the tensorial product of a triangle with Hesthaven points, and an edge with Gauss-Lobatto points;
pyramids joining triangular and quadrangular faces in a continuous way.
A new family of high-order nodal pyramidal finite element which can be used in hybrid meshes which include hexahedra, tetrahedra, wedges and pyramids has been proposed at any order, along with optimal approximate integration rules for the evaluation of the finite element matrices, so that the order of convergence of the method is conserved.Numerical results have been made to demonstrate the efficiency of hybrid meshes compared to pure tetrahedral meshes or hexahedral meshes obtained by splitting tetrahedra into hexahedra.
Coupling Retarded Potentials and Discontinuous Galerkin Methods for time dependent wave propagation
Participant : Patrick Joly.
This topic is developed in collaboration with J. Rodriguez (Santiago de Compostela), T. Abboud (IMACS) and I. Terrasse (EADS) in the framework of the contract ADNUMO with AIRBUS. Let us recall that our objective was to use time-domain integral equations (developed in particular at IMACS and EADS) as a tool for contructing transparent boundary conditions for wave problems in unbounded media. More precisely, we wished couple the use of (possibly high order) discontinuous Galerkin methods for the numerical approximation of the interior equations with a space-time Galerkin approximation of the integral equations that represent the transparent boundary conditions, using possibly a local time stepping for the interior equations (see the activity report of 2008 for the motivation).
This year was the last year of the contract ADNUMO. Two methods have been developed during the last two years for the case of the standard 3D wave equation: a conservative method (CM) (that preserves a discrete energy by construction - see the activity report of 2008) And a so-called post-processed method (PPM) that computes directly a solution that results from a post-processing of the solution given by (CM) (this second method was obtained by adapting the ideas developed in the past for volumic time stepping methods) . These two methods have been implemented and tested numerically. The emphasis has been put on designing an optimal implementation (two algorithms have been considered and compared). The numerical results we have obtained are quite good and have been compiled in a contract report. In short, these results confirm that (PPM) provides a better accuracy than (CM) for an equivalent computational cost. In particular, in the case of local time stepping, it reduces spurious phenomena induced by the aliasing process.
The theoretical aspects of our work have also be completed. This work has been widely presented in various conferences and a esearch article is being written.
The subject is now entering a more operational phase that should achieved in the framework of a new contract ADNUMO2. One additional objective is to extend the work to linearized Euler equations for the applications in aeroacoustics.
Numerical MicroLocal Analysis (NMLA)
Participants : Jean-David Benamou, Anne-Sophie Bonnet-Ben Dhia, Francis Collino, Patrick Joly, Simon Marmorat.
Part I: Application to seismic data
Speed analysis consists in the approximation of the smooth part of the underground wave speed using surface data called seismograms. It is one of the components of seismic imaging and and inverse problem : one tries to recover the seismogram using a numerical model that takes the sought speed as input. Asymptotic models (geometrical optics) are a popular choice because of their cheap cost but also as they correspond top the propagation regime for the class of speed considered. In this context one may wish to perform an assymptotic analysis of the seismograms. We try to apply the NMLA technique develloped at INRIA to this problem.
Part II: Application to Non-destructive control
The numerical simulation of acoustic wave propagation yields a scalar field at all points of the grid discretizing the domain. Focusing on the neighborhood of a fixed chosen point, we want to indentify a small number of significative rays (direction, amplitude, arrival time).
This is a difficult ill-posed non linear optimisation problem with many interesting applications. For instance in the coupling of geometric optics and Finite Difference or Finite Elements codes in Electromagnetics. Seismic imaging also need to extract these data for different processing purposes ...
We have obtained promising results on simple synthetic cases (point sources generated fields in slowly varying smooth media) by extending the numerical microlocal analysis technique originally working in the frequency domain to time domain wave propagation.
Efficient numerical solution of the Vlasov-Maxwell system of equations
Participant : Patrick Ciarlet.
A collaboration with Simon Labrunie (Nancy Univ.)
This is a follow-up to a theoretical work initiated with Eric Sonnendrücker and Régine Barthelmé (M3AS, 2007). The aim was to revisit existing approaches (see the works of Heintzé et al during the 90's) to compute the motion of charged particles, using an H(curl, div) conforming discretization method to compute the electromagnetic fields. The improvements over the existing methods are twofold: with the help of the Weighted Regularization Method, one can handle 3D singular geometries, that can generate intense fields; also, one can solve systems in which the charge conservation equation is not satisfied rigorously. The numerical analysis has been carried out for explicit and implicit discrete time schemes. Finally, the so-called correction methods, very popular among engineers and applied physicists, are covered by this study. A paper on this topic has been accepted for publication in M3AS.
Similar problems are currently investigated in 2 D axisymmetric geometries, using the Fourier Singular Complement Method, previously developed for scalar problems.
Numerical Methods for Vlasov-Maxwell Equations
Participants : Gary Cohen, Marc Duruflé, Alexandre Sinding.
There exist a large number of methods for approximating the motion of charged particles. They rely on a suitable discretization of Maxwell's equations, and a suitable discretization of Vlasov's equation.
A. Sinding's thesis, directed by G. Cohen, is devoted to the coupling of higher order hexaedral finite elements for Maxwell's equations with a Particle in Cell (PIC) method for Vlasov's equations. It is part of a joint work between ONERA Toulouse, CEG and INRIA rocquencourt.
Since continuous approximations seem more fitted to Vlasov'equations, an original implementation of continuous spectral finite elements for solving Maxwell's equations has been constructed. It is based on a mixed formulation of the system. In order to avoid the occurence of spurious modes, a dissipation term corresponding to the tangential jump of the discontinuous fields is added to the formulation. This penalization term is a good way to get rid of the parasitic modes, by sending them into the complex plane with a negative imaginary part (the parasitic waves becoming evanescent).
This solver (high order, continuous fields, numerical dissipation) has given good results in 2D and 3D . An important advantage of this penalization is that it eliminates spurious modes without introducing a divergence term in the formulation so that there will be no problem with non-convex geometries (see Costabel and Dauge). In particular we don't need to use correction techniques such as the Weighted Regularization (Costabel and Dauge) or the Singular Complement Method (Ciarlet et al.) which require the knowledge of the geometry of the computational domain.
So far this technique has been coupled with a Particle In Cell method using a charge conservative coupling such as described by Eric Sonnendrücker and Sebastien Jund for axisymmetric domains and comparisons remains to be done with alternative methods concerning charge and current conservation.
Modelisation of coupled non linear piano strings
Participants : Juliette Chabassier, Patrick Joly.
This work is developed in collaboration with Antoine Chaigne (UME, ENSTA).
The linear wave equation does not describe the complexity of the piano strings vibration well enough for physics based sound synthesis, and especially the modelization of a grand piano via numerical simulations. The nonlinear coupling between transversal and longitudinal vibrations has to be taken into account, as does the “geometrically exact” model described by Morse and Ingard. The string's stiffness can moreover be modeled by introducing an angle that represents the deviation of the string's section from the normal vector. This angle is coupled to the transversal vibration, leading finally to a nonlinear hamiltonian system of wave equations entering the second general class of systems of § 6.1.2 .
Measurements have been made on on a stringed soundboard, in order to compare velocities at a certain point of the string with numerical results in similar conditions. Experimental studies on the agraffe have confirmed the assumption that this end of the strings does not significantly move, accordingly with the Dirichlet conditions used in our numerical simulations. For piano strings, the boundary condition at the other end of the strings are not as simple as Dirichlet conditions, because of the bridge vibration. Moreover, most of the notes of a piano are performed with several (one, two or three) strings, between which the bridge is the coupling point. An experimental and modelization work is in progress concerning the bridge, which has been for now considered as a damped linear oscillator but seems to behave as an elastic beam.
Our numerical method is based on a variational formulation of the coupled problem, a high order finite element approximation for the space discretization and energy dissipating schemes, extending the ideas of § 6.1.2 , for the time discretization. The energy method allowed us to easily take into account the coupling with a nonlinear hysteretic felt hammer, as well as the coupling between strings and bridge.
Modeling Piezoelectric Sensors
Participants : Gary Cohen, Sébastien Impériale, Patrick Joly.
Non destructive sensors are built up with piezoelectric materials. These materials provide a small displacement when we apply a electric current and vice versa. To simulate their behaviour it is necessary to take into account the coupling between linear elastodynamics equations and Maxwells' equations. A first step has consisted in developing a simplified model. Because of the very large contrast between elastic waves velocity and electromagnetic waves velocity, at the time scale we are interseted in, the appropriate model for the electric part is an electrostatic model. Moreover we are able to restrict the computation of the electric potential on a relatively small subdomain. These approximate models have been justified theoritically, which enable to control the modelling error. For the numerical computation we couple an explicit scheme for the mechanical part of the model with an implicit scheme for the electric part. The different behaviours of the sensors (emission or reception) induce different boundary conditions so appropriate time discretizations have been studied and implemented. Theses discretizations are energy preserving, which ensure the stability of the resolution. Finally, we use a high order spectral finite element method to solve the elastodynamics equations, whereas the particular device's geometry enable us to use a modal resolution to solve the electrostatic problem.
High order time domain tetrahedral finite elements and GPU implementation
Participants : Alexandre Impériale, Sébastien Impériale.
Graphic processing unit offers a cheap way to do intensive computations. The trade off is that the on-chip memory space is small, this constraint motivates the developpment of particular algorithm to exploit GPU technologie for wave propagation problems. We have developped new high order tetrahedral elements. Whereas they do not provide a fully explicit scheme (no mass lumping) the stability properties of theses elements are optimized and nearly optimal. These elements are well suited for GPU since every computations can be made using small local matrices that fit into the GPU memory. Implementation of order 2 and 3 tetrahedral elements on relatively big problems provide encouraging results.
Evolution problems in locally perturbed infinite periodic media
Participants : Julien Coatléven, Patrick Joly.
This work is part of the PhD of J. Coatléven. The principles and theoretical basis of a numerical method has been developped rigorously for the treatement of linear evolution problems in locally perturbed periodic media, for the moment in 1D and in wave-guide geometries. The idea is to take advantage of the periodicity to reduce the effective computations to small neighbourhood of the perturbation, by constructing transparent boundary conditions. The method is based on a semi-discretization in time of the problem in the whole infinite media, and the transparent boundary conditions constructed are showed to provide an interior solution that is exactly the restriction to the computational domain of the solution of the semi-dicretized problem, thus providing automatically stability of the method. The transparent boundary conditions are constructed through the resolution of semi-discretized unperturbed half-guide problems. We show that the solutions of these problems can be computed using a family of recursively defined functions, for which we use the periodicity to restrict our computations to local-cell problems. So far, we have tested numerically the method on several equations : the heat equation, the linear Schrödinger equation and of course the wave equation. The case of a locally perturbed infinite periodic media in has been theoretically solved, and should be numericaly tested soon.
Ultrasonic Waves Propagation in a Bended Stratified Medium with Local Fold's Undulation Defect
Participants : Gary Cohen, Édouard Demaldent, Sébastien Impériale.
We were dedicated to the study and development of a finite element code that simulates the propagation of ultrasonic waves in a bended stratified elastic medium. This medium is embedded in fluid and can present a local fold's undulation defect.
The originality of this study lies in the behaviour of the elasticity tensor which is not expressed in terms of the canonical basis vectors, but in terms of the tangents to the surface layers (including a rotation between each layer) and the normal to their average surface. When the medium is free from defect, the model is based on the shell model of Kirchhof-Love which restricts the set of admissible geometries. First, the function that maps a reference (flat) medium into a bended one has to be a diffeomorphism and, secondly, its cotangent basis has to be orthogonal (this implies a constant thickness in each layer). We define a defect as a transgression of the second law and introduce it by a local variation of the thickness. Figure 2 shows us some computations with and without defect.
A code has been fully developed (written in Fortran 90, almost 30.000 lines): it is operational in the two dimensional space and its extension to the three dimensional space should be confined to pre and post processing tasks. In this code, the mixed spectral finite element method (MSFEM) has been implemented within an inhomogeneous order of approximation technique. This trick is motivated by the low thickness of the layers (that restrains the increase of the order within the classical homogeneous MSFEM), and avoids us to suffer from the gap of celerity between fluid and solid domains (the fluid-solid coupling is nonconformal).
A suitable absorbing layer method had to be introduce to bound the domain of simulation. Classical (cartesian) perfectly matched layer (PML) methods can leads to instabilities since the wavefront depends on the bending, so we have investigated a new procedure to handle this problem. It consists in choosing absorbing directions that coincide with the directions associated to (cotangent directions). Therefore, the change of variables that leads to the PML system is made with respect to the parametric coordinates. The procedure we use can treat the same range of anisotropic materials than classical PML and stands perfectly matched when the mapping function is locally linear. Another procedure that stands perfectly matched with any mapping function has also been tested but it increases computational time whereas results are not improved.
From a more general point of view, an alternative (not perfectly matched) absorbing method that can treat all elastic materials without corrupting too much the overall accuracy has been studied, and an original 'velocity-strain' formulation of linear elastodynamic equations that simplifies the PML equations system and reduces computational costs has been introduced.