Two aspects of solving the incompressible Navier-Stokes equations are described in the thesis. The preconditioning of the algebraic systems arising from the Finite Element Method discretization of the Navier-Stokes equations is complex due to the saddle point structure of the resulting algebraic problems. The Pressure Convection Diffusion Reaction and the Least Squares Commutator preconditioners constitute two possible choices studied in the thesis. Solving the flow problems in time-dependent domains requires special numerical methods, such as the Fictitious Boundary method and the Arbitrary Lagrangian Eulerian formulation of Navier-Stokes equations which are used in the thesis. The problems examined in the thesis are simulations of experiments conducted in liquid Helium at low temperatures. These simulations can be used to establish a relationship between vorticity and new quantity pseudovorticity in an experiment-like setting.

O presente trabalho consiste no uso do Método dos Elementos Finitos (MEF) para a análise de interação fluido-estruturas de barras com foco em problemas transientes envolvendo explosões ou outras ações com propagação de ondas de choque. Para isso é necessário o estudo de três diferentes aspectos: a dinâmica das estruturas computacional, a dinâmica dos fluidos computacional e o problema do acoplamento. No caso da dinâmica das estruturas computacional deve-se identificar em função da cinemática de deformações, quais são os requisitos para que um elemento seja adequado para analisar tais problemas, tendo em vista que a formulação deve admitir grandes deslocamentos. Para evitar problemas relacionados com aproximações de rotações finitas, opta-se por empregar uma formulação descrita em termos de posições e que leva em consideração os efeitos de empenamento da seção transversal. No caso da dinâmica dos fluidos computacional, busca-se uma formulação para escoamentos compressíveis que seja estável e ao mesmo tempo sensível ao movimento da estrutura, sendo empregado um algoritmo de integração temporal explícito baseado em características com as equações governantes descritas na forma Lagrangeana-Euleriana Arbitrária (ALE). No que se refere ao acoplamento, busca-se modularidade e versatilidade, empregando-se um modelo particionado fraco (explícito) de acoplamento e técnicas de transferência das condições de contorno (Dirichlet-Neummann), sendo estudados os efeitos de utilizar transferência bidirecional ou unidirecional dessas condições de contorno. / This work consists in the use of the Finite Element Method (FEM) for numerical analysis of fluid-bar structures, focusing on transient problems involving explosions or other actions with shock waves propagation. For this purpose, one needs to study three different aspects: the computational structural dynamics, the computational fluid dynamics and the coupling problem. Regarding computational structural dynamics, one need firstly to identify the requirements for an element to be adequate to analyze such problems, taking into account the fact that such element should admit large displacements. In order to avoid problems related to finite rotation approximations and to give a realist representation of a 3D bar structure, we chose a formulation defined in terms of positions and that considers the cross-section warping effects. Regarding computational fluid dynamics, we seek for a stable formulation for compressible flows, and at same time, sensitive to the movement of the structure, leading to an explicit time integration algorithm based on characteristics with governing equations described in the Arbitrary Lagrangian-Eulerian (ALE) form. Regarding to coupling, we chose to use a weak (explicit) partitioning coupling model in order to ensure modularity and versatility. The developed coupling scheme is bases on boundary conditions transfer techniques (Dirichlet-Neummann), and we study the effects of using bidirectional or unidirectional boundary conditions transfers.

Die künstliche Grundwasserabsenkung stellt eine wichtige Maßnahme für die Entwässerung von Baugruben und bergbaulich genutzten Flächen dar. Eine erfolgreiche und zielgerichtete Absenkung des Grundwasserspiegels setzt ein zweckmäßiges Design und die richtige Auswahl der genutzten Absenkungstechniken voraus. Dabei sind insbesondere die Dimension des abzusenkenden Bereichs, die Untergrundbeschaffenheit sowie zu erfüllende umweltschutzrechtliche Regelungen zu berücksichtigen. Zur Grundwasserabsenkung kommen üblicherweise verschiedene Ausführungen und Anordnungen von Pumpbrunnen zum Einsatz. Konventionelle Pumpbrunnen, welche für Absenkungsmaßnahmen eingesetzt werden, entnehmen Grundwasser aus dem Aquifer. Durch das fortwährende Abpumpen von in der Regel erheblichen Wassermengen können jedoch Umweltprobleme entstehen, und es ist mit zusätzlichen Entsorgungskosten für die Ableitung des geförderten Wassers zu rechnen. Im Gegensatz hierzu stellen vertikale Zirkulationsbrunnen (VCW) einen innovativen Ansatz dar, der eine lokale Grundwasserabsenkung ohne Nettowasserentnahme aus dem Aquifer erlaubt. Ein VCW kann als ein Einbohrlochsystem aufgefasst werden, bei dem im oberen Bereich eines Brunnens Wasser entnommen und dieses in einem separaten, weiter unten installierten Brunnenbereich wieder injiziert wird. Die erfolgreiche Anwendung dieser neuen Grundwasserabsenkungstechnik erfordert die genaue Kenntnis der Faktoren, welche für die Grundwasserströmungsverhältnisse relevant sind und somit die Absenkung bestimmen. Traditionelle Berechnungsansätze vernachlässigen oft vertikale Grundwasserbewegungen und sind deshalb für die Beschreibung der komplexen Strömungsverhältnisse in unmittelbarer Nähe eines VCW nicht geeignet. Aus diesem Grund steht die systematische Untersuchung der Grundwasserströmung unter Berücksichtigung vertikaler Strömungskomponenten im Hauptfokus dieser Arbeit. Die Untersuchungen beschäftigen sich in erster Linie mit der Entwicklung einer geeigneten Simulationsmethode, mit der Evaluierung des Einflusses relevanter hydrogeologischer Parameter sowie mit der Durchführung und Auswertung von Pumpversuchen an einem Feldstandort. Die hier vorgestellte neue Simulationsmethode koppelt den sogenannten Arbitrary‐Lagrangian‐Eulerian‐(ALE)‐Algorithmus mit der Grundwasserströmungsgleichung. Die Simulationsergebnisse werden mit mehreren analytischen Lösungen verglichen und verifiziert. Das entwickelte numerische Modell berücksichtigt auch Vertikalströmungen und eignet sich somit zur Simulation der Grundwasserströmung in der Nähe von VCW. Folglich kann nun die Lage des Grundwasserspiegels, vor allem für ungespannte Grundwasserleiter, präzise berechnet werden. Nach erfolgter Kalibrierung des numerischen Modells anhand von Felddaten wurde eine Sensitivitätsanalyse relevanter Parameter im Hinblick auf die erzielte Absenkung und deren Auswirkungen auf die Grundwasserströmungssituation durchgeführt. Die dabei erhaltenen Ergebnisse zeigen, dass die Grundwasserabsenkung proportional zur Pumprate, indirekt proportional zur hydraulischen Leitfähigkeit und fast unabhängig von der Anisotropie des Grundwasserleiters um den VCW ist. Des Weiteren zeigte sich, dass die Lage des oberen Entnahmepunktes einen größeren Einfluss auf die Absenkung als die Lage des darunter liegenden Injektionspunktes hat. Die Größe des von der Grundwasserzirkulation beeinflussten Bereiches hängt dagegen neben dem Abstand dieser beiden Punkte hauptsächlich auch von der Anisotropie des Aquifermaterials ab. Um den Einfluss der Hydrostratigraphie auf die Grundwasserströmung zu untersuchen, wurden die Eigenschaften der einzelnen Schichten genau charakterisiert. Hierfür wurden Direct‐Push‐, Pump‐, Injektions‐ sowie Zirkulationsversuche an einem Feldstandort durchgeführt. Zudem wurden Bohrkerne entnommen und mithilfe von Siebanalysen vertikale Korngrößenverteilungsprofile im Labor bestimmt. Die eingesetzten experimentellen Methoden stellen in Kombination mit numerischen Simulationsrechnungen eine gute Basis dar, um die Rolle der Schichtstruktur im Aquifer besser beurteilen zu können. Die Untersuchungen leisten somit einen wichtigen Beitrag für das zukünftige Design und den Betrieb von VCW für Grundwasserabsenkungszwecke in ungespannten Grundwasserleitern. Zudem zeigt die hier vorliegende Arbeit das große Potential dieser neuen Grundwasserabsenkungstechnik als vielversprechende Alternative zu konventionellen Absenkungsverfahren auf.

Title: Fluid-structure interaction of compressible flow Author: RNDr. Jaroslava Hasnedlová Department: Department of Numerical Mathematics, Institute of Applied Mathematics Supervisors: Prof. RNDr. Miloslav Feistauer, DrSc., Dr. h. c., Prof. Dr. Dr. h. c. Rolf Rannacher Supervisors' e-mail addresses: feist@karlin.mff.cuni.cz, rannacher@iwr.uni-heidelberg.de Abstract: The presented work is split into two parts. The first part is devoted to the theory of the discontinuous Galerkin finite element (DGFE) method for the space-time discretization of a nonstationary convection-diffusion initial-boundary value problem with nonlinear convection and linear diffusion. The DGFE method is applied sep- arately in space and time using, in general, different space grids on different time levels and different polynomial degrees p and q in space and time discretization. The main result is the proof of error estimates in L2 (L2 )-norm and in DG-norm formed by the L2 (H1 )-seminorm and penalty terms. The second part of the thesis deals with the realization of fluid-structure interaction problem of the compressible viscous flow with the elastic structure. The time-dependence of the domain occupied by the fluid is treated by the ALE (Arbitrary Lagrangian-Eulerian) method, when the compress- ible Navier-Stokes equations are formulated in...

O presente trabalho consiste no uso do Método dos Elementos Finitos (MEF) para a análise de interação fluido-estruturas de barras com foco em problemas transientes envolvendo explosões ou outras ações com propagação de ondas de choque. Para isso é necessário o estudo de três diferentes aspectos: a dinâmica das estruturas computacional, a dinâmica dos fluidos computacional e o problema do acoplamento. No caso da dinâmica das estruturas computacional deve-se identificar em função da cinemática de deformações, quais são os requisitos para que um elemento seja adequado para analisar tais problemas, tendo em vista que a formulação deve admitir grandes deslocamentos. Para evitar problemas relacionados com aproximações de rotações finitas, opta-se por empregar uma formulação descrita em termos de posições e que leva em consideração os efeitos de empenamento da seção transversal. No caso da dinâmica dos fluidos computacional, busca-se uma formulação para escoamentos compressíveis que seja estável e ao mesmo tempo sensível ao movimento da estrutura, sendo empregado um algoritmo de integração temporal explícito baseado em características com as equações governantes descritas na forma Lagrangeana-Euleriana Arbitrária (ALE). No que se refere ao acoplamento, busca-se modularidade e versatilidade, empregando-se um modelo particionado fraco (explícito) de acoplamento e técnicas de transferência das condições de contorno (Dirichlet-Neummann), sendo estudados os efeitos de utilizar transferência bidirecional ou unidirecional dessas condições de contorno. / This work consists in the use of the Finite Element Method (FEM) for numerical analysis of fluid-bar structures, focusing on transient problems involving explosions or other actions with shock waves propagation. For this purpose, one needs to study three different aspects: the computational structural dynamics, the computational fluid dynamics and the coupling problem. Regarding computational structural dynamics, one need firstly to identify the requirements for an element to be adequate to analyze such problems, taking into account the fact that such element should admit large displacements. In order to avoid problems related to finite rotation approximations and to give a realist representation of a 3D bar structure, we chose a formulation defined in terms of positions and that considers the cross-section warping effects. Regarding computational fluid dynamics, we seek for a stable formulation for compressible flows, and at same time, sensitive to the movement of the structure, leading to an explicit time integration algorithm based on characteristics with governing equations described in the Arbitrary Lagrangian-Eulerian (ALE) form. Regarding to coupling, we chose to use a weak (explicit) partitioning coupling model in order to ensure modularity and versatility. The developed coupling scheme is bases on boundary conditions transfer techniques (Dirichlet-Neummann), and we study the effects of using bidirectional or unidirectional boundary conditions transfers.

Problems governed by partial differential equations (PDEs) in deformable domains, t Rd; d = 2; 3; are of fundamental importance in science and engineering. They are of particular relevance in the design of many engineering systems e.g., aircrafts and bridges as well as to the analysis of several biological phenomena e.g., blood ow in arteries. However, developing numerical scheme for such problems is still very challenging even when the deformation of the boundary of domain is prescribed a priori. Possibility of excessive mesh distortion is one of the major challenge when solving such problems with numerical methods using boundary tted meshes. The arbitrary Lagrangian- Eulerian (ALE) approach is a way to overcome this difficulty. Numerical simulations of convection-dominated problems have for long been the subject to many researchers. Galerkin formulations, which yield the best approximations for differential equations with high diffusivity, tend to induce spurious oscillations in the numerical solution of convection dominated equations. Though such spurious oscillations can be avoided by adaptive meshing, which is computationally very expensive on ne grids. Alternatively, stabilization methods can be used to suppress the spurious oscillations. In this work, the considered equation is designed within the framework of ALE formulation. In the first part, Streamline Upwind Petrov-Galerkin (SUPG) finite element method with conservative ALE formulation is proposed. Further, the first order backward Euler and the second order Crank-Nicolson methods are used for the temporal discretization. It is shown that the stability of the semi-discrete (continuous in time) ALE-SUPG equation is independent of the mesh velocity, whereas the stability of the fully discrete problem is unconditionally stable for implicit Euler method and is only conditionally stable for Crank-Nicolson time discretization. Numerical results are presented to support the stability estimates and to show the influence of the SUPG stabilization parameter in a time-dependent domain. In the second part of this work, SUPG stabilization method with non-conservative ALE formulation is proposed. The implicit Euler, Crank-Nicolson and backward difference methods are used for the temporal discretization. At the discrete level in time, the ALE map influences the stability of the corresponding discrete scheme with different time discretizations, and it leads to schemes where conservative and non-conservative formulations are no longer equivalent. The stability of the fully discrete scheme, irrespective of the temporal discretization, is only conditionally stable. It is observed from numerical results that the Crank-Nicolson scheme induces high oscillations in the numerical solution compare to the implicit Euler and the backward difference time discretiza-tions. Moreover, the backward difference scheme is more sensitive to the stabilization parameter k than the other time discretizations. Further, the difference between the solutions obtained with the conservative and non-conservative ALE forms is significant when the deformation of domain is large, whereas it is negligible in domains with small deformation. Finally, the local projection stabilization (LPS) and the higher order dG time stepping scheme are studied for convection dominated problems. The analysis is based on the quadrature formula for approximating the integrals in time. We considered the exact integration in time, which is impractical to implement and the Radau quadrature in time, which can be used in practice. The stability and error estimates are shown for the mathematical basis of considered numerical scheme with both time integration methods. The numerical analysis reveals that the proposed stabilized scheme with exact integration in time is unconditionally stable, whereas Radau quadrature in time is conditionally stable with time-step restriction depending on the ALE map. The theoretical estimates are illustrated with appropriate numerical examples with distinct features. The second order dG(1) time discretization is unconditionally stable while Crank-Nicolson gives the conditional stable estimates only. The convergence order for dG(1) is two which supports the error estimate.

<p align="justify">This work deals with the development of numerical methods for compressible flow simulation with application to the interaction of fluid flows and structural bodies.</p> <p align="justify">First, we develop numerical methods based on multidimensional upwind residual distribution (RD) schemes. Theoretical results for the stability and accuracy of the methods are given. Then, the RD schemes for unsteady problems are extended for computations on moving meshes. As a second approach, cell centered and vertex centered finite volume (FV) schemes are considered. The RD schemes are compared to FV schemes by means of the 1D modified equation and by the comparison of the numerical results for scalar problems and system of Euler equations. We present a number of two and three dimensional steady and unsteady test cases, illustrating properties of the numerical methods. The results are compared with the theoretical solution and experimental data.</p> <p align="justify">In the second part, a numerical method for fluid-structure interaction problems is developed. The problem is divided into three distinct sub-problems: Computational Fluid Dynamics, Computational Solid Mechanics and the problem of fluid mesh movement. The problem of Computational Solid Mechanics is formulated as a system of partial differential equations for an anisotropic elastic continuum and solved by the finite element method. The mesh movement is determined using the pseudo-elastic continuum approach and solved again by the finite element method. The coupling of the problems is achieved by a simple sub-iterative approach. Capabilities of the methods are demonstrated on computations of 2D supersonic panel flutter and 3D transonic flutter of the AGARD 445.6 wing. In the first case, the results are compared with the theoretical solution and the numerical computations given in the references. In the second case the comparison with experimental data is presented.</p>

Dans ce manuscrit, un modèle complet pour la simulation de l'écoulement d'un fluide thermor éactif à travers un milieu poreux fortement compressible est présenté. Ce modèle est utilisé pour l'étude des procédés d'élaboration des matériaux composites par infusion à travers leur épaisseur (Liquid Resin Infusion-LRI et Resin Film Infusion-RFI ). Dans ces procédés, le mélange entre les renforts et la résine liquide est réalisé dans la direction transverse aux plans des préformes pendant la phase de mise en forme. Les coˆuts sont ainsi réduits et les problèmes de remplissage éliminés. Ces procédés sont néanmoins peu maîtrisés et les caractéristiques de la pièce finale difficilement prévisibles (principalement les épaisseurs et les porosités). La mise au point d'un modèle numérique constituerait un bon outil pour développer et finaliser de nouvelles solutions composites. D'un point de vue physique, l'infusion de la résine à travers l'épaisseur des préformes est une conséquence de la pression appliquée sur l'empilement résine/préforme. Dans cette analyse multi-physique deux types de problèmes sont rencontrés. Tout d'abord, on connait mal les conditions de couplage entre les zones liquides, gouvernées par les équations de Stokes, et les préformes imprégnées assimilées à des milieux poreux, gouvernées par une loi de Darcy et une loi de comportement mécanique non-linéaire. Par ailleurs, les interactions entre l'écoulement de la résine et la compression des préformes ne sont pas bien maîtrisées. Le modèle développé inclut donc une condition de Beaver-Joseph- Schaffman modifiée pour le couplage entre les zones de Darcy et de Stokes. Une formulation ALE pour l'écoulement de la résine dans un milieu poreux déformable subissant de fortes déformations est utilisée et couplée à une formulation Lagrangienne Réactualisée pour la partie solide. Ces deux mécanismes physiques sont couplés à des modèles thermo-chimiques pour traiter la réticulation de la résine sous l'action du cycle de température. Dans ce travail, un certain nombre d'outils numériques et de nouvelles formulations ont été développés en vue de simuler les procédés LRI et RFI. Chaque outil est étudié et validé analytiquement ou numériquement avant d'être intégré dans les modèles LRI /RFI. Des simulations numériques d'infusion sont ensuite présentées et commentées, puis une première comparaison avec des essais expérimentaux est proposée.

This thesis presents a finite element methodology for simulating fluid–solid interactions in the left ventricle (LV) under LVAD support. The developed model was utilised to study the passive and active characteristics of ventricular function in anatomically accurate LV geometries constructed from normal and patient image data. A non–conforming ALE Navier–Stokes/finite–elasticity fluid–solid coupling system formed the core of the numerical scheme, onto which several novel numerical additions were made. These included a fictitious domain (FD) Lagrange multiplier method to capture the interactions between immersed rigid bodies and encasing elastic solids (required for the LVAD cannula), as well as modifications to the Newton–Raphson/line search algorithm (which provided a 2 to 10 fold reduction in simulation time). Additional developments involved methods for extending the model to ventricular simulations. This required the creation of coupling methods, for both fluid and solid problems, to enable the integration of a lumped parameter representation of the systemic and pulmonary circulatory networks; the implementation and tuning of models of passive and active myocardial behaviour; as well as the testing of appropriate element types for coupling non–conforming fluid– solid finite element models under high interface tractions (finding that curvilinear spatial interpolations of the fluid geometry perform best). The behaviour of the resulting numerical scheme was investigated in a series of canonical test problems and found to be convergent and stable. The FD convergence studies also found that discontinuous pressure elements were better at capturing pressure gradients across FD boundaries. The ventricular simulations focused firstly on studying the passive diastolic behaviour of the LV both with and without LVAD support. Substantially different vortical flow features were observed when LVAD outflow was included. Additionally, a study of LVAD cannula lengths, using a particle tracking algorithm to determine recirculation rates of blood within the LV, found that shorter cannulas improved the recirculation of blood from the LV apex. Incorporating myocardial contraction, the model was extended to simulate the full cardiac cycle, converging on a repeating pressure–volume loop over 2 heart beats. Studies on the normal LV geometry found that LVAD implementation restricts the recirculation of early diastolic inflow, and that fluid–solid coupled models introduce greater heterogeneity of myocardial work than was observed in equivalent solid only models. A patient study was undertaken using a myocardial geometry constructed using image data from an LVAD implant recipient. A series of different LVAD flow regimes were tested. It was found that the opening of the aortic valve had a homogenising effect on the spatial variation of work, indicating that the synchronisation of LVAD outflow with the cardiac cycle is more important if the valve remains shut. Additionally, increasing LVAD outflow during systole and decreasing it during diastole led to improved mixing of blood in the ventricular cavity – compared with either the inverse, or holding outflow constant. Validation of these findings has the potential to impact the treatment protocols of LVAD patients.

Este trabalho consiste no desenvolvimento de ferramentas computacionais para análise não linear geométrica de interação fluido-casca utilizando o Método dos Elementos Finitos (MEF). O algoritmo para dinâmica dos fluidos é explícito e a integração temporal é baseada em linhas características. O código computacional é capaz de simular as equações de Navier-Stokes para escoamentos compressíveis tanto na descrição Euleriana como na descrição Lagrangeana-Euleriana arbitrária (ALE), na qual é possível prescrever movimentos para a malha do fluido. A estrutura é modelada em descrição Lagrangeana total através de uma formulação de MEF para análise dinâmica não linear geométrica de cascas baseada no teorema da mínima energia potencial total escrito em função das posições nodais e vetores generalizados e não em deslocamentos e rotações. Essa característica evita o uso de aproximações de grandes rotações. Dois modelos de acoplamentos são desenvolvidos. O primeiro modelo, ideal para problemas onde a escala de deslocamentos não é muito grande comparada com as dimensões do domínio do fluido, é baseado na descrição ALE e o acoplamento entre as duas diferentes malhas é feito através do mapeamento das posições locais dos nós do contorno do fluido sobre os elementos de casca e vice-versa, evitando a necessidade de coincidência entre os nós da casca e do fluido. A malha do fluido é adaptada dinamicamente usando um procedimento simples baseado nas posições e velocidades nodais da casca. O segundo modelo de acoplamento, ideal para problemas com grande escala de deslocamentos tais como estruturas infláveis, considera a casca imersa na malha do fluido e consiste em um procedimento robusto baseado em curvas de nível da função distância assinalada do contorno, o qual integra o algoritmo Lagrangeano de casca com o Fluido em descrição Euleriana, sem necessidade de movimentação da malha do fluido, onde a representação computacional do fluido se resume a uma malha não estruturada maior ou igual ao domínio inicial do fluido e a interface fluido-casca dentro da malha do fluido é identificada por meio de curvas de nível da função distância assinalada do contorno. Ambos os modelos são testados através de exemplos numéricos mostrando robustez e eficiência. Finalmente, como uma sugestão para o futuro desenvolvimento desta pesquisa, iniciaram-se estudos relativos a funções B-splines. O uso desse tipo de funções deverá resolver problemas de estabilidade relativos a oscilações espúrias devidas ao uso de polinômios de Lagrange para a representação de descontinuidades. / This work consists of the development of computational tools for nonlinear geometric fluid-shell interaction analysis using the Finite Element Method (FEM). The fluid solver is explicit and its time integration based on characteristics. The computational code is able to simulate the Navier-Stokes equations for compressible flows written in the Eulerian description as well as in the arbitrary Lagrangian-Eulerian (ALE) description, enabling movements prescription for the fluid mesh. The structure is modeled in a total Lagrangian description, using a FEM formulation to deal with geometrical nonlinear dynamics of shells based on the minimum potential energy theorem written regarding nodal positions and generalized unconstrained vectors, not displacements and rotations, avoiding the use of large rotation approximations. Two partitioned coupling models are developed. The first model, ideal for simulations where the displacements scale is not very large compared to the fluid domain, is based on the ALE description and the coupling between the two different meshes is done by mapping the fluid boundary nodes local positions over the shell elements and vice-versa, avoiding the need for matching fluid and shell nodes. The fluid mesh is adapted using a simple approach based on shell nodal positions and velocities. The second model, ideal for problems with large scales of displacements such as inflatable structures, is based on immersed boundary and consists of a robust level-set based approach that integrates the Lagrangian shell finite and the Eulerian finite element high speed fluid flow solver, with no need for mesh adaptation, where the fluid representation relies on a fixed unstructured mesh larger or equal to the initial fluid domain and the fluid-shell interface inside the fluid mesh is tracked with level sets of a boundary signed distance function. Both models are tested with numerical examples, showing efficiency and robustness. Finally, as a suggestion for future development of this research, we started studies relatives to B-Spline functions. The use of this kind of functions should solve stability problems related to spurious oscillations due to the use of Lagrange polynomials for representing discontinuities.