Time Step Refinement Research Articles

An implicit large-eddy simulation perspective on the flow over periodic hills

The periodic hills simulation case is a well-established benchmark for computational fluid dynamics solvers due to its complex features derived from the separation of a turbulent flow from a curved surface. We study the case with the open-source implicit large-eddy simulation (ILES) software Lethe. Lethe solves the incompressible Navier–Stokes equations by applying a stabilized continuous finite element discretization. The results are validated by comparison to experimental and computational data available in the literature for Re = 5600. We study the effect of the time step, averaging time, and global mesh refinement. The ILES approach shows good accuracy for average velocities and Reynolds stresses using less degrees of freedom than the reference numerical solution. The time step has a greater effect on the accuracy when using coarser meshes, while for fine meshes the results are rapidly time-step independent when using an implicit time-stepping approach. A good prediction of the reattachment point is obtained with several meshes and this value approaches the experimental benchmark value as the mesh is refined. We also run simulations at Reynolds equal to 10600 and 37000 and observe promising results for the ILES approach.

Fluid–structure interaction simulation of a flapping flag in a laminar jet

The flapping flag is a classical and important Fluid–Structure Interaction (FSI) problem with many applications in nature, industry, and biological systems. However, the interaction between the flow field and the flag response is not fully understood due to the complex and nonlinear nature of the problem. Although theoretical models qualitatively predict the flapping frequency and the onset velocity, there are still quantitative discrepancies between the theoretical predictions and the experimental data. In this study, numerical FSI simulations of a standard flag in a viscous uniform laminar flow are performed to characterize the flow field around the flapping flag and correlate the flow field with the flag response for six Reynolds numbers (Re= 4.5 ×104, Re= 5.3 ×104, Re= 5.7 ×104, Re= 5.9 ×104, Re= 6.6 ×104, and Re= 7.5 ×104). Two-way coupled FSI simulations are conducted using commercial software ANSYS®. Rigorous solution verification is conducted first on three Re to ensure that solutions are independent of mesh and time-step refinement. Then validation is achieved by comparing the simulation predictions with experimental data. Overall, results showed a good agreement with an error range between 0.7%–4.7% for oscillation amplitude (Am), 1.3%–6.41% for drag coefficient (Cd), and 1.3%–5.3% for frequency (f). Moreover, a comparison between inviscid and viscous models is performed to evaluate the accuracy of using inviscid theoretical models for the same flow. It is observed that the inviscid model accurately predicted f but underpredicted Am and Cd, and overpredicted lift coefficient Cl. In the case of the viscous model, it is observed that a flow separation bubble and a thin flow circulation region exist on the suction side of the flag during the limit cycle oscillation (LCO). Finally, an instantaneous negative drag is observed for a short time during the LCO, suggesting that flag reconfiguration and flow separation significantly affects the drag. These observations suggest that viscous models may be needed to improve the theoretical predictions.

CFD-DEM simulation of fluorination reaction in fluidized beds with local grid and time refinement method

The gas-solid reaction process with wide particle size distribution is extensively used in the chemical engineering field, especially the particle reacts with the gas gradually, such as fluorination reactions in fluidized beds. When the computational fluid dynamics-discrete element method (CFD-DEM) is used for the coupling simulation of multiphase and polydisperse particle reaction system, the grid size directly affects the accuracy of flow field information and simulation of chemical reaction. Furthermore, particle calculation time step will directly affect the efficiency of coupling calculation. In this work, a local grid and time step refinement method is proposed to simulate multiphase and polydisperse particle fluidization reaction system. In this method, the refined DEM grids are automatically generated in the computational domain around the fine particles, and the detailed fluid phase information is obtained with the interpolation algorithm. In the two-phase coupling process, particles are divided into different groups based on physical properties, each group has its own independent time step. The multistage conical-cylindrical spouted bed is proposed for the fluorination reaction process; the operating gas velocity range of the polydisperse particle system is extended by the new design while the particle size distribution changes with the gas-solid reaction process. It is demonstrated that the local grid and time step refinement method can improve the accuracy and efficiency of the traditional CFD-DEM method in the reaction process simulation, which describes a polydisperse particle system with wide particle size distribution. Aimed at improving the simulation accuracy and efficiency, this paper will be helpful for simulating the particle reaction process in the gas-solid fluidized bed and beneficial for the development of the CFD-DEM method.

Adaptive Finite Element Simulation of Double-Diffusive Convection

Double-diffusive convection plays an important role in many physical phenomena of practical importance. However, the numerical simulation of these phenomena is challenging since fine meshes are often required to capture the flow physics. Hence, several different numerical methods have been employed in the past. This work reports the development and application of an adaptive finite element method for the simulation of these phenomena, thereby avoiding the need for the use of very fine meshes over the whole domain. The weak formulation of the conservation equations for mass, momentum, energy and species concentration is used. The Boussinesq approximation relates the density of the fluid to the temperature and/or the species concentration. A second-order backward difference method is used for time discretization and the Galerkin method is employed for spatial discretization. Both adaptive time step and grid refinement techniques are employed, and the code is parallelized using MPI. Three different stabilization methods of the convective-diffusion equations are compared; namely, the streamline upwind Petrov–Galerkin (SUPG) method, and two modified methods aimed at diminishing spurious oscillations that include an artificial diffusion term. This diffusion term may be either isotropic or orthogonal to the streamlines. The addition of artificial isotropic diffusion to the SUPG method provides enhanced stability. The method is applied to double-diffusive finger convection in a sucrose-salt aqueous mixture and a stratified salt solution heated from below. The method accurately reproduces the experimentally observed temporal evolution of the salt fingers in the former case and the location of the interfaces between convective and non-convective zones in the latter.

A general framework to derive linear, decoupled and energy-stable schemes for reversible-irreversible thermodynamically consistent models

This paper presents a general numerical platform for designing accurate, efficient, and stable numerical algorithms for incompressible hydrodynamic models that obey thermodynamical laws. The obtained numerical schemes are automatically linear in time. It decouples the hydrodynamic variable and other state variables such that only small-size linear problems need to be solved at each time marching step. Furthermore, if the classical velocity projection method is utilized, the velocity field and pressure field can be decoupled. In the end, only a few elliptic-type equations shall be solved in each time step. This strategy is made possible through a sequence of model reformulations by exploring the models' thermodynamic structures. The generalized Onsager principle directly guides these reformulation procedures. In the reformulated but equivalent models, the reversible and irreversible components can be identified, guiding the numerical platform to decouple the reversible and irreversible dynamics. This eventually leads to decoupled numerical algorithms, given that the coupling terms only involve irreversible dynamics. To further demonstrate the numerical platform's power, we apply it to several specific incompressible hydrodynamic models. The energy stability of the proposed numerical schemes is shown in detail. The second-order accuracy in time is verified numerically through time step refinement tests. Several benchmark numerical examples are presented to further illustrate the proposed numerical framework's accuracy, stability, and efficiency.

A computational examination of large-scale pool fires: variations in crosswind velocity and pool shape

A high-fidelity large-eddy simulation and unsteady flamelet combustion model construct is deployed to numerically investigate the effects of crosswind magnitude and pool fire shape on large-scale pool fire attributes. These include general flame dynamics, flame shape and radiative flux magnitude in and around the fire. Three pool fire shapes at a nominal length scale of 10 m are subjected to four crosswind magnitudes between 0 and 20 m s $^{-1}$ . The pool shapes studied are circular, square and rectangular. The study includes the sensitivity of parameters to mesh and time step refinement. Results demonstrate that the rectangular shape, under crosswind, has low-levels of vertical velocity induction, resulting in a plume that is closer to the ground. In the quiescent regime, under-resolved meshes provide a higher radiative heat flux prediction compared with the most refined mesh. However, as crosswind increases, low mesh resolutions underpredicted radiative flux. This is due to the coarse mesh resolution not capturing small-scale vortical features that increased mixing and combustion efficiency. A transition of peak radiative flux with respect to crosswind occurs from the leeward- to windward-side of the pool, while sharp pool features result in larger radiative heat fluxes concentrated in regions of high scalar dissipation rate.

An arbitrary Lagrangian-Eulerian method for simulating interfacial dynamics between a hydrogel and a fluid

Hydrogels are crosslinked polymer networks swollen with an aqueous solvent, and play central roles in biomicrofluidic devices. In such applications, the gel is often in contact with a flowing fluid, thus setting up a fluid-hydrogel two-phase system. Using a recently proposed model (Young et al. [41] 2019), we treat the hydrogel as a poroelastic material consisting of a Saint Venant-Kirchhoff polymer network and a Newtonian viscous solvent, and develop a finite-element method for computing flows involving a fluid-hydrogel interface. The interface is tracked by using a fixed-mesh arbitrary Lagrangian-Eulerian method that maps the interface to a reference configuration. The interfacial deformation is coupled with the fluid and solid governing equations into a monolithic algorithm using the finite-element library deal.II. The code is validated against available analytical solutions in several non-trivial flow problems: one-dimensional compression of a gel layer by a uniform flow, two-layer shear flow, and the deformation of a Darcy gel particle in a planar extensional flow. In all cases, the numerical solutions are in excellent agreement with the analytical solutions. Numerical tests show second-order convergence with respect to mesh refinement, and first-order convergence with respect to time-step refinement.

Predicting large-scale pool fire dynamics using an unsteady flamelet- and large-eddy simulation-based model suite

A low-Mach, unstructured, large-eddy-simulation-based, unsteady flamelet approach with a generalized heat loss combustion methodology (including soot generation and consumption mechanisms) is deployed to support a large-scale, quiescent, 5-m JP-8 pool fire validation study. The quiescent pool fire validation study deploys solution sensitivity procedures, i.e., the effect of mesh and time step refinement on capturing key fire dynamics such as fingering and puffing, as mesh resolutions approach O(1) cm. A novel design-order, discrete-ordinate-method discretization methodology is established by use of an analytical thermal/participating media radiation solution on both low-order hexahedral and tetrahedral mesh topologies in addition to quadratic hexahedral elements. The coupling between heat losses and the flamelet thermochemical state is achieved by augmenting the unsteady flamelet equation set with a heat loss source term. Soot and radiation source terms are determined using flamelet approaches for the full range of heat losses experienced in fire applications including radiative extinction. The proposed modeling and simulation paradigm are validated using pool surface radiative heat flux, maximum centerline temperature location, and puffing frequency data, all of which are predicted within 10% accuracy. Simulations demonstrate that under-resolved meshes predict an overly conservative radiative heat flux magnitude with improved comparisons as compared to a previously deployed hybrid Reynolds-averaged Navier–Stokes/eddy dissipation concept-based methodology.

A refined two-scale model for Newtonian and non-Newtonian fluids in fractured poroelastic media

A refined model is presented which numerically resolves the fluid flow within a fracture inside a saturated poroelastic material. At the discontinuity, the mass balance of the fluid is solved using the velocity profile inside the fracture, with the velocity profile being determined numerically from the momentum balance in the integration points at the discontinuity. The resolution of the mass balance and the velocity profile in the integration points at the discontinuity is coupled to surrounding poroelastic material model in a two-scale approach. The resulting monolithic scheme allows for complex fluid behaviour to be included, which is demonstrated via the inclusion of the fluid inertia and the use of a Carreau fluid. The governing equations are discretised using T-splines, while the fracture is represented by spline-based interface elements, and an implicit time discretisation scheme is used. Mesh refinement studies are carried out for a typical case, which contains a pressurised and propagating fracture. A coarse macro-scale mesh is sufficient to obtain the correct propagation velocities, but a finer mesh is needed to prevent pressure oscillations. Time step refinement studies show the capabilities of the model to capture pressure waves within the fracture. Finally, it is shown that if inertial effects are limited and a Newtonian fluid is used, a fracture scale discretisation using a single element is sufficient. However, for a Carreau fluid or when the problem is inertia-dominated, smaller elements are needed to correctly represent the velocity profile.

Hydrodynamic Study and Performance Analysis of the OC4-DeepCWind Platform by CFD Method

(1) The RAOs of OC4-DeepCWind platform motions are more sensitive to the low-frequency wave than the high-frequency wave. The nonlinear motion responses for platform heave and pitch motions are comparatively remarkable. (2) The pitch motion of OC4-DeepCWind platform is much more apparently influenced by the height of center of gravity (COG) than surge and heave motions. The lower COG height within a suitable range leads to a smaller fluctuation amplitude of platform pitch motion in waves. (3) A large horizontal displacement abruptly occurs to the OC4-DeepCWind platform when one mooring line is failure. The risk of failure for the other mooring lines significantly increases. To better understand the hydrodynamic performance of a floating support platform in various wave environments, a two-phase CFD solver naoe-FOAM-SJTU based on the open source CFD toolbox OpenFOAM is applied to investigate the hydrodynamic characteristics and motion performance of the OC4-DeepCWind platform. Moreover, the restoring force and moment of mooring lines are simulated using the solver in time domain. The studies of grid sensitivity and time step refinement are first conducted to determine an appropriate time step and mesh size. Then hydrodynamic responses of the floater in free-decay tests are analyzed and compared with experimental data, and the motion performance of the platform in regular waves with different parameters is also investigated. In addition, the platform motion responses with one mooring line broken and different heights of center of gravity are explored. It is shown that simulation results have good agreement with published data, and several conclusions can be drawn through the study. The RAOs of platform motions are found to be more sensitive to the low-frequency wave than the high-frequency wave. Nonlinear motion responses are comparatively remarkable in platform heave and pitch motions. Besides, the lower height of center of gravity within a suitable range is benefit to the stability of floating platform. Survival condition with broken mooring line should be paid enough attention to avoid the failure of other mooring lines.

Comparison of fully-coupled and sequential solution methodologies for enhanced geothermal systems

The simulation of enhanced geothermal systems (EGS) requires finding the solution to a highly coupled nonlinear set of partial differential equations. Verification of EGS modelling has been difficult due to a lack of analytical or semi-analytical solutions and is typically limited to a subset of processes. A comparison of alternative solution schemes is one method to numerically verify the coupling between all processes simultaneously and provide confidence in the solutions produced by the different solution schemes. This article presents the first such comparison of monolithic and sequential solution schemes for the modelling of EGS wells. A custom thermo-hydro-mechanical finite element model for an EGS well connected by planar fractures is developed along with two monolithic schemes; one using the Newton–Raphson iterative method and another using the Newton–Raphson iterative method modified by Aitken’s Δ2 relaxation method. The two monolithic schemes are compared against two sequential schemes and one loosely-coupled scheme in terms of accuracy and computational efficiency. The monolithic schemes are shown to have an optimal rate of convergence with respect to mesh and timestep refinement. The monolithic Newton–Raphson scheme is shown to require small timesteps that violate the minimum timestep size of the β-method of time integration, introducing spurious oscillations into the solution. Aitken’s Δ2 relaxation method is shown to improve the ability of the Newton–Raphson scheme to find a converged solution at larger timesteps, avoiding spurious oscillations and reducing overall computation time. The sequential solution schemes are shown to provide optimal rates of convergence with respect to mesh and timestep refinement and converge to the same solution as the monolithic solution schemes, verifying the results of both methodologies and previous modelling efforts. It is demonstrated that the sequential solution schemes are computationally more expensive than the monolithic schemes when seeking well-converged solutions. The loosely-coupled scheme demonstrates an optimal rate of convergence with mesh and timestep refinement, but is also shown to be much less accurate than the other schemes. The recommended solution scheme for the efficient simulation of EGSs is thus the Newton–Raphson scheme with Aitken’s Δ2 relaxation.

Conservative multirate multiscale simulation of multiphase flow in heterogeneous porous media

Accurate and efficient simulation of multiphase flow in heterogeneous porous media motivates the development of space-time multiscale strategies for the coupled nonlinear flow (pressure) and saturation transport equations. The flow equation entails heterogeneous high-resolution (fine-scale) coefficients and is global (elliptic or parabolic). The time-dependent saturation profile, on the other hand, may exhibit sharp local gradients or discontinuities (fronts) where the solution accuracy is highly sensitive to the time-step size. Therefore, accurate flow solvers need to address the multiscale spatial scales, while advanced transport solvers need to also tackle multiple time scales. This paper presents the first multirate multiscale method for space-time conservative multiscale simulation of sequentially coupled flow and transport equations. The method computes the pressure equation at the coarse spatial scale with a multiscale finite volume technique, while the transport equation is solved by taking variable time-step sizes at different locations of the domain. At each coarse time step, the developed local time-stepping technique employs an adaptive recursive time step refinement to capture the fronts accurately. The applicability (accuracy and efficiency) of the method is investigated for a wide range of two-phase flow simulations in heterogeneous porous media. For the studied cases, the proposed method is found to provide 3 to 4 times faster simulations. Therefore, it provides a promising strategy to minimize the tradeoff between accuracy and efficiency for field-scale applications.

Code verification in computational geomechanics: Method of manufactured solutions (MMS)

Code verification is the process of ensuring that there are no coding mistakes or algorithm deficiencies in a scientific computing code. In this paper, code verification in computational geomechanics based on the method of manufactured solutions (MMS) is discussed. The application of this technique, which is the most rigorous method among all code verification techniques, is presented through verification of an in-house geomechanics code based on the Edge-based Smoothed Point Interpolation Method (ESPIM) for spatial, and three-point scheme for temporal discretisation of the governing equations. The MMS, as a generic approach to construct all-inclusive and complex test solutions for code verification, is explained. Two manufactured solutions along with an order of accuracy study are employed in several different setups to obtain the numerical and formal order of accuracies of the code. Results from the mesh and time step refinement studies confirm the verification of our code, and demonstrate the application of the MMS and order of convergence study in coupled problems in geomechanics.

Transient 3D finite element method for predicting extrudate swell of domains containing sharp edges

A new transient 3D finite element method for predicting extrudate swell of domains containing sharp edges is proposed. Here, the sharp edge is maintained over a large distance in the extrudate by describing the corner lines as material lines. The positions of these lines can be used to describe the transverse swelling of the 2D free surfaces and expand the domain over which a 2D height function on the free surfaces is applied. Solving the 2D height functions gives the positions of the free surfaces. First a 2D axisymmetric case was tested for comparison, using three different constitutive models. The Giesekus, linear Phan-Thien Tanner (PTT) and exponential PTT constitutive models all showed convergence upon mesh- and time-step refinement. It was found that convergence remains challenging due to the singularity at the die exit. The new method is validated by comparing the final volume change of the extrudate of a 3D cylinder to the final volume change of a reference mesh of the 2D axisymmetric case. Finally, simulations were performed for different, complex, die shapes for a viscous fluid and for viscoelastic fluids. The results compared favorably with literature. Viscoelastic results, using the Giesekus model and the exponential PTT model, were compared for different Weissenberg numbers and different values for the non-linear parameters of the constitutive models. It was found that the swell is highly dependent on the rheological parameters and the constitutive model used.

Numerical Simulations of Breaking Waves and Steep Waves Past a Vertical Cylinder at Different Keulegan–Carpenter Numbers

A three-dimensional (3D) numerical two-phase flow model based on solving unsteady Reynolds-averaged Navier–Stokes (URANS) equations has been used to simulate breaking waves and steep waves past a vertical cylinder on a 1:10 slope. The volume of fluid (VOF) method is employed to capture the free surface and the k–ω shear–stress transport (k–ω SST) turbulence model is used to simulate the turbulence effects. Mesh and time-step refinement studies have been conducted. The numerical results of wave forces on the structure are compared with the experimental data (Irschik et al., 2004, “Breaking Wave Loads on a Slender Pile in Shallow Water,” Coastal Engineering, Vol. 4, World Scientific, Singapore, pp. 568–581) to validate the numerical model, and the numerical results are in good agreement with the measured data. The wave forces on the structure at different Keulegan–Carpenter (KC) numbers are discussed in terms of the slamming force. The secondary load cycles are observed after the wave front past the structure. The dynamic pressure and velocity distribution, as well as the characteristics of the vortices around the structure at four important time instants, are studied.

Uncertainty estimation of a CFD-methodology for the performance analysis of a collective and cyclic pitch propeller

Estimation and analysis of the uncertainty introduced by using a numerical model for the investigation and study of any type of flow problem have become common industry practice. Through understanding and evaluation of the uncertainty introduced by a numerical model, the accuracy and applicability of the model itself are evaluated. In this paper, the numerical uncertainty of a CFD-methodology developed to analyse the hydrodynamic performance of a collective and cyclic pitch propeller (CCPP) is estimated and analysed. The CCPP is a novel propulsion and manoeuvring concept for autonomous underwater vehicles, aimed to generate both propulsion and manoeuvring forces through advanced control of the propeller's blade pitch. The numerical uncertainty is established for three performance parameters, the generated propulsive force, the side-force magnitude, and the side-force orientation, by conducting a grid and time-step refinement study over three operational conditions. Additionally, the influence of the oscillatory uncertainty, introduced by the periodic nature of the problem, is investigated although shown to have a minimal effect when properly monitored. Based on a least-squares regression analysis of the refined simulation results, the numerical uncertainty is proven to be dominated by the introduced discretisation errors. In the case of the propulsive and side-force magnitude, the total uncertainty is dictated by the time discretisation uncertainty under bollard pull conditions, while the total uncertainty of the captive cases is mainly a result of the spatial discretisation uncertainty. The total uncertainty in the side-force orientation is observed to be primarily a consequence of the time discretisation uncertainty for all simulated cases. Overall, the total uncertainty for captive cases can be considered satisfactory for all three performance parameters, while further work is needed to reduce the observed uncertainty of the simulations under bollard pull conditions.

Characteristics of higher-harmonic breaking wave forces and secondary load cycles on a single vertical circular cylinder at different Froude numbers

A 3D numerical two-phase flow model based on solving Unsteady Reynolds-averaged Navier-Stokes (URANS) equations has been used to simulate spilling breaking waves past a single vertical cylinder installed at the edge of a 1:10 slope. The volume of fluid (VOF) method is employed to capture the free surface, and the k−ω Shear-Stress Transport (k−ω SST) turbulence model is used to simulate turbulence effects. Mesh and time-step refinement studies have been conducted by examining the free surface elevations in front of the cylinder and the total horizontal wave forces on the cylinder. These two quantities have also been compared with experimental data from Irschik et al. [21] to validate the present numerical model. The present numerical results are in good agreement with the measured data. The process of breaking waves past the cylinder is described and explained. Moreover, the ‘secondary load cycle’ phenomenon which may lead to higher-harmonic structural response is observed. The characteristics of higher-harmonic wave forces are discussed further at different Froude numbers by changing cylindrical diameters and incident wave properties for both breaking and non-breaking wave cases. For non-breaking wave cases, pronounced secondary load cycles are observed for the Froude number exceeding about 0.4, and the magnitudes of the secondary loads are about 8.1%–12.6% of the peak-to-peak wave force (crest to trough value). However, for breaking wave cases, the critical value of the Froude number to observe the pronounced secondary load cycles reduces to 0.35. The relative magnitude of the secondary load to the peak-to-peak wave force reduces to a range of 2.4%–7.6% due to the high slamming force. The duration of the secondary load cycle is not greater than 0.242T for all the cases, where T is incident wave period. For a vertical cylinder with a larger diameter, the pronounced secondary load cycle occurs with a larger wave steepness.

Numerical study of wave run-up on a fixed and vertical surface-piercing cylinder subjected to regular, non-breaking waves using OpenFOAM

In wave-structure interaction, wave run-up is an important phenomenon that needs to be considered in the design of offshore structures. A thorough understanding of the physics of the nonlinear flow phenomena is necessary for the better insight into the run-up phenomenon. The present work, primarily, is focused on the hydrodynamic simulation of wave run-up and mainly seeks to evaluate the importance of wave scattering types identified by Swanet al. (2005) and lateral progressive edge waves on nonlinear wave amplification around a single fixed cylinder. The analysis is performed numerically using Navier-Stokes solver based on OpenFOAM over a range of wave steepnesses and wavelengths. In order to study wave-cylinder interaction, a numerical wave tank is utilized which includes IHFOAM toolbox for wave generation. For the simulations, the cylinder is assumed vertical and surface piercing with a circular cross-section and the incident wave is considered fifth-order Stokes, regular and non-breaking in deep water. The obtained numerical results that are mainly analyzed in this study include local surface elevation, the scattered wave field around the cylinder and the wave forces. A grid-space and time step refinement study are carried out to evaluate the performance and accuracy of numerical wave tank. The numerical results are compared with experimental data provided by ITTC (OEC), (2013) and a good agreement was obtained which indicates that OpenFOAM is very capable of accurate modeling of nonlinear wave interaction with offshore structures. The harmonic analysis for both short and long wave cases indicates that the scattered wave field around the cylinder involves high harmonics wave run-up. It is confirmed that even the wave-structure interaction caused by the linear incident wave, can lead to weakly nonlinear wave amplification around the cylinder. In the case of long waves cases, the harmonic loading involves at-least third harmonics due to the occurrence of the second wave run-up at the front part of the cylinder.

Linear, second order and unconditionally energy stable schemes for the viscous Cahn–Hilliard equation with hyperbolic relaxation using the invariant energy quadratization method

In this paper, we consider numerical approximations for the viscous Cahn–Hilliard equation with hyperbolic relaxation. This type of equations processes energy-dissipative structure. The main challenge in solving such a diffusive system numerically is how to develop high order temporal discretization for the hyperbolic and nonlinear terms, allowing large time-marching step, while preserving the energy stability, i.e. the energy dissipative structure at the time-discrete level. We resolve this issue by developing two second-order time-marching schemes using the recently developed “Invariant Energy Quadratization” approach where all nonlinear terms are discretized semi-explicitly. In each time step, one only needs to solve a symmetric positive definite (SPD) linear system. All the proposed schemes are rigorously proven to be unconditionally energy stable, and the second-order convergence in time has been verified by time step refinement tests numerically. Various 2D and 3D numerical simulations are presented to demonstrate the stability, accuracy, and efficiency of the proposed schemes.

A new numerical methodology for simulation of unstable crack growth in time independent brittle materials

This paper focuses on a new algorithm for the quasi-static simulation of unstable crack propagation in time-independent brittle material. The proposed approach is based on a combination of the standard Newton’s algorithm with a fictive path loading algorithm. First, a time step refinement process, associated to the Newton’s algorithm, enables us to efficiently detect the beginning of unstable crack extension. Then, a fictive path loading algorithm is performed to assess an estimation of the post instability solution. Finally, this estimated solution is used as initial guess in the standard Newton’s algorithm for the same time increment. The final converged solution satisfies both the mechanical equilibrium and the mechanical behaviour integration. This paper highlights several advantages of the proposed algorithm such as robustness, easy implementation on an existing Newton’s type solver, genericity for different unstable nonlinear physical problems. An application of this new method is presented for a smeared crack model, based on a local continuous damage formulation, devoted to the simulation of brittle nuclear fuel rupture. The efficiency of the proposed algorithm is demonstrated through a set of 2d and 3d simulation results.