Search

Mathematical models of the cardiovascular system have come a long way since they were first introduced in the early 19th century. Driven by a rapid development of experimental techniques, numerical methods, and computer hardware, detailed models that describe physical scales from the molecular level up to organs and organ systems have been derived and used for physiological research. Mathematical and computational models can be seen as condensed and quantitative formulations of extensive physiological knowledge and are used for formulating and testing hypotheses, interpreting and directing experimental research, and have contributed substantially to our understanding of cardiovascular physiology. However, in spite of the strengths of mathematics to precisely describe complex relationships and the obvious need for the mathematical and computational models to be informed by experimental data, there still exist considerable barriers between experimental and computational physiological research. In this review, we present a historical overview of the development of mathematical and computational models in cardiovascular physiology, including the current state of the art. We further argue why a tighter integration is needed between experimental and computational scientists in physiology, and point out important obstacles and challenges that must be overcome in order to fully realize the synergy of experimental and computational physiological research.

The triple decomposition of a velocity gradient tensor provides an analysis tool in fluid mechanics by which the flow can be split into a sum of irrotational straining flow, shear flow, and rigid body rotational flow. In 2007, Kolář formulated an optimization problem to compute the triple decomposition [V. Kolář, “Vortex identification: New requirements and limitations,” Int. J. Heat Fluid Flow 28, 638–652 (2007)], and more recently, the triple decomposition has been connected to the Schur form of the associated matrix. We show that the standardized real Schur form, which can be computed by state of the art linear algebra routines, is a solution to the optimization problem posed by Kolář. We also demonstrate why using the standardized variant of the real Schur form makes computation of the triple decomposition more efficient. Furthermore, we illustrate why different structures of the real Schur form correspond to different alignments of the coordinate system with the fluid flow and may, therefore, lead to differences in the resulting triple decomposition. Based on these results, we propose a new, simplified algorithm for computing the triple decomposition, which guarantees consistent results.

There is today a significant interest in harvesting renewable energy, specifically wind energy, in offshore and urban environments. Vertical axis wind turbines get increasing attention since they a ...

In this paper we address the challenges of adjoint based mesh adaptation in the setting of industrial CFD, characterized by complex geometry and high Reynolds numbers. We use a space-time finite element discretization of the NS equations, and we focus on the challenges of massively parallel algorithms and the mathematical formulation of the NS equations.

To date, it is unclear how fluid dynamics stimulate mechanosensory cells to induce an osteoprotective or osteodestructive response. We investigated how murine hematopoietic progenitor cells respond to 2 minutes of dynamic fluid flow stimulation with a precisely controlled sequence of fluid shear stresses. The response was quantified by measuring extracellular adenosine triphosphate (ATP), immunocytochemistry of Piezo1, and sarcoplasmic/endoplasmic Ca2+ reticulum ATPase 2 (SERCA2), and by the ability of soluble factors produced by mechanically stimulated cells to modulate osteoclast differentiation. We rejected our initial hypothesis that peak wall shear stress rate determines the response of hematopoietic progenitor cells to dynamic fluid shear stress, as it had only a minor correlation with the abovementioned parameters. Low stimulus amplitudes corresponded to activation of Piezo1, SERCA2, low concentrations of extracellular ATP, and inhibition of osteoclastogenesis and resorption area, while high amplitudes generally corresponded to osteodestructive responses. At a given amplitude (3 Pa) and waveform (square), the duration of individual stimuli (duty cycle) showed a strong correlation with the release of ATP and osteoclast number and resorption area. Collectively, our data suggest that hematopoietic progenitor cells respond in a viscoelastic manner to loading, since a combination of high shear stress amplitude and prolonged duty cycle is needed to trigger an osteodestructive response. Plain Language Summary: In case of painful joints or missing teeth, the current intervention is to replace them with an implant to keep a high-quality lifestyle. When exercising or chewing, the cells in the bone around the implant experience mechanical loading. This loading generally supports bone formation to strengthen the bone and prevent breaking, but can also stimulate bone loss when the mechanical loading becomes too high around orthopedic and dental implants. We still do not fully understand how cells in the bone can distinguish between mechanical loading that strengthens or weakens the bone. We cultured cells derived from the bone marrow in the laboratory to test whether the bone loss response depends on (i) how fast a mechanical load is applied (rate), (ii) how intense the mechanical load is (amplitude), or (iii) how long each individual loading stimulus is applied (duration). We mimicked mechanical loading as it occurs in the body, by applying very precisely controlled flow of fluid over the cells. We found that a mechanosensitive receptor Piezo1 was activated by a low amplitude stimulus, which usually strengthens the bone. The potential inhibitor of Piezo1, namely SERCA2, was only activated by a low amplitude stimulus. This happened regardless of the rate of application. At a constant high amplitude, a longer duration of the stimulus enhanced the bone-weakening response. Based on these results we deduce that a high loading amplitude tends to be bone weakening, and the longer this high amplitude persists, the worse it is for the bone.

We present a high performance computing framework for finite element simulation of blood flow in the left ventricle of the human heart. The mathematical model is described together with the discretization method and the parallel implementation in Unicorn which is part of the open source software framework FEniCS-HPC. We show results based on patient-specific data that capture essential features observed with other computational models and imaging techniques, and thus indicate that our framework possesses the potential to provide relevant clinical information for diagnosis and medical treatment. Several other studies have been conducted to simulate the three dimensional blood flow in the left ventricle of the human heart with prescribed wall movement. Our contribution to the field of cardiac research lies in establishing an open source framework modular both in modelling and numerical algorithms.

The International Energy Agency Technology Collaboration Programme for Ocean Energy Systems (OES) initiated the OES Wave Energy Conversion Modelling Task, which focused on the verification and validation of numerical models for simulating wave energy converters (WECs). The long-term goal is to assess the accuracy of and establish confidence in the use of numerical models used in design as well as power performance assessment of WECs. To establish this confidence, the authors used different existing computational modelling tools to simulate given tasks to identify uncertainties related to simulation methodologies: (i) linear potential flow methods; (ii) weakly nonlinear Froude–Krylov methods; and (iii) fully nonlinear methods (fully nonlinear potential flow and Navier–Stokes models). This article summarizes the code-to-code task and code-to-experiment task that have been performed so far in this project, with a focus on investigating the impact of different levels of nonlinearities in the numerical models. Two different WECs were studied and simulated. The first was a heaving semi-submerged sphere, where free-decay tests and both regular and irregular wave cases were investigated in a code-to-code comparison. The second case was a heaving float corresponding to a physical model tested in a wave tank. We considered radiation, diffraction, and regular wave cases and compared quantities, such as the WEC motion, power output and hydrodynamic loading. The Danish partners acknowledge the support from the Danish Energy Agency through project 374 64017-05197. The Swedish partners were supported by the Swedish Energy Agency under Grants P44423-1 and 375 P44432-1. J.V.R. and G.G. acknowledge the support by Science Foundation Ireland under Grant 13/IA/1886. This research was made possible by support from U.S. the Department of Energy’s EERE Water Power Technologies Office. Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This work was authored (in part) by the National Renewable Energy Laboratory, operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. Funding was provided by the U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Wind Energy Technologies Office. The views expressed in this article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledge that the U.S. Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes.

For problems involving large deformations of thin structures, simulating fluid-structure interaction (FSI) remains challenging largely due to the need to balance computational feasibility, efficiency, and solution accuracy. Overlapping domain techniques have been introduced as a way to combine the fluid-solid mesh conformity, seen in moving-mesh methods, without the need for mesh smoothing or re-meshing, which is a core characteristic of fixed mesh approaches. In this work, we introduce a novel overlapping domain method based on a partition of unity approach. Unified function spaces are defined as a weighted sum of fields given on two overlapping meshes. The method is shown to achieve optimal convergence rates and to be stable for steady-state Stokes, Navier-Stokes, and ALE Navier-Stokes problems. Finally, we present results for FSI in the case of a 2D mock aortic valve simulation. These initial results point to the potential applicability of the method to a wide range of FSI applications, enabling boundary layer refinement and large deformations without the need for re-meshing or user-defined stabilization., 34 pages, 15 figure

International audience; The Bloch-Torrey partial differential equation can be used to describe the evolution of the transverse magnetization of the imaged sample under the influence of diffusion-encoding magnetic field gradients inside the MRI scanner. The integral of the magnetization inside a voxel gives the simulated diffusion MRI signal. This paper proposes a finite element discretization on manifolds in order to efficiently simulate the diffusion MRI signal in domains that have a thin layer or a thin tube geometrical structure. The variable thickness of the three-dimensional domains is included in the weak formulation established on the manifolds. We conducted a numerical study of the proposed approach by simulating the diffusion MRI signals from the extracellular space (a thin layer medium) and from neurons (a thin tube medium), comparing the results with the reference signals obtained using a standard three-dimensional finite element discretization. We show good agreements between the simulated signals using our proposed method and the reference signals for a wide range of diffusion MRI parameters. The approximation becomes better as the diffusion time increases. The method helps to significantly reduce the required simulation time, computational memory, and difficulties associated with mesh generation, thus opening the possibilities to simulating complicated structures at low cost for a better understanding of diffusion MRI in the brain.

The numerical simulation of the diffusion MRI signal arising from complex tissue micro-structures is helpful for understanding and interpreting imaging data as well as for designing and optimizing MRI sequences. The discretization of the Bloch-Torrey equation by finite elements is a more recently developed approach for this purpose, in contrast to random walk simulations, which has a longer history. While finite element discretization is more difficult to implement than random walk simulations, the approach benefits from a long history of theoretical and numerical developments by the mathematical and engineering communities. In particular, software packages for the automated solutions of partial differential equations using finite element discretization, such as FEniCS, are undergoing active support and development. However, because diffusion MRI simulation is a relatively new application area, there is still a gap between the simulation needs of the MRI community and the available tools provided by finite element software packages. In this paper, we address two potential difficulties in using FEniCS for diffusion MRI simulation. First, we simplified software installation by the use of FEniCS containers that are completely portable across multiple platforms. Second, we provide a portable simulation framework based on Python and whose code is open source. This simulation framework can be seamlessly integrated with cloud computing resources such as Google Colaboratory notebooks working on a web browser or with Google Cloud Platform with MPI parallelization. We show examples illustrating the accuracy, the computational times, and parallel computing capabilities. The framework contributes to reproducible science and open-source software in computational diffusion MRI with the hope that it will help to speed up method developments and stimulate research collaborations., La Caixa 2016

The Bloch–Torrey equation describes the evolution of the spin (usually water proton) magnetization under the influence of applied magnetic field gradients and is commonly used in numerical simulations for diffusion MRI and NMR. Microscopic heterogeneity inside the imaging voxel is modeled by interfaces inside the simulation domain, where a discontinuity in the magnetization across the interfaces is produced via a permeability coefficient on the interfaces. To avoid having to simulate on a computational domain that is the size of an entire imaging voxel, which is often much larger than the scale of the microscopic heterogeneity as well as the mean spin diffusion displacement, smaller representative volumes of the imaging medium can be used as the simulation domain. In this case, the exterior boundaries of a representative volume either must be far away from the initial positions of the spins or suitable boundary conditions must be found to allow the movement of spins across these exterior boundaries. Many approaches have been taken to solve the Bloch–Torrey equation but an efficient high performance computing framework is still missing. In this paper, we present formulations of the interface as well as the exterior boundary conditions that are computationally efficient and suitable for arbitrary order finite elements and parallelization. In particular, the formulations are based on the partition of unity concept which allows for a discontinuous solution across interfaces conforming with the mesh with weak enforcement of real (in the case of interior interfaces) and artificial (in the case of exterior boundaries) permeability conditions as well as an operator splitting for the exterior boundary conditions. The method is straightforward to implement and it is available in FEniCS for moderate-scale simulations and in FEniCS-HPC for large-scale simulations. The order of accuracy of the resulting method is validated in numerical tests and a good scalability is shown for the parallel implementation. We show that the simulated dMRI signals offer good approximations to reference signals in cases where the latter are available and we performed simulations for a realistic model of a neuron to show that the method can be used for complex geometries.

This work presents a direct comparison of unsteady, turbulent flow simulations with measurements performed using a Gulfstream G550 nose landing gear model. The experimental campaign, which was carried out by researchers from the NASA Langley Research Center, provided a series of detailed, well documented wind-tunnel measurements for comparison and validation of computational fluid dynamics (CFD) and computational aeroacoustics (CAA) methodologies. Several computational efforts were collected and presented at the Benchmark for Airframe Noise Computation workshops, BANC-I and II. For our simulations, we used a General Galerkin finite element method (G2), where no explicit subgrid model is used, and where the computational mesh is adaptively refined with respect to a posteriori estimates of the error in a quantity of interest, here the source term in Lighthill’s equation. The mesh is fully unstructured and the solution is time-resolved, which are key ingredients for solving problems of industrial relevance in the field of aeroacoustics. Moreover, we choose to model the boundary layers on the landing gear geometry with a free-slip condition for the velocity, which we previously observed to produce good results for external flows at high Reynolds numbers, and which considerably reduces the amount of cells required in the mesh. The comparisons presented here are an attempt to quantify the accuracy of our models, methods and assumptions; to that end, several results containing both time-averaged and unsteady flow quantities, always side by side with corresponding experimental values, are reported. The main finding is that we are able to simulate a complex, unsteady flow problem using a parameter-free methodology developed for high Reynolds numbers, external aerodynamics and aeroacoustics applications.

Due to advances in medical imaging, computational fluid dynamics algorithms and high performance computing, computer simulation is developing into an important tool for understanding the relationship between cardiovascular diseases and intraventricular blood flow. The field of cardiac flow simulation is challenging and highly interdisciplinary. We apply a computational framework for automated solutions of partial differential equations using Finite Element Methods where any mathematical description directly can be translated to code. This allows us to develop a cardiac model where specific properties of the heart such as fluid-structure interaction of the aortic valve can be added in a modular way without extensive efforts. In previous work, we simulated the blood flow in the left ventricle of the heart. In this paper, we extend this model by placing prototypes of both a native and a mechanical aortic valve in the outflow region of the left ventricle. Numerical simulation of the blood flow in the vicinity of the valve offers the possibility to improve the treatment of aortic valve diseases as aortic stenosis (narrowing of the valve opening) or regurgitation (leaking) and to optimize the design of prosthetic heart valves in a controlled and specific way. The fluid-structure interaction and contact problem are formulated in a unified continuum model using the conservation laws for mass and momentum and a phase function. The discretization is based on an Arbitrary Lagrangian-Eulerian space-time finite element method with streamline diffusion stabilization, and it is implemented in the open source software Unicorn which shows near optimal scaling up to thousands of cores. Computational results are presented to demonstrate the capability of our framework.

The locally adapted parametric finite element method for interface problems on triangular meshes

Our simulation methodology is referred to as Direct FEM Simulation (DFS), or General Galerkin (G2) and uses a finite element method (FEM) with piecewise linear approximation in space and time, and with numerical stabilization in the form of a weighted least squares method based on the residual. The incompressible Navier-Stokes Equations (NSE) are discretized directly, without applying any filter. Thus, the method does not result in Large Eddy Simulation (LES) filtered solutions, but is instead an approximation of a weak solution satisfying the weak form of the NSE. In G2 we have a posteriori error estimates for quantities of interest that can be expressed as functionals of a weak solution. These a posteriori error estimates, which form the basis for our adaptive mesh refinement algorithm, are based on the solution of an associated adjoint problem with a goal quantity (the aerodynamic forces in this work) as data, similarly to an optimal control problem. We provide references to related work below. The methodology and software have been previously validated for a number of turbulent flow benchmark problems, including one of the HiLiftPW-2 high Reynolds number cases. The DFS method is implemented in the Unicorn solver, which uses the open source software framework FEniCS-HPC, designed for automated solution of partial differential equations on massively parallel architectures using the FEM. In this chapter we present adaptive results from the Third AIAA High Lift Prediction Workshop in Denver, Colorado based on our DFS methodology and Unicorn/FEniCS-HPC software. We show that the methodology quantitavely and qualitatively captures the main features of the experiment - aerodynamic forces and the stall mechanism with a novel numerical tripping, with a much coarser mesh resolution and cheaper computational cost than the standard in the field.

Echocardiography is the most commonly used image modality in cardiology, assessing several aspects of cardiac viability. The importance of cardiac hemodynamics and 4D blood flow motion has recently been highlighted, however such assessment is still difficult using routine echo-imaging. Instead, combining imaging with computational fluid dynamics (CFD)-simulations has proven valuable, but only a few models have been applied clinically. In the following, patient-specific CFD-simulations from transthoracic dobutamin stress echocardiography have been used to analyze the left ventricular 4D blood flow in three subjects: two with normal and one with reduced left ventricular function. At each stress level, 4D-images were acquired using a GE Vivid E9 (4VD, 1.7MHz/3.3MHz) and velocity fields simulated using a presented pathway involving endocardial segmentation, valve position identification, and solution of the incompressible Navier-Stokes equation. Flow components defined as direct flow, delayed ejection flow, retained inflow, and residual volume were calculated by particle tracing using 4th-order Runge-Kutta integration. Additionally, systolic and diastolic average velocity fields were generated. Results indicated no major changes in average velocity fields for any of the subjects. For the two subjects with normal left ventricular function, increased direct flow, decreased delayed ejection flow, constant retained inflow, and a considerable drop in residual volume was seen at increasing stress. Contrary, for the subject with reduced left ventricular function, the delayed ejection flow increased whilst the retained inflow decreased at increasing stress levels. This feasibility study represents one of the first clinical applications of an echo-based patient-specific CFD-model at elevated stress levels, and highlights the potential of using echo-based models to capture highly transient flow events, as well as the ability of using simulation tools to study clinically complex phenomena. With larger patient studies planned for the future, and with the possibility of adding more anatomical features into the model framework, the current work demonstrates the potential of patient-specific CFD-models as a tool for quantifying 4D blood flow in the heart.

We present a framework for coupled multiphysics in computational fluid dynamics, targeting massively parallel systems. Our strategy is based on general problem formulations in the form of partial differential equations and the finite element method, which open for automation, and optimization of a set of fundamental algorithms. We describe these algorithms, including finite element matrix assembly, adaptive mesh refinement and mesh smoothing; and multiphysics coupling methodologies such as unified continuum fluid-structure interaction (FSI), and aeroacoustics by coupled acoustic analogies. The framework is implemented as FEniCS open source software components, optimized for massively parallel computing. Examples of applications are presented, including simulation of aeroacoustic noise generated by an airplane landing gear, simulation of the blood flow in the human heart, and simulation of the human voice organ.

The massive computational cost for resolving all scales in a turbulent flow makes a direct numerical simulation of the underlying Navier-Stokes equations impossible in most engineering applications. Recent advances in adaptive finite element methods offer a new powerful tool in Computational Fluid Dynamics (CFD). The computational cost for simulating turbulent flow can be minimized by adaptively resolution of the mesh, based on a posteriori error estimation. Such adaptive methods have previously been implemented for efficient serial computations, but the extension to an efficient parallel solver is a challenging task. This work concerns the development of an adaptive finite element method that enables efficient computation of time resolved approximations of turbulent flow for complex geometries with a posteriori error control. We present efficient data structures and data decomposition methods for distributed unstructured tetrahedral meshes. Our work also concerns an efficient parallelization of local mesh refinement methods such as recursive longest edge bisection, and the development of an a priori predictive dynamic load balancing method, based on a weighted dual graph. We also address the challenges of emerging supercomputer architectures with the development of new hybrid parallel programming models, combining traditional message passing with lightweight one-sided communication. Our implementation has proven to be both general and efficient, scaling up to more than twelve thousands cores.

We present a framework for adaptive finite element computation of turbulent flow and fluid-structure interaction, with focus on general algorithms that allow for complex geometry and deforming domains. We give basic models and finite element discretization methods, adaptive algorithms and strategies for e cient parallel implementation. To illustrate the capabilities of the computational framework, we show a number of application examples from aerodynamics, aero-acoustics, biomedicine and geophysics. The computational tools are free to download open source as Unicorn, and as a high performance branch of the finite element problem solving environment DOLFIN, both part of the FEniCS project

High-Reynolds number turbulent incompressible multiphase flow represents a large class of engineering problems of key relevance to society. Here we describe our work on modeling two such problems: 1. The Consorcio de Aguas Bilbao Bizkaia is constructing a new storm tank system with an automatic cleaning system, based on periodically pushing tank water out in a tunnel 2. In the framework of the collaboration between BCAM - Basque Center for Applied Mathematics and Tecnalia R & I, the interaction of the sea flow with a semi submersible floating offshore wind platform is computationally investigated. We study the MARIA' benchmark modeling breaking waves over objects in marine environments. Both of these problems are modeled in the the Direct FEM/General Galerkin methodology for turbulent incompressible variable-densitv flow 1,2 in the FEniCS software framework.

The combination of refined medical imaging techniques and computational fluid dynamics (CFD) models has enabled the study of complex flow behavior on a highly regional level. Recently, we have developed a platform for patient-specific CFD modelling of blood flow in the left ventricle (LV), with input data and required boundary conditions acquired from 4D echocardiography. The platform robustness has been evaluated with respect to input variable variations, but for any clinical implementation model flow validation is essential. Therefore, the aim of this study is to evaluate the accuracy of the patient-specific CFD model against multimodal image-based flow measurements. For the validation, 4D echocardiography was acquired from two healthy subjects, from which LV velocity fields were simulated. In-vivo flows from the same two subjects were then acquired by pulsed wave (PW) Doppler imaging over both LV-valves, and by cine phase-contract magnetic resonance imaging (PC-MRI) at eight defined anatomical planes in the LV. By fusing PC-MRI and the ultrasound acquisitions using a three-chamber alignment algorithm, simulated and measured flows were quantitatively compared. General flow pattern correspondence was observed, with a mean error of 1.4 cm/s and root mean square deviation of 5.7 cm/s for all measured PC-MRI LV-planes. For the PW-Doppler comparison, a mean error of 3.6 cm/s was reported. Overall, the following work represents a validation of the proposed patient-specific CFD platform, and the agreement with clinical data highlight the potential for future clinical use of the models.

In this paper, we address a 3D fluid-structure interaction benchmark problem that represents important characteristics of biomedical modeling. We present a goal-oriented adaptive finite element methodology for incompressible fluid-structure interaction based on a streamline diffusion-type stabilization of the balance equations for mass and momentum for the entire continuum in the domain, which is implemented in the Unicorn/FEniCS software framework. A phase marker function and its corresponding transport equation are introduced to select the constitutive law, where the mesh tracks the discontinuous fluid-structure interface. This results in a unified simulation method for fluids and structures. We present detailed results for the benchmark problem compared with experiments, together with a mesh convergence study.

Developing multiphysics finite element methods (FEM) and scalable HPC implementations can be very challenging in terms of software complexity and performance, even more so with the addition of goal-oriented adaptive mesh refinement. To manage the complexity we in this work present general adaptive stabilized methods with automated implementation in the FEniCS-HPC automated open source software framework. This allows taking the weak form of a partial differential equation (PDE) as input in near-mathematical notation and automatically generating the low-level implementation source code and auxiliary equations and quantities necessary for the adaptivity. We demonstrate new optimal strong scaling results for the whole adaptive framework applied to turbulent flow on massively parallel architectures down to 25000 vertices per core with ca. 5000 cores with the MPI-based PETSc backend and for assembly down to 500 vertices per core with ca. 20000 cores with the PGAS-based JANPACK backend. As a demonstration of the power of the combination of the scalability together with the adaptive methodology allowing prediction of gross quantities in turbulent flow we present an application in aerodynamics of a full DLR-F11 aircraft in connection with the HiLift-PW2 benchmarking workshop with good match to experiments.

We present a new mathematical theory explaining the fluid mechanics of subsonic flight, which is fundamentally different from the existing boundary layer-circulation theory by Prandtl–Kutta–Zhukovsky formed 100 year ago. The new theory is based on our new resolution of d’Alembert’s paradox showing that slightly viscous bluff body flow can be viewed as zero-drag/lift potential flow modified by 3d rotational slip separation arising from a specific separation instability of potential flow, into turbulent flow with nonzero drag/lift. For a wing this separation mechanism maintains the large lift of potential flow generated at the leading edge at the price of small drag, resulting in a lift to drag quotient of size 15–20 for a small propeller plane at cruising speed with Reynolds number $Re \approx 10^7$ and a jumbojet at take-off and landing with $Re \approx 10^8$, which allows flight at affordable power. The new mathematical theory is supported by computed turbulent solutions of the Navier–Stokes equations with a slip boundary condition as a model of observed small skin friction of a turbulent boundary layer always arising for $Re > 10^6$, in close accordance with experimental observations over the entire range of angle of attacks including stall using a few millions of mesh points for a full wing-body configuration.

In diffusion nuclear magnetic resonance (NMR) and diffusion magnetic resonance imaging (MRI), the multi-compartment Bloch-Torrey equation plays an important role in probing the diffusion characteristics from a nanometer scale to a macroscopic scale. The signal attenuation can be computed by solving the equation. If the volume of interest is composed by multiple compartments, interface conditions with permeability are imposed. Depending on applications, different gradient strengths can be used to capture the signal attenuation. In probing porous media, for instance, high gradient strengths are used. In diffusion MRI, since water molecules enter and exit the computational domain in realistic cases, pseudo-periodic boundary conditions are used. These conditions cause difficulties in solving the equation efficiently and many efforts have been made to develop an efficient numerical method. However, large-scale problems for supercomputers with realistic applications have not been considered yet. We propose a framework for the multi-compartment Bloch-Torrey models based on the FEniCS-HPC platform, a part of the FEniCS project that allows for automated discretization, automated error control with mesh adaptivity and high performance computing. The framework runs on supercomputers with near optimal weak and strong scaling. Our work includes two parts. First, we simplify the multi-compartment Bloch-Torrey model used in diffusion MRI by proposing an approximation to the pseudo-periodic boundary conditions to derive a general form for the interface and boundary conditions. The second part includes implementation and numerical validation of our method on the FEniCS-HPC platform. This simplified model is straightforward to implement and to parallelize and shows promise in validation against more realistic models.

This work develops finite element methods with high order stabilization, and robust and efficient adaptive algorithms for Large Eddy Simulation of turbulent compressible flows. The equations are approximated by continuous piecewise linear functions in space, and the time discretization is done in implicit/explicit fashion: the second order Crank-Nicholson method and third/fourth order explicit Runge-Kutta methods. The full residual of the system and the entropy residual, are used in the construction of the stabilization terms. These methods are consistent for the exact solution, conserves all the quantities, such as mass, momentum and energy, is accurate and very simple to implement. We prove convergence of the method for scalar conservation laws in the case of an implicit scheme. The convergence analysis is based on showing that the approximation is uniformly bounded, weakly consistent with all entropy inequalities, and strongly consistent with the initial data. The convergence of the explicit schemes is tested in numerical examples in 1D, 2D and 3D. To resolve the small scales of the flow, such as turbulence fluctuations, shocks, discontinuities and acoustic waves, the simulation needs very fine meshes. In this thesis, a robust adjoint based adaptive algorithm is developed for the time-dependent compressible Euler/Navier-Stokes equations. The adaptation is driven by the minimization of the error in quantities of interest such as stresses, drag and lift forces, or the mean value of some quantity. The implementation and analysis are validated in computational tests, both with respect to the stabilization and the duality based adaptation.

In this paper we present a computational study of the stability of time dependent dual problems for compressible flow at high Reynolds numbers in two dimensions. The dual problem measures the sensitivity of an output functional with respect to numerical errors and is a key part of goal oriented a posteriori error estimation. Our investigation shows that the dual problem associated with the computation of the drag force for the compressible Euler/Navier--Stokes equations, which are approximated numerically using different temporal discretization and stabilization techniques, is unstable and exhibits blow-up for several Mach regimes considered in this paper.

In this paper we first review our recent work on a new framework for adaptive turbulence simulation: we model turbulence by weak solutions to the Navier-Stokes equations that are wellposed with respect to mean value output in the form of functionals, and we use an adaptive finite element method to compute approximations with a posteriori error control based on the error in the functional output. We then derive a local energy estimate for a particular finite element method, which we connect to related work on dissipative weak Euler solutions with kinetic energy dissipation due to lack of local smoothness of the weak solutions. The ideas are illustrated by numerical results, where we observe a law of finite dissipation with respect to a decreasing mesh size.

In this paper, we describe an incompressible Unified Continuum (UC) model in Euler (laboratory) coordinates with a moving mesh for tracking the fluid-structure interface as part of the discretization, allowing simple and general formulation and efficient computation. The model consists of conservation equations for mass and momentum, a phase convection equation and a Cauchy stress and phase variable θ as data for defining material properties and constitutive laws. We target realistic 3D turbulent fluid-structure interaction (FSI) applications, where we show simulation results of a flexible flag mounted in the turbulent wake behind a cube as a qualitative test of the method, and quantitative results for 2D benchmarks, leaving adaptive error control for future work. We compute piecewise linear continuous discrete solutions in space and time by a general Galerkin (G2) finite element method (FEM). We introduce and compensate for mesh motion by defining a local arbitrary Euler–Lagrange (ALE) map on each space-time slab as part of the discretization, allowing a sharp phase interface given by θ on cell facets. The Unicorn implementation is published as part of the FEniCS Free Software system for automation of computational mathematical modeling. Simulation results are given for a 2D stationary convergence test, indicating quadratic convergence of the displacement, a simple 2D Poiseuille test for verifying correct treatment of the fluid-structure interface, showing quadratic convergence to the exact drag value, an established 2D dynamic flag benchmark test, showing a good match to published reference solutions and a 3D turbulent flag test as indicated above.

We present an adaptive finite element method for the compressible Euler equations, based on a posteriori error estimation of a quantity of interest in terms of a dual problem for the linearized equations. Continuous piecewise linear approximation is used in space and time, with componentwise weighted least-squares stabilization of convection terms and residual-based shock-capturing. The adaptive algorithm is demonstrated numerically for the quantity of interest being the drag force on a body.

In this paper. we identify and propose solutions for several issues encountered when designing a mesh adaptation package, such as mesh-to-mesh projections and mesh database design, and we describe an algorithm to integrate a mesh adaptation procedure in a physics solver. The open-source MAdLib package is presented as an example of such a mesh adaptation library. A new technique combining global node repositioning and mesh optimization in order to perform arbitrarily large deformations is also proposed. We then present several test cases to evaluate the performances of the proposed techniques and to show their applicability to fluid-structure interaction problems with arbitrarily large deformations. Copyright (C) 2009 John Wiley & Sons, Ltd.

General Galerkin (G2) is a new computational method for turbulent flow, where a stabilized Galerkin finite element method is used to compute approximate weak solutions to the Navier–Stokes equations directly, without any filtering of the equations as in a standard approach to turbulence simulation, such as large eddy simulation, and thus no Reynolds stresses are introduced, which need modelling. In this paper, G2 is used to compute the drag coefficient cD for the flow past a circular cylinder at Reynolds number Re=3900, for which the flow is turbulent. It is found that it is possible to approximate cD to an accuracy of a few percent, corresponding to the accuracy in experimental results for this problem, using less than 105 mesh points, which makes the simulations possible using a standard PC. The mesh is adaptively refined until a stopping criterion is reached with respect to the error in a chosen output of interest, which in this paper is cD. Both the stopping criterion and the mesh-refinement strategy are based on a posteriori error estimates, in the form of a space–time integral of residuals times derivatives of the solution of a dual problem, linearized at the approximate solution, and with data coupling to the output of interest. Copyright © 2008 John Wiley & Sons, Ltd.

We present analytical and computational evidence of blowup of initially smooth solutions of the incompressible Euler equations into non-smooth turbulent solutions. We detect blowup by observing inc ...

This article is a review of our work towards a parameter-free method for simulation of turbulent flow at high Reynolds numbers. In a series of papers we have developed a model for turbulent flow in the form of weak solutions of the Navier-Stokes equations, approximated by an adaptive finite element method, where: (i) viscous dissipation is assumed to be dominated by turbulent dissipation proportional to the residual of the equations, and (ii) skin friction at solid walls is assumed to be negligible compared to inertial effects. The result is a computational model without empirical data, where the only model parameter is the local size of the finite element mesh. Under adaptive refinement of the mesh based on a posteriori error estimation, output quantities of interest in the form of functionals of the finite element solution converge to become independent of the mesh resolution, and thus the resulting method has no adjustable parameters. No ad hoc design of the mesh is needed, instead the mesh is optimized based on solution features, in particular no boundary layer mesh is needed. We connect the computational method to the mathematical concept of a dissipative weak solution of the Euler equations, as a model of high Reynolds number turbulent flow, and we highlight a number of benchmark problems for which the method is validated. The purpose of the article is to present the computational framework in a concise form, to report on recent progress, and to discuss open problems that are subject to ongoing research.

Offshore floating platforms for wind turbines represent challenging concepts for design- ers trying to combine an optimal compromise between cost effectiveness and performance. Modelling of the hydrodynamic behaviour of the structure is still the subject of wide de- bate in the technical communities. The assessment of the hydrodynamics of the support structure is not an easy task as the floaters consist of an assembly of columns, braces and pontoons, commonly also with heave plates: each of these components corresponds to a different hydrodynamic model and it further interacts with the other elements. This results in very complex non-linear modeling, which makes it necessary to resort to computational fluid dynamics (CFD) methods for the evaluation of the combined hydrodynamics. In the framework of the collaboration between the Basque Centre for Applied Mathe- matics (BCAM) and Tecnalia R&I, the interaction of the sea flow with a semisubmersible floating offshore wind platform have been calculated by using the open source solver Uni- corn in the FEniCS-HPC framework when subject to a steady inflow. The prototype of the platform consists in a semi-submersible 4-columns column stabilized platform - NAUTILUS Floating Solutions concept-; columns are connected by a rigid ring pontoon provided with heave damping plates at the bottom. The novelty of the approach in FEniCS-HPC hinges upon an implicit formulation for the turbulence, a cheap free slip model of the boundary layer and goal-oriented mesh adaptivity [8, 6, 9, 20, 1]. We find that the results are consistent with experimental results for cylinders at high Reynolds number.

We present a new approach to Computational Fluid Dynamics CFD using adaptive stabilized Galerkin finite element methods with duality based a posteriori error control for chosen output quantities of interest. We address the basic question of computability in CFD: For a given flow, what quantity is computable to what tole rance to what cost? We focus on incompressible Newtonian flow with medium to large Rey nolds numbers involving both laminar and turbulent flow features. We estimate a poste riori the output of the computed solution with the output based on the exact solution to the Navier-Stokes equations, thus circumventing introducing and modeling Reynolds stresses in averaged Navier-Stokes equations. Our basic tool is a representation formula for th e error in the quantity of interest in terms of a space-time integral of the residual of a computed solution multiplied by weights related to derivatives of the solution of an associa ted dual problem with data connected to the output. We use the error representation formula to derive an a posterori error estimate combining residuals with computed dual weights, which is used for mesh adaptivity in space-time with the objective of satisfying a give n error tolerance with minimal computational effort. We show in concrete examples that outputs such as mean values in time of drag and lift of turbulent flow around a bluff body are c omputable on a PC with a tolerance of a few percent using a few hundred thousand mesh points in space. We refer to our methodology as Adaptive DNS/LES, where automatically by adaptivity certain features of the flow are resolved in a Direct Numerical Simulation DNS , while certain other small scale turbulent features are left unresolved in a Large Eddy Simulation LES . The stabilization of the Galerkin method giving a weighted least square control of the residual acts as the subgrid model in the LES. The a posteriori error estimate takes into account both the error from discretization and the error from the subgrid model. We pay particular attention to the stability of the dual solution from (i) perturbations repla cing the exact convection velocity by a computed velocity, and (ii) computational solution of the dual problem, which are the crucial aspects entering by avoiding using averaged Navier-Stokes equations including Reynolds stresses. A crucial observation is that the contri bution from subgrid modeling in the a posteriori error estimation is small, making it possib le to simulate aspects of turbulent flow without accurate modeling of Reynolds stresses.

In recent years adaptive stabilized finite element methods, here referred to as General Galerkin (G2) methods, have been developed as a general methodology for the computation of mean value output in turbulent flow. In earlier work, in the setting of bluff body flow, the use of no slip boundary conditions has been shown to accurately capture the separation from a laminar boundary layer, in a number of benchmark problems. In this paper we extend the G2 method to problems with turbulent boundary layers, by including a simple wall-model in the form of a friction boundary condition, to account for the skin friction of the unresolved turbulent boundary layer. In particular, we use G2 to simulate drag crisis for a circular cylinder, by adjusting the friction parameter to match experimental results. By letting the Reynolds number go to infinity and the skin friction go to zero, we get a G2 method for the Euler equations with slip boundary conditions, which we here refer to as an EG2 method. The only parameter in the EG2 method is the discretization parameter, and we present computational results indicating that EG2 may be used to model very high Reynolds numbers flow, such as geophysical flow.