Fluid flow network simulation codes are commonly used as a design and analysis tool for many engineering problems such as gas distribution networks, power plants and heat pumps. Two formulations of conservation of momentum have been widely applied in fluid flow network simulation models namely those based on static pressure and those based on total pressure. The total pressure formulations are convenient in that they eliminate the difficulties associated with the calculation of the convective terms and components such as pipe junctions are treated in a straightforward manner based on total pressure losses. However, the different formulations of total pressure for compressible and incompressible flow require different formulations of the momentum conservation equation, which is inconvenient for implementation in a generic network simulation code. In this thesis a united total pressure formulation is first derived which is valid for all fluids and therefore eliminates the inconvenience of switching between the compressible and incompressible formulations. A non-iterative method for the solution of the non-isothermal discretised equations based on the total pressure formulation is then introduced and consistency is illustrated. The method appears to be very stable for subsonic flows, while rapid steady state convergence is observed. A systematic comparison is also done with traditional static pressure based methods and the similarities and differences between the two formulations are illuminated. The different time scales involved in the simulation of transient flow in fluid networks are problematic when conventional fixed time step methods are used for time-wise integration. The time scales associated with acoustic and kinematic wave phenomena as well as storage effects can differ by orders in magnitude. This thesis also presents a simple adaptive time step algorithm which can be readily used in conjunction with all the commonly used first order methods for fluid flow networks. Two test problems are selected to demonstrate the efficiency and savings obtained with this procedure. The adaptive time step algorithm correctly selects appropriate time steps for all phenomena and significant computational savings are observed for accurate integration. In addition, a procedure is implemented which automatically selects the appropriate integration method. The resulting algorithm is a fully adaptive algorithm which switches between a fully implicit method and a semi-implicit method. Two test problems are once again used to demonstrate the efficiency and savings. The fully adaptive algorithm correctly selects appropriate methods for all phenomena and significant additional computational savings are observed. / Thesis (PhD (Mechanical Engineering))--North-West University, Potchefstroom Campus, 2013

Fluid flow network simulation codes are commonly used as a design and analysis tool for many engineering problems such as gas distribution networks, power plants and heat pumps. Two formulations of conservation of momentum have been widely applied in fluid flow network simulation models namely those based on static pressure and those based on total pressure. The total pressure formulations are convenient in that they eliminate the difficulties associated with the calculation of the convective terms and components such as pipe junctions are treated in a straightforward manner based on total pressure losses. However, the different formulations of total pressure for compressible and incompressible flow require different formulations of the momentum conservation equation, which is inconvenient for implementation in a generic network simulation code. In this thesis a united total pressure formulation is first derived which is valid for all fluids and therefore eliminates the inconvenience of switching between the compressible and incompressible formulations. A non-iterative method for the solution of the non-isothermal discretised equations based on the total pressure formulation is then introduced and consistency is illustrated. The method appears to be very stable for subsonic flows, while rapid steady state convergence is observed. A systematic comparison is also done with traditional static pressure based methods and the similarities and differences between the two formulations are illuminated. The different time scales involved in the simulation of transient flow in fluid networks are problematic when conventional fixed time step methods are used for time-wise integration. The time scales associated with acoustic and kinematic wave phenomena as well as storage effects can differ by orders in magnitude. This thesis also presents a simple adaptive time step algorithm which can be readily used in conjunction with all the commonly used first order methods for fluid flow networks. Two test problems are selected to demonstrate the efficiency and savings obtained with this procedure. The adaptive time step algorithm correctly selects appropriate time steps for all phenomena and significant computational savings are observed for accurate integration. In addition, a procedure is implemented which automatically selects the appropriate integration method. The resulting algorithm is a fully adaptive algorithm which switches between a fully implicit method and a semi-implicit method. Two test problems are once again used to demonstrate the efficiency and savings. The fully adaptive algorithm correctly selects appropriate methods for all phenomena and significant additional computational savings are observed. / Thesis (PhD (Mechanical Engineering))--North-West University, Potchefstroom Campus, 2013

Modelling in cardiac electrophysiology results in a complex system of partial differential equations (PDE) describing the propagation of the electrical wave in the heart muscle coupled with a highly nonlinear system of ordinary differential equations (ODE) describing the ionic activity in the cardiac cells. This system forms the widely accepted bidomain model or its slightly simpler version, the monodomain model. To a large extent, the stiffness of the whole model depends on the choice of the ionic model, which varies in terms of complexity and realism. These simulations require accurate and, depending on the ionic model used, possibly very stable numerical methods. At this time, solving these models numerically requires CPU time of around one day per heartbeat. Therefore, it is necessary to use the most efficient method for these simulations. This research focuses on the comparison and analysis of several time-stepping methods: explicit or semi-implicit, operator splitting, deferred correction and Rush-Larsen methods. The goal is to find the optimal method for the ionic model used. For our analysis, we used the monodomain model but our results apply to the bidomain model as well. We compare the methods for three ionic models of varying complexity and stiffness: the Mitchell-Schaeffer models with only 2 variables, the more realistic Beeler-Reuter model with 8 variables, and the stiff and very complex ten Tuscher-Noble-Noble-Panfilov (TNNP) models with 17 variables. For each method, we derived absolute stability criteria of the spatially discretized monodomain model and verified that the theoretical critical time steps obtained closely match the ones in numerical experiments. Convergence tests were also conducted to verify that the numerical methods achieve an optimal order of convergence on the model variables and derived quantities (such as speed of the wave, depolarization time), and this in spite of the local non-differentiability of some of the ionic models. We looked at the efficiency of the different methods by comparing computational times for similar accuracy. Conclusions are drawn on the methods to be used to solve the monodomain model based on the model stiffness and complexity, measured respectively by the most negative eigenvalue of the model's Jacobian and the number of variables, and based on strict stability and accuracy criteria.

In this thesis study, implementation, numerical stability and convergence rate issues of turbulence modeling are explored. For this purpose, a one equation turbulence model, Spalart-Allmaras, and a two-equation turbulence model, SST k-w, are adapted to an explicit, cell centered, finite volume method based, structured / hybrid multi grid flow solver, SENSE2D, developed at TUBITAK-SAGE. Governing equations for both the flow and the turbulence are solved in a loosely coupled manner, however, each set of equations are solved using a coupled, semi-implicit solution algorithm. In multigrid solutions, the semi-implicit solution algorithm and the turbulence model equations are employed only in the finest level grid. As a result, stable and convergent numerical solutions are obtained. In order to validate the turbulence models and the semi-implicit solution algorithm implemented, turbulent flow solutions over a flat plate, RAE2822 airfoil and NLR7301 multi element airfoil are performed. The results are compared with the experimental data and the numerical results of the commercial CFD package FLUENT. It is shown that the numerical results obtained by SENSE2D are in good agreement with the experimental data and the FLUENT results. In addition to the turbulence modeling studies, convergence rate studies are also performed by multigrid and semi-implicit solution methods. It is shown that, the convergence rates of the semi-implicit solutions are increased about 5 times for single grid and 35% for multigrid solutions in comparison to the explicit solutions.

The construction of compact high-order Residual Distribution schemes for the discretizationof steady multidimensional advection-diffusion problems on unstructuredgrids is presented. Linear and non-linear scheme are considered. A piecewise continuouspolynomial approximation of the solution is adopted and a gradient reconstructionprocedure is used in order to have a continuous representation of both thenumerical solution and its gradient. It is shown that the gradient must be reconstructedwith the same accuracy of the solution, otherwise the formal accuracy ofthe numerical scheme is lost in applications in which diffusive effects prevail overthe advective ones, and when advection and diffusion are equally important. Thenthe method is extended to systems of equations, with particular emphasis on theNavier-Stokes and RANS equations. The accuracy, efficiency, and robustness of theimplicit RD solver is demonstrated using a variety of challenging aerodynamic testproblems.

Sea ice plays a key role in the global climate system. Indeed, through the albedo effect it reflects significant solar radiation away from the oceans, while it also plays a key role in the momentum and heat transfer between the atmosphere and ocean by acting as an insulating layer between the two. Furthermore, as more sea ice melts due to climate change, additional fresh water is released into the upper oceans, affecting the global circulation of the ocean as a whole. While there has been significant effort in recent decades, the ability to simulate sea ice has lagged behind other components of the climate system and most Earth System Models fail to capture the observed losses of Arctic sea ice, which is largely attributed to our inability to resolve sea ice dynamics. The most widely accepted model for sea ice dynamics is the Viscous- Plastic (VP) rheology, which leads to a very non-linear set of partial differential equations that are known to be intrinsically difficult to solve numerically. This work builds on recent advances in solving these equations with a Jacobian-Free Newton- Krylov (JFNK) solver. We present an improved JFNK solver, where a fully second order discretization is achieved via the Crank Nicolson scheme and consistency is improved via a novel approach to the rheology term. More importantly, we present a significant improvement to the Jacobian approximation used in the Newton iterations, and partially form the action of the matrix by expressing the linear and nearly linear terms in closed form and approximating the remaining highly non-linear term with a second order approximation of its Gateaux derivative. This is in contrast with the previous approach which used a first order approximation for the Gateaux derivative of the whole functional. Numerical tests on synthetic equations confirm the theoretical convergence rate and demonstrate the drastic improvements seen by using a second order approximation in the Gateaux derivative. To produce a fast and efficient solver for VP sea ice dynamics, the improved JFNK solver is then coupled with a non- oscillatory, central differencing scheme for transporting sea ice as well as a novel method for tracking the ice domain using a level set method. Two idealized test cases are then presented and simulation results discussed, demonstrating the solver’s ability to efficiently produce Viscous-Plastic, physically motivated solutions. / Graduate

Partial differential equations (PDEs) and stochastic partial differential equations (SPDEs) are powerful tools in modeling real-world phenomena in many fields such as geo-engineering. For instance processes such as oil or gas recovery from hydrocarbon reservoirs and mining heat from geothermal reservoirs can be modeled by PDEs or SPDEs. An important task is to understand the behavior of such phenomena. This can be achieved through explicit solutions of equations. Since explicit solutions of many PDEs and SPDEs are rarely known, developing numerical schemes is a good alternative to provide approximations of these explicit solutions. The study of numerical solutions of PDEs and SPDEs is therefore an active research area and has attracted a lot of attentions since at least two decades. The aims of this dissertation is to develop numerical schemes to approximate semilinear parabolic PDEs and SPDEs in space and in time. The approximation in space is done via the standard Galerkin finite element method and the approximation in time, which is the core of our work is done via various numerical integrators. This dissertation consists of two general parts. The first part of this thesis deals with autonomous PDEs and SPDEs. Here, our main interest is on semilinear PDEs and SPDEs where the nonlinear part is stronger than the linear part also called (stochastic) reactive dominated transport equations, or stiff problems. For such problems, many numerical techniques in the current scientific literature lose their good stability properties. We develop a new explicit exponential integrator (called exponential Rosenbrock-type method) and a new semi-implicit method (called linear implicit Rosenbrock-type method), appropriate for such PDEs and SPDEs. We analyze the strong convergence of our novel fully discrete schemes towards the mild solution of the (S)PDE and obtain convergence orders similar to existing ones in the literature. The second part of this thesis focuses on numerical approximations of semilinear non-autonomous parabolic PDEs and SPDEs. Such equations are more realistic than the autonomous ones and find applications in many fields such as fluid mechanics, quantum field theory, electromagnetism, etc. Numerics of non-autonomous semilinear parabolic PDEs and SPDEs are far from being well understood in the literature. We fill that gap in this thesis by developing a new exponential integrator (called Magnus-type method) and the semi-implicit method for such problems and provide their strong convergence towards the mild solution. Moreover, for both autonomous and non-autonomous SPDEs driven by additive noise, we achieve optimal convergence order in time 1 or approximately 1. Numerical simulations are provided to illustrate our theoretical findings in both autonomous and non-autonomous cases.

Este é o primeiro trabalho que apresenta simulações tridimensionais completamente adaptativas de um modelo de campo de fase para um fluido incompressível com densidade de massa constante e viscosidade variável, conhecido como Modelo H. Solucionando numericamente as equações desse modelo em malhas refinadas localmente com a técnica AMR, simulamos computacionalmente escoamentos bifásicos tridimensionais. Os modelos de campo de fase oferecem uma aproximação física sistemática para investigar fenômenos que envolvem sistemas multifásicos complexos, tais como fluidos com camadas de mistura, a separação de fases sob forças de cisalhamento e a evolução de micro-estruturas durante processos de solidificação. Como as interfaces são substituídas por delgadas regiões de transição (interfaces difusivas), as simulações de campo de fase requerem muita resolução nessas regiões para capturar corretamente a física do problema em estudo. Porém essa não é uma tarefa fácil de ser executada numericamente. As equações que caracterizam o modelo de campo de fase contêm derivadas de ordem elevada e intrincados termos não lineares, o que exige uma estratégia numérica eficiente capaz de fornecer precisão tanto no tempo quanto no espaço, especialmente em três dimensões. Para obter a resolução exigida no tempo, usamos uma discretização semi-implícita de segunda ordem para solucionar as equações acopladas de Cahn-Hilliard e Navier-Stokes (Modelo H). Para resolver adequadamente as escalas físicas relevantes no espaço, utilizamos malhas refinadas localmente que se adaptam dinamicamente para recobrir as regiões de interesse do escoamento, como por exemplo, as vizinhanças das interfaces do fluido. Demonstramos a eficiência e a robustez de nossa metodologia com simulações que incluem a separação dos componentes de uma mistura bifásica, a deformação de gotas sob cisalhamento e as instabilidades de Kelvin-Helmholtz. / This is the first work that introduces 3D fully adaptive simulations for a phase field model of an incompressible fluid with matched densities and variable viscosity, known as Model H. Solving numerically the equations of this model in meshes locally refined with AMR technique, we simulate computationally tridimensional two-phase flows. Phase field models offer a systematic physical approach to investigate complex multiphase systems phenomena such as fluid mixing layers, phase separation under shear and microstructure evolution during solidification processes. As interfaces are replaced by thin transition regions (diffuse interfaces), phase field simulations need great resolution in these regions to capture correctly the physics of the studied problem. However, this is not an easy task to do numerically. Phase field model equations have high order derivatives and intricate nonlinear terms, which require an efficient numerical strategy that can achieve accuracy both in time and in space, especially in three dimensions. To obtain the required resolution in time, we employ a semi-implicit second order discretization scheme to solve the coupled Cahn-Hilliard/Navier-Stokes equations (Model H). To resolve adequatly the relevant physical scales in space, we use locally refined meshes which adapt dynamically to cover special flow regions, e.g., the vicinity of the fluid interfaces. We demonstrate the efficiency and robustness of our methodology with simulations that include spinodal decomposition, the deformation of drops under shear and Kelvin-Helmholtz instabilities.

This dissertation concerns the model reduction of linear periodic descriptor systems both in continuous and discrete-time case. In this dissertation, mainly the projection based approaches are considered for model order reduction of linear periodic time varying descriptor systems. Krylov based projection method is used for large continuous-time periodic descriptor systems and balancing based projection technique is applied to large sparse discrete-time periodic descriptor systems to generate the reduce systems. For very large dimensional state space systems, both the techniques produce large dimensional solutions. Hence, a recycling technique is used in Krylov based projection methods which helps to compute low rank solutions of the state space systems and also accelerate the computational convergence. The outline of the proposed model order reduction procedure is given with more details. The accuracy and suitability of the proposed method is demonstrated through different examples of different orders. Model reduction techniques based on balance truncation require to solve matrix equations. For periodic time-varying descriptor systems, these matrix equations are projected generalized periodic Lyapunov equations and the solutions are also time-varying. The cyclic lifted representation of the periodic time-varying descriptor systems is considered in this dissertation and the resulting lifted projected Lyapunov equations are solved to achieve the periodic reachability and observability Gramians of the original periodic systems. The main advantage of this solution technique is that the cyclic structures of projected Lyapunov equations can handle the time-varying dimensions as well as the singularity of the period matrix pairs very easily. One can also exploit the theory of time-invariant systems for the control of periodic ones, provided that the results achieved can be easily re-interpreted in the periodic framework. Since the dimension of cyclic lifted system becomes very high for large dimensional periodic systems, one needs to solve the very large scale periodic Lyapunov equations which also generate very large dimensional solutions. Hence iterative techniques, which are the generalization and modification of alternating directions implicit (ADI) method and generalized Smith method, are implemented to obtain low rank Cholesky factors of the solutions of the periodic Lyapunov equations. Also the application of the solvers in balancing-based model reduction of discrete-time periodic descriptor systems is discussed. Numerical results are given to illustrate the effciency and accuracy of the proposed methods.

Coordenação de Aperfeiçoamento de Pessoal de Nível Superior / Em uma grande gama de problemas físicos, governados por equações diferenciais, muitas vezes é de interesse obter-se soluções para o regime transiente e, portanto, deve-se empregar técnicas de integração temporal. Uma primeira possibilidade seria a de aplicar-se métodos explícitos, devido à sua simplicidade e eficiência computacional. Entretanto, esses métodos frequentemente são somente condicionalmente estáveis e estão sujeitos a severas restrições na escolha do passo no tempo. Para problemas advectivos, governados por equações hiperbólicas, esta restrição é conhecida como a condição de Courant-Friedrichs-Lewy (CFL). Quando temse a necessidade de obter soluções numéricas para grandes períodos de tempo, ou quando o custo computacional a cada passo é elevado, esta condição torna-se um empecilho. A fim de contornar esta restrição, métodos implícitos, que são geralmente incondicionalmente estáveis, são utilizados. Neste trabalho, foram aplicadas algumas formulações implícitas para a integração temporal no método Smoothed Particle Hydrodynamics (SPH) de modo a possibilitar o uso de maiores incrementos de tempo e uma forte estabilidade no processo de marcha temporal. Devido ao alto custo computacional exigido pela busca das partículas a cada passo no tempo, esta implementação só será viável se forem aplicados algoritmos eficientes para o tipo de estrutura matricial considerada, tais como os métodos do subespaço de Krylov. Portanto, fez-se um estudo para a escolha apropriada dos métodos que mais se adequavam a este problema, sendo os escolhidos os métodos Bi-Conjugate Gradient (BiCG), o Bi-Conjugate Gradient Stabilized (BiCGSTAB) e o Quasi-Minimal Residual (QMR). Alguns problemas testes foram utilizados a fim de validar as soluções numéricas obtidas com a versão implícita do método SPH. / In a wide range of physical problems governed by differential equations, it is often of interest to obtain solutions for the unsteady state and therefore it must be employed temporal integration techniques. One possibility could be the use of an explicit methods due to its simplicity and computational efficiency. However, these methods are often only conditionally stable and are subject to severe restrictions for the time step choice. For advective problems governed by hyperbolic equations, this restriction is known as the Courant-Friedrichs-Lewy (CFL) condition. When there is the need to obtain numerical solutions for long periods of time, or when the computational cost for each time step is high, this condition becomes a handicap. In order to overcome this restriction implicit methods can be used, which are generally unconditionally stable. In this study, some implicit formulations for time integration are used in the Smoothed Particle Hydrodynamics (SPH) method to enable the use of larger time increments and obtain a strong stability in the time evolution process. Due to the high computational cost required by the particles tracking at each time step, the implementation will be feasible only if efficient algorithms were applied for this type of matrix structure such as Krylov subspace methods. Therefore, we carried out a study for the appropriate choice of methods best suited to this problem, and the methods chosen were the Bi-Conjugate Gradient (BiCG), the Bi-Conjugate Gradient Stabilized (BiCGSTAB) and the Quasi-Minimal Residual(QMR). Some test problems were used to validate the numerical solutions obtained with the implicit version of the SPH method.