Finite Element Discretization Research Articles

Variational data assimilation with finite-element discretization for second-order parabolic interface equation

Abstract In this paper, we propose and analyze a finite-element method of variational data assimilation for a second-order parabolic interface equation on a two-dimensional bounded domain. The Tikhonov regularization plays a key role in translating the data assimilation problem into an optimization problem. Then the existence, uniqueness and stability are analyzed for the solution of the optimization problem. We utilize the finite-element method for spatial discretization and backward Euler method for the temporal discretization. Then based on the Lagrange multiplier idea, we derive the optimality systems for both the continuous and the discrete data assimilation problems for the second-order parabolic interface equation. The convergence and the optimal error estimate are proved with the recovery of Galerkin orthogonality. Moreover, three iterative methods, which decouple the optimality system and significantly save computational cost, are developed to solve the discrete time evolution optimality system. Finally, numerical results are provided to validate the proposed method.

Stokes problem with slip boundary conditions using stabilized finite elements combined with Nitsche

We discuss how slip conditions for the Stokes equation can be handled using Nitsche method, for a stabilized finite element discretization. Emphasis is made on the interplay between stabilization and Nitsche terms. Well-posedness of the discrete problem and optimal convergence rates, in natural norm for the velocity and the pressure, are established, and illustrated with various numerical experiments. The proposed method fits naturally in the context of a finite element implementation while being accurate, and allows an increased flexibility in the choice of the finite element pairs.

Well‐posedness and error analysis of wave equations with Markovian switching

Compared to traditional partial differential equation modeling methods, Markov switching models can accurately capture the abrupt changes or jumps that complex systems often experience in the real world. In this paper, we propose a novel wave equation model with Markovian switching to represent complex systems with state‐jumping phenomena better, and the well‐posedness of the model is proved. In addition, a numerical method with non‐uniform grids is also proposed for the proposed model to simulate the data in realistic situations, which is based on the use of finite element discretization in space and central difference discretization in time. Finally, we conduct several experiments to analyze the errors and stability of the proposed model and the traditional model. The results show that the Markov switching model proposed in this paper has a smaller error than the traditional models while ensuring stability and can more accurately simulate the state jump phenomena of real‐world systems.

A Parallel Finite Element Discretization Algorithm Based on Grad-Div Stabilization for the Navier–Stokes Equations

A Parallel Finite Element Discretization Algorithm Based on Grad-Div Stabilization for the Navier–Stokes Equations

Multilevel Monte Carlo Methods for Stochastic Convection–Diffusion Eigenvalue Problems

We develop new multilevel Monte Carlo (MLMC) methods to estimate the expectation of the smallest eigenvalue of a stochastic convection–diffusion operator with random coefficients. The MLMC method is based on a sequence of finite element (FE) discretizations of the eigenvalue problem on a hierarchy of increasingly finer meshes. For the discretized, algebraic eigenproblems we use both the Rayleigh quotient (RQ) iteration and implicitly restarted Arnoldi (IRA), providing an analysis of the cost in each case. By studying the variance on each level and adapting classical FE error bounds to the stochastic setting, we are able to bound the total error of our MLMC estimator and provide a complexity analysis. As expected, the complexity bound for our MLMC estimator is superior to plain Monte Carlo. To improve the efficiency of the MLMC further, we exploit the hierarchy of meshes and use coarser approximations as starting values for the eigensolvers on finer ones. To improve the stability of the MLMC method for convection-dominated problems, we employ two additional strategies. First, we consider the streamline upwind Petrov–Galerkin formulation of the discrete eigenvalue problem, which allows us to start the MLMC method on coarser meshes than is possible with standard FEs. Second, we apply a homotopy method to add stability to the eigensolver for each sample. Finally, we present a multilevel quasi-Monte Carlo method that replaces Monte Carlo with a quasi-Monte Carlo (QMC) rule on each level. Due to the faster convergence of QMC, this improves the overall complexity. We provide detailed numerical results comparing our different strategies to demonstrate the practical feasibility of the MLMC method in different use cases. The results support our complexity analysis and further demonstrate the superiority over plain Monte Carlo in all cases.

A robust finite strain isogeometric solid-beam element

In this work, an efficient and robust isogeometric three-dimensional solid-beam finite element is developed for large deformations and finite rotations with merely displacements as degrees of freedom. The finite strain theory and hyperelastic constitutive models are considered and B-Spline and NURBS are employed for the finite element discretization. Similar to finite elements based on Lagrange polynomials, also NURBS-based formulations are affected by the non-physical phenomena of locking, which constrains the field variables and negatively impacts the solution accuracy and deteriorates convergence behavior. To avoid this problem within the context of a solid-beam formulation, the Assumed Natural Strain (ANS) method is applied to alleviate membrane and transversal shear locking and the Enhanced Assumed Strain (EAS) method against Poisson thickness locking. Furthermore, the Mixed Integration Point (MIP) method is used to make the formulation more efficient and robust. The proposed novel isogeometric solid-beam element is tested on several single-patch and multi-patch benchmark problems, and it is validated against classical solid finite elements and first-order Lagrange solid-beam element. The results show that the proposed formulation can alleviate the locking effects and significantly improve the performance of the isogeometric solid-beam element. With the developed element, efficient and accurate predictions of mechanical properties of lattice-based structured materials can be achieved. The proposed solid-beam element inherits both the merits of solid elements e.g. flexible boundary conditions and of the beam elements i.e. higher computational efficiency.

A MEEVC discretization for two-dimensional incompressible Navier-Stokes equations with general boundary conditions

In this work, we present a mass, energy, enstrophy and vorticity conserving (MEEVC) mixed finite element discretization for two-dimensional incompressible Navier-Stokes equations as an alternative to the original MEEVC scheme proposed in [A. Palha and M. Gerritsma, J. Comput. Phys., 2017]. The present method can incorporate general boundary conditions. Conservation properties are proven. Supportive numerical experiments with both exact and inexact quadrature are provided.

Influence of weak interlayer thickness on mechanical response and failure behavior of rock under true triaxial stress condition

In this study, we meticulously investigate the influence of the thickness of weak interlayer on mechanical response and failure behavior of rock mass near the free surface of tunnel, with the aid of a true triaxial test system. We simulated the true triaxial stress condition of rock mass containing weak interlayer proximate to the excavation boundary, where one face remains free, while others are subjected to loading stress and increasing tangential stress. The experimental results reveal a distinct change in the fracture morphology near the free surface of rock mass as the thickness of weak interlayer increases, concomitant with a gradual decrease in the fractal dimension. With the increase of the thickness of weak interlayer, the failure characteristics of rock mass change from a “binary” to a “unitary” feature, i.e. mixed tension-shear failure near the free surface and shear failure far away from the free surface changes to mixed tension-shear failure that only occurs near the free surface. Moreover, an inverse relationship is observed between the bearing capacity of the rock mass and the thickness of weak interlayer. Conversely, the peak strain in the tangential stress direction increases with the increase of the thickness of weak interlayer, illustrating that the slip shear strain of rock block positioned above the weak interlayer is gradually formed, and governs the strain behavior of rock mass in the tangential stress direction with the increases of the thickness of weak interlayer. By employing the finite-discrete element method (FDEM), we replicated the failure process of surrounding rocks containing weak interlayers near the free surface of tunnel, demonstrating considerable congruity with our physical test results. In light of the failure characteristics of surrounding rocks containing weak interlayers, we advocated a localized area priority reinforcement support strategy, prioritizing regions with weak interlayers. Numerical simulations verified the effectiveness of this approach in controlling fracture proliferation in surrounding rocks and inhibiting the splitting failure of surrounding rocks containing weak interlayers near the free surface of tunnel.

Modeling and simulation of pollution transport in the Mediterranean Sea using enriched finite element method

This paper presents a novel numerical method for simulating the transport and dispersion of pollutants in the Mediterranean sea. The governing mathematical equations consist of a barotropic ocean model with friction terms, bathymetric forces, Coriolis and wind stresses coupled to an advection–diffusion equation with anisotropic dispersion tensor and source terms. The proposed numerical solver uses a multilevel adaptive semi-Lagrangian finite element method that combines various techniques, including the modified method of characteristics, finite element discretization, coupled projection scheme based on a rotational pressure correction algorithm, and an adaptive L2-projection. The approach employs the gradient of the concentration as an error indicator for enrichment adaptations and increasing the number of quadrature points where needed without refining the mesh. The method is shown to provide accurate and efficient simulations for pollution transport in the Mediterranean sea. The proposed approach distinguishes itself from the well-established adaptive finite element methods for incompressible viscous flows by retaining the same structure and dimension of linear systems during the adaptation process.

An enhanced polar-based GPGPU-parallelized contact detection algorithm for 3D FDEM and its application to cracking analysis of shield tunnel segmental linings

This study introduces the three-dimensional combined finite-discrete element method (3D FDEM) to perform cracking analysis of segmental linings under extreme conditions. Considering the complex contact interactions of segments, this study first proposes a polar-based GPGPU-parallelized contact detection algorithm to handle memory issues confronted by existing FDEM algorithms. Initially, a spatial decomposition approach based on the polar coordinate system is implemented during the broad search phase. Each tetrahedral element is positioned within suitable search cells according to the axis-aligned bounding box. Subsequently, element pairs within each search cell are iteratively traversed, and potential contact pairs are identified using the judge cell criteria. After the broad search, a narrow search phase is executed to determine all real contacts. Following the above implementation, a load-structure model encompassing soil spring and external pressure calculations is then proposed. Based on the shape functions and coordinate transformation, the formulations of soil springs and external loads are derived and parallelized in the 3D FDEM framework. Three numerical tests are presented to validate the effectiveness of the proposed approach. Simulation results confirm the suitability of the proposed method for cracking analysis of segmental linings. Compared to the existing methods, it reduces GPU memory usage by 56 ∼ 76 % without extending time, which enhances the computational scale of 3D FDEM simulations. Furthermore, the proposed method is applied to two engineering scenarios, i.e., straight and curved segmental linings, both considering the absence and presence of contact defects between segments. The results show that contact defects can significantly reduce the resistance of the structural system comprising bolts and concrete segments for straight and curved scenarios.

Port-Hamiltonian discontinuous Galerkin finite element methods

Abstract A port-Hamiltonian (pH) system formulation is a geometrical notion used to formulate conservation laws for various physical systems. The distributed parameter port-Hamiltonian formulation models infinite dimensional Hamiltonian dynamical systems that have a nonzero energy flow through the boundaries. In this paper, we propose a novel framework for discontinuous Galerkin (DG) discretizations of pH-systems. Linking DG methods with pH-systems gives rise to compatible structure preserving semidiscrete finite element discretizations along with flexibility in terms of geometry and function spaces of the variables involved. Moreover, the port-Hamiltonian formulation makes boundary ports explicit, which makes the choice of structure and power preserving numerical fluxes easier. We state the Discontinuous Finite Element Stokes–Dirac structure with a power preserving coupling between elements, which provides the mathematical framework for a large class of pH discontinuous Galerkin discretizations. We also provide an a priori error analysis for the port-Hamiltonian discontinuous Galerkin Finite Element Method (pH-DGFEM). The port-Hamiltonian discontinuous Galerkin finite element method is demonstrated for the scalar wave equation showing optimal rates of convergence.

Differentiable solver for time-dependent deformation problems with contact

We introduce a general differentiable solver for time-dependent deformation problems with contact and friction. Our approach uses a finite element discretization with a high-order time integrator coupled with the recently proposed incremental potential contact method for handling contact and friction forces to solve ODE- and PDE-constrained optimization problems on scenes with complex geometry. It supports static and dynamic problems and differentiation with respect to all physical parameters involved in the physical problem description, which include shape, material parameters, friction parameters, and initial conditions. Our analytically derived adjoint formulation is efficient, with a small overhead (typically less than 10% for nonlinear problems) over the forward simulation, and shares many similarities with the forward problem, allowing the reuse of large parts of existing forward simulator code. We implement our approach on top of the open-source PolyFEM library and demonstrate the applicability of our solver to shape design, initial condition optimization, and material estimation on both simulated results and physical validations.

Tucker Tensor Approach for Accelerating Fock Exchange Computations in a Real-Space Finite-Element Discretization of Generalized Kohn-Sham Density Functional Theory.

The evaluation of Fock exchange is often the computationally most expensive part of hybrid functional density functional theory calculations in a systematically improvable, complete basis. In this work, we employ a Tucker tensor based approach that substantially accelerates the evaluation of the action of Fock exchange by transforming three-dimensional convolutional integrals into a tensor product of one-dimensional convolution integrals. Our numerical implementation uses a parallelization strategy that balances the memory and communication bottlenecks, alongside overlapping compute and communication operations to enhance computational efficiency and parallel scalability. The accuracy and computational efficiency are demonstrated on various systems, including Pt clusters of various sizes and a TiO2 cluster with 3684 electrons.

Numerical investigation of mechanical responses and failure features of coal measure strata using combined finite-discrete element method

Investigating the mechanical responses and failure features of the deep coal measure strata (CMS) is essential for designing hydraulic fracturing treatments and preventing geological dynamic hazards in underground engineering. In this paper, a three-dimensional CMS model is established by the combined finite-discrete element method (FDEM). The mechanical properties, stress evolutions and failure features of the CMS with different lithological combinations (sandstone-coal-sandstone (SCS), sandstone-coal-mudstone (SCM) and mudstone-coal-mudstone (MCM)) and thickness ratios are investigated under the uniaxial compressive condition. Results indicate that the SCS combination exhibits the highest uniaxial compressive strength (UCS) and elastic modulus, followed by the SCM and MCM combinations. The mechanical properties of the CMS are predominantly determined by the softer coal layer. The mechanical contrasts in elastic modulus and Poisson’s ratio between adjacent layers cause the generation of additional compressive stress in the softer layer and additional tensile stress in the stiffer layer. Therefore, the softer coal layer exhibits the shear-dominated failure pattern while the stiffer sandstone/mudstone layer displays the tensile-dominated failure pattern in the CMS. Furthermore, as the thickness ratio of the coal layer increases, the UCS, elastic modulus, and interfacial constraint stress decrease in the CMS. Additionally, the UCS and stress–strain relationship in the post-peak stage exhibit the noticeable size effect, whereas the mechanical behaviors in the pre-peak stage and the ultimate fracture patterns show the opposite. The elastic modulus and interfacial constraint stress obtained by the FDEM simulation match well with those calculated by theoretical models. The key findings of this paper can provide fundamental understanding and reliable parameters for multi-scale geomechanical modeling in petroleum and mining engineering.

Stabilized bases for high-order, interpolation semi-Lagrangian, element-based tracer transport

In a computational fluid model of the atmosphere, the advective transport of trace species, or tracers, can be computationally expensive. For efficiency, models often use semi-Lagrangian advection methods. High-order interpolation semi-Lagrangian (ISL) methods, in particular, can be extremely efficient, if the problem of property preservation specific to them can be addressed. Atmosphere models often use geometrically and logically nonuniform grids for efficiency and, as a result, element-based discretizations. Such grids and discretizations make stability a particular problem for ISL methods. Generally, high-order, element-based ISL methods that use the natural polynomial interpolant associated with a nodal finite-element discretization are unstable. We derive new bases having order of accuracy up to nine, with positive nodal weights, that stabilize the element-based ISL method. We use these bases to construct the linear advection operator in the property-preserving Interpolation Semi-Lagrangian Element-based Transport (Islet) method. Then we discuss key software implementation details. Finally, we show performance results for the Energy Exascale Earth System Model's atmosphere dynamical core, comparing the original and new transport methods. These simulations used up to 27,600 Graphical Processing Units (GPU) on the Oak Ridge Leadership Computing Facility's Summit supercomputer.

Estimating the macro strength of rock based on the determined mechanical properties of grains and grain-to-grain interfaces

Obtaining complete rock cores is exceptionally challenging in certain extreme environments, such as deep earth and deep space; consequently, it is difficult to obtain the macro strengths of rock, which are the key indexes for engineering design, via standard rock mechanical tests. This research indicates that by determining and leveraging the mechanical properties of grains and grain-to-grain interfaces, the rock macro strength can be effectively estimated. Nanoindentation tests were first conducted to determine the elastic modulus of the mineral grains and their interfaces. Subsequently, symmetrical and anti-symmetrical four-point bending tests were executed to establish the statistical law governing interfacial tensile and shear fracture toughnesses. Furthermore, by employing the Dugdale–Barenblatt model and considering the oblique angular distribution of interfacial cracks, the corresponding tensile strength, shear strength, and friction parameter were derived. These meso‑mechanical parameters were then inputted into the 3D Finite-Discrete Element Method to estimate rock macro strength. And the numerical results were then compared with the results of standard uniaxial compression, direct tension, Brazilian splitting, and direct shear tests. The consistency observed between the predictions and experimental results attests to the validity of the proposed method.

Efficient dynamic topology optimization of 2D metamaterials based on a complementary energy formulation

The advent of additive manufacturing has revolutionized the design and development of hierarchical structures, with potential applications in compliant, auxetic, and band-gap structures. This paper presents an innovative approach to developing a dynamic Topology Optimization (TO) framework for designing printable lattice structures that exhibit specific dynamic properties. Utilizing parametrically defined filament-based unit cell structures for topology optimization, we achieve desired natural frequency bandgaps in the structures composed of these unit cells. To enhance computational efficiency, we employ a complementary energy-based formulation to (semi)analytically derive the flexibility and stiffness matrices of the unit cell structure, thus, eliminating extensive finite element discretization. Consequently, a wide variety of parametrically defined filament-based meso-structures can be mathematically explored. We apply this innovative framework specifically for band-gap maximization of 2D lattice structures. By tuning the geometry within each cell using TO, we maximize the band gap. Our results show the potential of this approach to create more efficient and effective hierarchical structures with desired band-gap properties.

Advanced tensile fracture analysis of alumina ceramics: Integrating hybrid finite-discrete element modeling with experimental insights

This study offers novel insights into the Mode I tensile response of an alumina ceramic through the use of computational modeling and the flattened Brazilian disk (FBD) experiments. A modified hybrid finite-discrete element method (HFDEM) is developed, integrating a coupled damage and friction cohesive model and a microscopic stochastic fracture model with a Weibull strength distribution by Monte Carlo simulation. The model is used to simulate direct tensile failure processes under quasi-static loading conditions, providing qualitative and quantitative predictions of direct tensile failure processes of an alumina ceramic. Concurrently, quasi-static flattened Brazilian disk tests (indirect tensile tests) are performed on a standard MTS machine coupled with a high-speed camera. The modified HFDEM model is also applied to reproduce the FBD experiments, and our simulated tensile strength is consistent with the experimental results. The results of the modified HFDEM model show three kinds of phenomena (i.e., “underestimation”, “reasonable estimation”, and “overestimation” of the indirect tensile strength) and four different types of associated fracture and fragment patterns of FBD testing. The integration of simulation and experimental results reveal relationships between fracture patterns, fragment geometry, tensile strength, and indirect tensile strength. The fracture and fragmentation patterns derived from our modified HFDEM model can be utilized to analyze the “tensile strength” measured in BD testing. Overall, this research offers important insights and direction for future Brazilian disk experiments and tensile strength assessments, enhancing our comprehension of ceramic failure mechanisms under tensile loadings.

Finite Element Discretizations of a Convective Brinkman–Forchheimer Model Under Singular Forcing

Finite Element Discretizations of a Convective Brinkman–Forchheimer Model Under Singular Forcing

Effects of Geo-Stress and Natural Fracture Strength on Hydraulic Fracture Propagation in a Deep Fractured Tight Sandstone Reservoir.

Hydraulic fracturing technology has been widely used in the tight reservoir reconstruction. Unfortunately, with the deepening of mining depth and the increase of geo-stress, the propagation mechanism of medium-pressure fractures in the reservoir is significantly different from that of conventional shallow reservoirs. Based on the combined finite discrete element method, this paper conducts numerical simulation research on deep tight sandstone reservoirs in the west. The discrete fracture network modeling method is used to establish a tight sandstone reservoir model with natural bedding, and the influence of geo-stress difference and natural fracture strength on hydraulic fracture propagation law in a high geo-stress environment is discussed in detail. The results show that the difference between geo-stress and the strength of natural fractures has a significant effect on the shape and expansion of hydraulic fractures under the high geo-stress conditions. The greater the difference in ground stress, the more obvious the tendency of the main fractures of the reservoir, and the shorter the branch fractures. With the increase of natural fracture strength, the changes in propagation pressure, fracture length, area, and width, which can be fitted with a linear function with a goodness of fit as high as 0.99. In addition, the morphological results of hydraulic fractures in the simulation are not only affected by the constitutive parameters of the model but also may be affected by the randomness of the natural fracture network, thus, showing a certain degree of dispersion. Therefore, it is extremely necessary to build a reservoir fracturing model in a specific area based on more detailed geological monitoring data to guide actual construction. The above achievements have certain reference significance for the field operation of deep tight sandstone reservoirs.