Mesh Refinement Research Articles

Finite element solution of coupled multiphysics reaction-diffusion equations for fracture healing in hard biological tissues

This study proposes a computational framework to investigate the multi-stage process of fracture healing in hard tissues, e.g., long bone, based on the mathematical Bailon-Plaza and Van der Meulen formulation. The goal is to explore the influence of critical biological factors by employing the finite element method for more realistic configurations. The model integrates a set of variables, including cell densities, growth factors, and extracellular matrix contents, managed by a coupled system of partial differential equations. A weak finite element formulation is introduced to enhance the numerical robustness for coarser mesh grids, complex geometries, and more accurate boundary conditions. This formulation is less sensitive to mesh quality and converges smoothly with mesh refinement, exhibiting superior numerical stability compared to previously available strong-form solutions. The model accurately reproduces various stages of healing, including soft cartilage callus formation, endochondral and intramembranous ossification, and hard bony callus development for various sizes of fracture gap. Model predictions align with the existing research and are logically coherent with the available experimental data. The developed multiphysics simulation clarifies the coordination of cellular dynamics, extracellular matrix alterations, and signaling growth factors during fracture healing.

Application of polyhedral meshes in large-eddy simulation of building array flow fields within neutral and unstable boundary layers

Urbanization has intensified the urban heat island phenomenon, a critical concern for human habitation in recent decades. Analyzing flow fields around building arrays is vital for mitigating issues such as heatstroke and air pollution, thereby enhancing human comfort in urban environments. However, challenges in accuracy and computational efficiency persist in simulating urban thermal environments. This study evaluates how boundary layer meshes in polyhedral meshes affect wind and temperature distributions in neutral and unstable boundary layers, utilizing large-eddy simulation (LES). A total of ten cases were examined, including eight with local mesh refinement and two without. The study demonstrates that local mesh refinement is an effective strategy to reduce mesh count by approximately 60 %, thereby reducing computational time by 60 % using the same computational configuration. The performance of the polyhedral meshes was found to be comparable to that of the hexahedral meshes. For the mean velocity and temperature, the discrepancies between the hexahedral and polyhedral cases were within 10 %, while for the second-order statistics, the discrepancies were within 20 %. Within the neutral boundary layer, boundary layer meshes on the ground and building surfaces had minimal impact on the wind field. However, these meshes enhanced the accuracy of low percentile wind speed estimations by approximately 20 % within the unstable boundary layer. Moreover, this study lays the foundation for research on the thermal environment in actual urban areas and contributes to the guidelines for LES with non-isothermal inflow conditions.

High-fidelity computational study of aerodynamic noise of side-by-side rotor in full configuration

This paper presents a numerical study of side-by-side rotor aircraft in full configuration during hover using high-fidelity computational fluid dynamics (CFD) and aeroacoustic simulations. CFD simulations are performed using HPCMP CREATE™-AV Helios, while acoustic calculations are conducted with PSU-WOPWOP. Overset structured grids are applied using the SA-DES turbulence model along with adaptive mesh refinement in the near- and off-body meshes. Four overlap cases at a specific collective pitch angle of 8° are considered. Results reveal a substantial influence of the fuselage on rotor performance, with improvement in the figure of merit being observed for each overlap configuration due to a partial ground effect. There is considerable increase in blade sectional thrust, pressure fluctuations, and noise levels seen in the 0% overlap configuration compared to the 25% overlap configuration. This increase is primarily due to the absence of rotor overlap, resulting in a lower induced velocity allowing for flow to travel upward, and reaching the rotor disk plane. The acoustic spectrum appears at all harmonics of the blade passing frequency for the 0% and 5% overlap configurations. In contrast, the acoustic spectrum is distinctly present only at even harmonics for the 15% and 25% overlap cases. A noise reduction benefit of 3-5 dBA is noticed in the full configuration cases having rotor overlap compared to the 0% overlap configuration. Acoustic destructive interference is observed in the directivity plot only for the 15% and 25% overlap configurations. Overall, this paper reveals that the fuselage has a significant impact on the performance, aerodynamics, and aeroacoustics of the side-by-side UAM vehicle in hover.

A block-spectral adaptive H-/p-refinement strategy for shock-dominated problems

An adaptive H-/p-refinement strategy using a novel sensor is devised and tested in a block-spectral compressible Euler code equipped with adaptive-mesh refinement (AMR) and high-order flux-reconstruction numerics. At each Gauss quadrature point (or solution point) within each spectral block (or mesh element) the discrete velocity jump ΔU = ∂U/∂y1Δy1+∂V/∂y2Δy2+∂W/∂y3Δy3 is calculated and normalized by the local speed of sound, a. The grid spacing, Δxi, is calculated in each direction as the distance between auxiliary Gauss-Lobatto points, staggered relative to the solution points. The polynomial order is increased from p=0 to p=pmax in regions of weak compression, (ΔU/a)crit<ΔU/a<0 and kept at p=pmax in regions of flow expansion ΔU/a≥0, while staying at the H=0 base mesh level. Regions experiencing strong compressions, i.e. ΔU/a<(ΔU/a)crit, are H-refined up to H=Hmax where Hmax is applied at the location of maximum compression, ΔU/a=min(ΔU/a) in the domain, while keeping p=0 to guarantee robustness and monotonicity of the solution in the H refined region. The critical value of (ΔU/a)crit= -0.06 is found to effectively separate smooth and non-smooth solution regions, supported by a 1D detonation initiation test case in ideal gas and a shock-to-detonation transition in high explosives. Using this value, the Sod shock tube, Shu-Osher problem, double Mach reflection and a 2D detonation in a high-explosive are simulated with the proposed adaptive H-/p-refinement. In the Sod shock tube case, p-refinement resolves the (weak) contact discontinuity while H-refinement enhances the grid resolution in the shock exploiting the monotonicity of the p=0 reconstruction. For the Shu-Osher problem, p-refinement captures the small-scale oscillations trailing the shock that would be otherwise attenuated, while H-refinement triggered by the ΔU-sensor appropriately tracks the shock. In the double Mach reflection problem, H-refinement confines the numerical diffusion around the reflected shock while p-refinement recaptures many physical features trailing the shock. Finally, in the 2D high-explosive detonation case, H-refinement follows the leading shock and resolves the curvature of the detonation wave, while p-refinement adds resolution to the trailing reaction zone. Finally, the proposed methodology is tested in a detonation-wave propagation test case in high-explosives with numerical predictions comparing favorably against experiments.

Vortex-p: A Helmholtz-Hodge and Reynolds decomposition algorithm for particle-based simulations

Astrophysical turbulent flows display an intrinsically multi-scale nature, making their numerical simulation and the subsequent analyses of simulated data a complex problem. In particular, two fundamental steps in the study of turbulent velocity fields are the Helmholtz-Hodge decomposition (compressive+solenoidal; HHD) and the Reynolds decomposition (bulk+turbulent; RD). These problems are relatively simple to perform numerically for uniformly-sampled data, such as the one emerging from Eulerian, fix-grid simulations; but their computation is remarkably more complex in the case of non-uniformly sampled data, such as the one stemming from particle-based or meshless simulations. In this paper, we describe, implement and test vortex-p, a publicly available tool evolved from the vortex code, to perform both these decompositions upon the velocity fields of particle-based simulations, either from smoothed particle hydrodynamics (SPH), moving-mesh or meshless codes. The algorithm relies on the creation of an ad-hoc adaptive mesh refinement (AMR) set of grids, on which the input velocity field is represented. HHD is then addressed by means of elliptic solvers, while for the RD we adapt an iterative, multi-scale filter. We perform a series of idealised tests to assess the accuracy, convergence and scaling of the code. Finally, we present some applications of the code to various SPH and meshless finite-mass (MFM) simulations of galaxy clusters performed with OpenGadget3, with different resolutions and physics, to showcase the capabilities of the code.

Non-uniform knot (NUK) SIAC post-processing of flow fields produced through unstructured grid adaptation and optimization

As the finite element method (FEM) and the finite volume method (FVM), both their traditional and high-order variants, continue their proliferation into various applied engineering disciplines, adaptive mesh refinement and optimization strategies have increased in their importance when solving real-world computational fluid mechanics applications. The post-processing and visualization of the resulting flow fields present two significant analysis and visualization challenges. The first challenge is the handling of elemental continuity, which is often only C0 continuous (in continuous Galerkin methods) or piecewise discontinuous (in discontinuous Galerkin methods). The second challenge is that, depending on the flow regime and the geometric configurations for which adaptive meshing strategies are used, the meshes generated are often highly anisotropic. The (uniform knot) line-SIAC (L-SIAC) filter has been proposed as a way of handling elemental continuity issues in an accuracy-conserving manner with the added benefit of casting the data in a smooth context even if the representation is element discontinuous. In this paper, we demonstrate that the state-of-the-art adaptive L-SIAC filter, designed for mildly anisotropic meshes, suffers degradation in the quality of the post-processed solution when applied to the types of highly anisotropic meshes produced through adaptive mesh refinement and optimization. Hence, a new Non-Uniform Knot (NUK) L-SIAC filter is proposed that automatically conforms to the underlying mesh anisotropy. We demonstrate that the new filter behaves similarly to the adaptive L-SIAC filter when applied to uniform and mildly anisotropic meshes, and furthermore we show the superiority of the NUK L-SIAC filter when applied to highly anisotropic meshes. The newly formulated filter is applied to 2D canonical scalar fields and used to visualize 2D and 3D fluid flow simulation results.

Guidelines for element size and type selection for the finite element simulation of laser-induced elastic waves in thermoelastic laser ultrasonic testing

This paper explores spatial discretization within finite element simulations of laser-induced elastic waves within the context of Laser Ultrasonic Testing (LUT). Motivated by discrepancies and oscillations detected in temperature and displacement results in the literature, we traced these issues back to spatial discretization challenges. These challenges originate from rapid localized heating and the generation and propagation of high-frequency waves across a relatively large domain. This study effectively addresses and rectifies these inaccuracies, offering guidance for selecting the appropriate element size and type. We examined two element types: four-node quadrilaterals (Q4) employing first-order Lagrange and nine-node quadrilaterals (Q9) using second-order Lagrange shape functions. Our analysis encompasses mesh refinement strategies, exploration of time and frequency domain plots for temperature and displacement, as well as an evaluation of different pulse durations. Our findings demonstrate that Q9 elements attain accuracy with grids four times larger than Q4 elements for temperature and wave propagation analyses. Furthermore, we observe that lower frequency waves exhibit reduced sensitivity to element size, emphasizing the relationship between element size and elastic wave frequency. Pulse durations in the 6 to 30 ns range affect the required element size in the heat-affected zone but exert minimal influence on wave frequency and spatial discretization in the remainder of the domain. Finally, we present a new formula for element size selection based on the dominant frequency. This study provides a comprehensive guideline for selecting element size and type, enabling the attainment of accurate results while effectively managing computational costs.

Effects of reservoir mechanical properties on induced seismicity during subsurface hydrogen storage.

The intermittent storage of hydrogen in subsurface porous media such as depleted gas fields could be pivotal to a successful energy transition. Numerical simulations investigate the intermittent storage of hydrogen in a porous, depleted subsurface reservoir. Various parametric studies are performed to assess the effect of mechanical properties of the reservoir (i.e. Young's modulus, Poisson's ratio, Biot coefficient and permeability) on the induced fault slip of a single through-going fault that transverses the entire reservoir. Simulations are run using a three-dimensional, finite element, fully coupled hydromechanical code with explicit representations of layers and faults. The effect of the domain mesh refinement and fault mesh refinement on the fault slip versus operation time solution is investigated. The fault is observed to slip in two distinct events, one during the second injection period and one in the third injection period. The fault is not observed to slip during the storage or withdrawal periods. It is found that in order to minimize seismic risk, a reservoir rock with high Young's modulus (>40 GPa), high Poisson's ratio (>0.30) and high Biot coefficient (>0.65) would be preferable for hydrogen storage. Reservoir rocks of low Young's modulus (10-30 GPa), intermediate Poisson's ratio (0.00-0.30) and low-to-intermediate Biot coefficient (0.25-0.65), at high injection rates, were found to have higher potential of inducing large seismic events.This article is part of the theme issue 'Induced seismicity in coupled subsurface systems'.

Formation dynamics of the satellite droplet in the breakup of a symmetrical liquid bridge

Inkjet printing technology has played an irreplaceable role in life science, precision manufacturing, and other frontier fields in recent years. However, the further development of this technology is limited by the fact that its printing resolution is difficult to raise to a higher level. The emerging satellite droplet printing technology offers a new approach for inkjet printing to break through the bottleneck of printing resolution limitations. In this paper, a symmetrical satellite droplet printing strategy is proposed. The effects of the geometric parameters of the satellite droplet generating device, the physical properties of the ink, and the operating parameters on the liquid bridge breakup process and the size of the satellite droplet are systematically studied. The phase field method and adaptive mesh refinement strategy are applied to solve the two-dimensional symmetrical model. The results indicate that the length of the liquid bridge, the radius of the bridge, the viscosity of the ink, and the drainage velocity are all positively correlated with the satellite droplet size, while the surface tension coefficient has a negative correlation with the satellite droplet size. Furthermore, the three-phase contact line at the orifice end will slip toward the center if the initial radius of the liquid bridge is quite large. Based on these investigations and discussions, a corresponding effective working space for satellite droplet printing is obtained, which lays the foundation for the popularization and further development of satellite droplet printing technology.

Numerical investigation of the interaction between an interface and a decaying Lamb–Oseen vortex

The present study investigates the dynamics of the interface in the presence of a decaying Lamb–Oseen vortex, and four distinct wave patterns are observed: non-breaking waves with small periodic oscillations, plunging breakers, depression breakers, and gravity–capillary waves. The deformation of the interface is induced by a two-dimensional Lamb–Oseen vortex, and the study examines the influence of vortex strength and surface tension on the resulting flow. The wave dynamics are characterized as a function of the Reynolds and Weber numbers, and a phase diagram is presented in terms of (Re, We) to distinguish the different wave patterns. To ensure accurate reconstruction of the interface, the numerical methods used in this study feature a mass and momentum consistent advection method, high-order interpolation schemes, and a block-structured adaptive mesh refinement strategy. The study presents the characteristics of the air cavity entrained at the moment of wave impact for each wave pattern. Furthermore, the results provide insight into the nature of bubble entrainment within a vortex and reveals the bubble entrainment process via a breakup cascade. Meanwhile, it is also shown that the entrainment of bubble results in significant vortex distortion. Overall, this research contributes to enhance our understanding of wave dynamics and the intricate interaction between vortices and interfaces.

Multirate time-integration based on dynamic ODE partitioning through adaptively refined meshes for compressible fluid dynamics

In this paper, we apply the Paired-Explicit Runge-Kutta (P-ERK) schemes by Vermeire et al. [1,2] to dynamically partitioned systems arising from adaptive mesh refinement. The P-ERK schemes enable multirate time-integration with no changes in the spatial discretization methodology, making them readily implementable in existing codes that employ a method-of-lines approach.We show that speedup compared to a range of state of the art Runge-Kutta methods can be realized, despite additional overhead due to the dynamic re-assignment of flagging variables and restricting nonlinear stability properties. The effectiveness of the approach is demonstrated for a range of simulation setups for viscous and inviscid convection-dominated compressible flows for which we provide a reproducibility repository.In addition, we perform a thorough investigation of the nonlinear stability properties of the Paired-Explicit Runge-Kutta schemes regarding limitations due to the violation of monotonicity properties of the underlying spatial discretization. Furthermore, we present a novel approach for estimating the relevant eigenvalues of large Jacobians required for the optimization of stability polynomials.

An adaptive dynamic phase-field modeling with variable-node elements for thermoelastic fracture in orthotropic media

In this work, a novel adaptive algorithm integrated into a dynamic phase-field method is developed to model the fracture of orthotropic materials under thermoelastic loading. The procedure of local refinement is synchronized with crack tip advancement, with the phase-field value serving as the criterion for marking elements. To address mesh refinement inconsistencies, variable-node elements are employed. Hughes–Hilbert–Taylor (HHT) and backward difference method are utilized for time discretization of displacement and temperature, respectively. A penalized structural matrix is incorporated into crack surface density to realize the orthotropy of crack growth. The performance of the proposed method is validated through several numerical benchmarks, i.e., the proposed method can effectively simulate the dynamic crack propagation of orthotropic materials under thermoelastic loading, and can assuage computational overhead while keeping acceptable accuracy.

Incorporating Airspace Constraints in Multi-phase 4D Trajectory Optimization Based on Improved Adaptive Mesh Refinement Method

Incorporating Airspace Constraints in Multi-phase 4D Trajectory Optimization Based on Improved Adaptive Mesh Refinement Method

Optimal Battery Energy Storage Dispatch for the Day-Ahead Electricity Market

This work presents an innovative application of optimal control theory to the strategic scheduling of battery storage in the day-ahead electricity market, focusing on enhancing profitability while factoring in battery degradation. This study incorporates the effects of battery degradation on the dynamics in the optimisation framework. Considering this cost in economic analysis and operational strategies is essential to optimise long-term performance and economic viability. Neglecting degradation costs can lead to suboptimal operation and dispatch strategies. We employ a continuous-time representation of the dynamics, in contrast with many other studies that use a discrete-time approximation with rather coarse intervals. We adopt an equivalent circuit model coupled with empirical degradation parameters to simulate a battery cell’s behaviour and degradation mechanisms with good support from experimental data. Utilising direct collocation methods with mesh refinement allows for precise numerical solutions to the complex, nonlinear dynamics involved. Through a detailed case study of Belgium’s day-ahead electricity market, we determine the optimal charging and discharging schedules under varying objectives: maximising net revenues, maximising profits considering capacity degradation, and maximising profits considering both capacity degradation and internal resistance increase due to degradation. The results demonstrate the viability of our approach and underscore the significance of integrating degradation costs into the market strategy for battery operators, alongside its effects on the battery’s dynamic behaviour. Our methodology extends previous work by offering a more comprehensive model that empirically captures the intricacies of battery degradation, including a fine and adaptive time domain representation, focusing on the day-ahead market, and utilising accurate direct methods for optimal control. This paper concludes with insights into the potential of optimal control applications in energy markets and suggestions for future research avenues.

Isogeometric analysis of functionally graded panels using Bézier triangles

Isogeometric Analysis is a numerical method that integrates the concepts of geometric modeling and structural analysis. It approximates the displacement field using the same basis functions employed by CAD systems to describe the structure’s geometry. This work proposes an isogeometric formulation for analysis of functionally graded panels based on rational Bézier triangles, allowing the exact geometry representation and automatic discretization of topologically complex models. The formulation is applied to the free vibration and stability analysis of functionally graded plates and curved panels. Monotonic convergence under mesh refinement was observed in all examples. Furthermore, results show that curved functionally graded panels display a complex nonlinear behavior and can present bifurcation buckling before reaching the limit load.

Mother-leaf-method accelerated parallel-GPU AMR phase-field simulations of dendrite growth

Parallel computing with multiple graphics processing units (GPUs) for adaptive mesh refinement (AMR) is a powerful method that improves the computational efficiency of phase-field (PF) simulations of dendrite growth. Octree-type block-structured AMR is often used in parallel computation because of its good parallelism; however, the process of filling the buffer region with physical data from adjacent blocks constitutes a major overhead cost that reduces computational performance. In this study, the mother-leaf (ML) method, which can reduce this overhead, was applied to accelerate parallel-GPU AMR PF simulations of dendrite growth. To evaluate the computational performance of the ML method, single- and parallel-GPU AMR PF simulations of dendrite growth were performed. The results show that single- and parallel-GPU AMR computations with the ML method accelerate the PF simulations of dendrite growth approximately 1.2–1.6 times and 1.2–1.4 times, respectively, compared to the conventional method without the ML method. In large-scale dendrite-growth PF simulations, such accelerations are significant for investigating various dendrite growth phenomena.

Adaptive meshing strategies for nanophotonics using a posteriori error estimation.

As nanophotonic devices become increasingly complex, computer simulations of such devices are becoming ever more important. Unfortunately, computer simulations of nanophotonic devices are computationally expensive, especially if many simulations are necessary, e.g., when optimizing or inverse designing a device. Here we study adaptive mesh refinement for finite-element method simulations using an a posteriori error estimation method. We demonstrate that the use of adaptive meshing leads to faster convergence with lower memory footprint for complex three-dimensional nanophotonic structures. Nevertheless, one needs to be careful to avoid a mesh propagation effect for adaptive mesh refinement to be a successful strategy.

Multiscale graph neural networks with adaptive mesh refinement for accelerating mesh-based simulations

Mesh-based Graph Neural Networks (GNNs) have recently shown capabilities to simulate complex multiphysics problems with accelerated performance times. However, mesh-based GNNs require a large number of message-passing (MP) steps and suffer from over-smoothing for problems involving very fine mesh. In this work, we develop a multiscale mesh-based GNN framework mimicking a conventional iterative multigrid solver, coupled with adaptive mesh refinement (AMR), to mitigate challenges with conventional mesh-based GNNs. We use the framework to accelerate phase field (PF) fracture problems involving coupled partial differential equations with a near-singular operator due to near-zero modulus inside the crack. We define the initial graph representation using all mesh resolution levels. We perform a series of downsampling steps using Transformer MP GNNs to reach the coarsest graph followed by upsampling steps to reach the original graph. We use skip connectors from the generated embedding during coarsening to prevent over-smoothing. We use Transfer Learning (TL) to significantly reduce the size of training datasets needed to simulate different crack configurations and loading conditions. The trained framework showed accelerated simulation times, while maintaining high accuracy for all cases compared to physics-based PF fracture model. Finally, this work provides a new approach to accelerate a variety of mesh-based engineering multiphysics problems.

Properties of molecular clumps and cores in colliding magnetized flows

ABSTRACT We simulate the formation of molecular clouds in colliding flows of warm neutral medium with the adaptive mesh refinement code flash in eight simulations with varying initial magnetic field strength, between 0.01–5 μG. We include a chemical network to treat heating and cooling and to follow the formation of molecular gas. The initial magnetic field strength influences the fragmentation of the forming cloud because it prohibits motions perpendicular to the field direction and hence impacts the formation of large-scale filamentary structures. Molecular clump and core formation occurs anyhow. We identify 3D clumps and 3D cores, which are defined as connected, CO-rich regions. Additionally, 3D cores are heavily shielded. While we do not claim those 3D objects to be directly comparable to observations, this enables us to analyse their full virial state. With increasing field strength, we find more fragments with a smaller average mass; yet the dynamics of the forming clumps and cores only weakly depends on the initial magnetic field strength. The molecular clumps are mostly unbound, probably transient objects, which are weakly confined by ram pressure or thermal pressure, indicating that they are swept up by the turbulent flow. They experience significant fluctuations in the mass flux through their surface, such that the Eulerian reference frame shows a dominant time-dependent term due to their indistinct nature. We define the cores to encompass highly shielded molecular gas. Most cores are in gravitational-kinetic equipartition and are well described by the common virial parameter $\alpha _\mathrm{vir}$, while some undergo minor dispersion by kinetic surface effects.

Multi-spacecraft study with the Icarus model. Modelling the propagation of CMEs to Mercury and Earth

Coronal mass ejections (CMEs) are the main drivers of the disturbances in interplanetary space. Earth-directed CMEs can be dangerous, and understanding the CME interior magnetic structure is crucial for advancing space weather studies. It is important to assess the capabilities of a numerical heliospheric model, as a firm understanding of the nature and extent of its limitations can be used to improve the model and the space weather predictions based on it. The aim of the present study is to test the capabilities of the recently developed heliospheric model Icarus and the linear force-free spheromak model that has been implemented in it. To validate the Icarus space weather modelling tool, two CME events were selected that were observed by two spacecraft located near Mercury and Earth, respectively. This enables us to test the heliospheric model computed with Icarus at two distant locations. The source regions for the CMEs were identified, and the CME parameters were determined and later optimised. Different adaptive mesh refinement levels were applied in the simulations to assess its performance by comparing the simulation results to in situ measurements. The first CME event erupted at 15:25 on July 9, 2013. The modelled time series were in good agreement with the observations both at MESSENGER and ACE. The second CME event started at 10:25 on February 16, 2014, and was more complicated, as three CME interactions occurred in this event. It was impossible to recover the observed profiles without modelling the other two CMEs that were observed, one before the main CME and one afterward. The parameters for the three CMEs were identified and the three CMEs were modelled in Icarus. For both CME studies, AMR level 3 was sufficient to reconstruct small-scale features near Mercury, while at Earth, AMR level 4 was necessary due to the radially stretched grid that was used. The profiles obtained at both spacecraft resemble the in situ measurements well. The current limitations of the space weather modelling tool result in an excessively small deceleration of the CME propagation during the CME--CME interaction as measured by MESSENGER and ACE.