Robust Numerical Implementation Research Articles

Mastering the art: Proficient finite element implementation and robust evaluation of a strain-hardening porous ductile material crack growth prediction model at finite strain

Accurate prediction of crack extension is crucial for ensuring the structural integrity of metal components exposed to diverse loads. While the Gurson model and its extensions are widely accepted for delineating ductile fracture stages, especially in porous materials with a rigid-perfectly plastic matrix, their efficacy diminishes in the presence of strain-hardening matrices. To address this limitation, our study introduces a robust numerical implementation and assessment of Perrin's model. In contrast to Gurson's approach, Perrin's model incorporates two distinct strain-hardening parameters derived from an intricate analysis of a strain-hardening hollow sphere subjected to axisymmetric loading.Moving beyond theoretical discussions, this paper actively engages in the rigorous evaluation of Perrin's model under various scenarios, encompassing isotropic, kinematic, and mixed isotropic-kinematic hardening conditions. The model's effectiveness is convincingly demonstrated through meticulous comparisons between numerical simulations and real-world experimental observations conducted on pre-cracked specimens.Facing challenges related to achieving global elastic-plastic convergence in large-scale simulations, our research introduces a sophisticated approach. Leveraging stiffness tangent moduli ensure the maintenance of quadratic convergence in global Newton method iterations. This not only enhances the reliability of our simulations but also underscores the practical applicability of our findings.The study also focuses on the necking and crack propagation of a cylindrical specimen through meticulous computational modeling. Our mesh captures intricate details and addresses stress gradients, revealing coplanar crack propagation contrary to prior beliefs. The confirmation of the cup-cone fracture effect challenges previous interpretations, emphasizing the critical role of void nucleation. This reinforces the credibility and applicability of our approach in advancing the understanding of material fracture behavior.In summary, this study significantly contributes to existing knowledge in fracture prediction models, offering a robust framework for predicting ductile fracture under substantial deformations. The presented findings pave the way for advancements in structural engineering, providing invaluable insights for optimizing the performance and resilience of metal structures in real-world applications.

Mapping of multiphase field model parameters to physical factors in order to simulate desired phase transformations

The multiphase field modelling is a powerful tool to understand sufficiently complicated material transformation phenomena. The numerical implementation of this modelling technique faced many challenges and finally Ohno and Matsuura (2009) provided a consistent model which overcame some of the numerical artifacts of its predecessors. However, the Ohno and Matsuura model has not had a robust numerical implementation. This difficulty in implementation is due to the complicated nature of equations that arise from the model. Furthermore, fitting the model for the desired problem needs a thorough understanding of the behaviour of the model parameters. Hence, this work uses an FEM based solver to understand the use of parameters in a multiphase field model and their effect on interface characteristics. Finally, it is shown how the parameters of the model can be tuned to simulate common types of phase transformation.

Non-local large-strain FFT-based formulation and its application to interface-dominated plasticity of nano-metallic laminates

This paper presents a novel formulation and its robust numerical implementation of strain-gradient (SG) crystal plasticity within a large-strain (LS) elasto-viscoplastic (EVP) fast Fourier transform (FFT)-based micromechanical model. The resulting non-local SG-LS-EVPFFT formulation is used to model and understand the process of kink band formation during layer-parallel compression of nano-metallic laminates (NMLs). NMLs are layered composites with nanoscale thicknesses, thus requiring consideration of the interaction between dislocations and interfaces within the micromechanical model. The length-scale parameter of the SG model is calibrated by simulating a double pile-up and comparing predictions to analytical solution. This required new expressions for the defect energy, resulting in more accurate double pile-up predictions. The calibrated SG-LS-EVPFFT model is then used to simulate layer-parallel compression of copper–niobium NML. Formation of kink bands is predicted, and the model is used to rationalize the microscopic mechanisms enabling the formation process. It is found that accumulation of dislocations at interfaces leads to activation of layer-parallel slip, which in turn leads to kink band formation.

A robust algorithm to test the observability of large linear systems with unknown parameters

This work proposes a robust algorithm to examine the observability of linear systems whose dynamic states and parameters are to be identified. The observability of a dynamical system plays a fundamental role in predicting whether it would be successful to use system identification methods to estimate the dynamic states and parameters of the system from a given set of input–output measurements. The motivation of the development of the suggested algorithm arises from the need to address the significant physical memory requirements of the standard implementation of the Observability Rank Condition (ORC). The high computational cost of the ORC substantially limits its applicability to real-world engineering systems, even when the underlying dynamics can be reasonably approximated as linear. The framework of the algorithm is obtained through the derivation of a recursive formula for the computation of the observability matrix of linear systems with unknown parameters. To further improve the efficiency, robust numerical implementations of the algorithm are achieved through random realizations of the dynamic states and parameters and the use of modular operations. The superior performance of the algorithm is demonstrated using several examples of large linear dynamical models of engineering systems containing up to thousands of dynamic states and parameters.

Numerical Simulation of Wave Runup and Overtopping for Short and Long Waves Using Staggered Grid Variational Boussinesq

In designing a numerical tool for simulating a wide variety of water waves, i.e. short to long waves, an accurate and robust wave model and numerical implementation are needed. Dispersion and nonlinearity are the two most important physical aspects that should be modeled accurately. To be applicable to simulate many coastal engineering applications, the numerical scheme should be capable of simulating wave runup and overtopping. In this paper, we extend the capability of a Boussinesq-type model called Variational Boussinesq (VB) model for simulating the runup and overtopping of water waves. To that end, the vertical layer of the fluid is modeled continuously by a linear combination of three functions. If two of these three functions have been incorporated in the previous numerical approximation called the SVB model, this paper discusses the improvement of SVB model by incorporating all the three functions. This approach improve the dispersive property of the SVB model due to its ability to simulate short waves up to kd = 20, compared to the previous model which was only up to kd = 7, where k denotes wave number and d water depth. Furthermore, the model is implemented numerically by using the staggered conservative scheme. In the new implementation, the model is switched to the non-dispersive Shallow Water Equations (SWE) when dealing with a dry area for runup and overtopping phenomena. The new implementation is tested against analytical solutions of soliton propagation and standing wave phenomenon; moreover, it is also tested against experimental data from hydrodynamic laboratories for simulating solitary wave breaking above a sloping bottom, composite beach, and in a structure for simulating overtopping phenomenon. The implementation is also tested against experimental data for simulating irregular wave propagation and runup above a fringing reef. The results of numerical simulation agree quite well with experimental data.

Staggered Conservative Scheme for 2-Dimensional Shallow Water Flows

Simulating discontinuous phenomena such as shock waves and wave breaking during wave propagation and run-up has been a challenging task for wave modeller. This requires a robust, accurate, and efficient numerical implementation. In this paper, we propose a two-dimensional numerical model for simulating wave propagation and run-up in shallow areas. We implemented numerically the 2-dimensional Shallow Water Equations (SWE) on a staggered grid by applying the momentum conserving approximation in the advection terms. The numerical model is named MCS-2d. For simulations of wet–dry phenomena and wave run-up, a method called thin layer is used, which is essentially a calculation of the momentum deactivated in dry areas, i.e., locations where the water thickness is less than the specified threshold value. Efficiency and robustness of the scheme are demonstrated by simulations of various benchmark shallow flow tests, including those with complex bathymetry and wave run-up. The accuracy of the scheme in the calculation of the moving shoreline was validated using the analytical solutions of Thacker 1981, N-wave by Carrier et al., 2003, and solitary wave in a sloping bay by Zelt 1986. Laboratory benchmarking was performed by simulation of a solitary wave run-up on a conical island, as well as a simulation of the Monai Valley case. Here, the embedded-influxing method is used to generate an appropriate wave influx for these simulations. Simulation results were compared favorably to the analytical and experimental data. Good agreement was reached with regard to wave signals and the calculation of moving shoreline. These observations suggest that the MCS method is appropriate for simulations of varying shallow water flow.

Fast Least-Squares Padé approximation of problems with normal operators and meromorphic structure

In this work, we consider the approximation of Hilbert space-valued meromorphic functions that arise as solution maps of parametric PDEs whose operator is the shift of an operator with normal and compact resolvent, e.g. the Helmholtz equation. In this restrictive setting, we propose a simplified version of the Least-Squares Pad\'e approximation technique introduced in [6] following [11]. In particular, the estimation of the poles of the target function reduces to a low-dimensional eigenproblem for a Gramian matrix, allowing for a robust and efficient numerical implementation (hence the "fast" in the name). Moreover, we prove several theoretical results that improve and extend those in [6], including the exponential decay of the error in the approximation of the poles, and the convergence in measure of the approximant to the target function. The latter result extends the classical one for scalar Pad\'e approximation to our functional framework. We provide numerical results that confirm the improved accuracy of the proposed method with respect to the one introduced in [6] for differential operators with normal and compact resolvent.

Robust numerical implementation of non-standard phase-field damage models for failure in solids

Recently several non-standard phase-field models different from the standard one for brittle fracture (Bourdin et al. 2000, 2008; Miehe et al. 2010a) have been proposed for the modeling of damage and fracture in solids. On the one hand, linear elastic behavior before the onset of damage and even cohesive cracking induced quasi-brittle failure can be considered by such non-standard phase-field models. On the other hand, they present great challenges to the numerical implementation since the damage boundedness is no longer automatically fulfilled. In this work, the numerical implementation of the unified phase-field damage model (Wu 2017, 2018) is addressed, though it also applies to the standard and other non-standard ones. In particular, an iterative alternate minimization (AM) algorithm enhanced with path-following strategies is presented. The phase-field or damage sub-problem is solved by the bound-constrained optimization solver in which the boundedness and irreversibility conditions of the crack phase-field are exactly enforced. Moreover, material softening induced snap-backs typically for localized failure in solids can also be effectively dealt with. The equilibrium paths in terms of the fracture surface and of the indirect displacement (e.g., the crack mouth opening or sliding displacement) are discussed in the context of the AM algorithm. The AM algorithm with the indirect displacement control is advocated, since both prescribed external forces and displacements can be naturally dealt with, which is more versatile than the one with the fracture surface control applicable only for prescribed displacements. Representative examples of benchmark tests show that the AM algorithm enhanced with the indirect displacement control is extremely robust even for large increment sizes. The insensitivity of the unified phase-field damage model to the incorporated length scale and mesh size is also confirmed.

Predicted Land Carbon Dynamics Are Strongly Dependent on the Numerical Coupling of Nitrogen Mobilizing and Immobilizing Processes: A Demonstration with the E3SM Land Model

AbstractWhile coupling carbon and nitrogen processes is critical for Earth system models to accurately predict future climate and land biogeochemistry feedbacks, it has not yet been analyzed how numerical methods that land biogeochemical models apply to couple soil mineral nitrogen mobilizing and immobilizing processes affect predicted ecosystem carbon and nitrogen cycling. These effects were investigated here by using the E3SM land model and an evaluation of three plausible and widely used numerical couplings: 1) the mineral nitrogen–based limitation scheme, 2) the net uptake–based limitation scheme, and 3) the proportional nitrogen flux–based limitation scheme. It was found that these three schemes resulted in large differences (with a range of 316 PgC) in predicted cumulative land–atmosphere carbon exchanges under the RCP4.5 atmospheric CO2 simulations. This large uncertainty is without accounting for the different representations of the many land biogeochemical processes, but is about 73% of the range (434 PgC) reported for CMIP5 RCP4.5 simulations. These results help explain the large uncertainty found in various model intercomparison experiments and suggest that more robust numerical implementations are needed to improve carbon–nutrient cycle coupling in terrestrial ecosystem models.

Numerical implementation of a cold-ion, Boltzmann-electron model for nonplanar plasma-surface interactions

Plasma-surface interactions are ubiquitous in the field of plasma science and technology. Much of the physics of these interactions can be captured with a simple model comprising a cold ion fluid and electrons which satisfy the Boltzmann relation. However, this model permits analytical solutions in a very limited number of cases. This paper presents a versatile and robust numerical implementation of the model for arbitrary surface geometries in cartesian and axisymmetric cylindrical coordinates. Specific examples of surfaces with sinusoidal corrugations, trenches, and hemi-ellipsoidal protrusions verify this numerical implementation. The application of the code to problems involving plasma-liquid interactions, plasma etching, and electron emission from the surface is discussed.

Nonlinear mechanics of non-rigid origami: an efficient computational approach.

Origami-inspired designs possess attractive applications to science and engineering (e.g. deployable, self-assembling, adaptable systems). The special geometric arrangement of panels and creases gives rise to unique mechanical properties of origami, such as reconfigurability, making origami designs well suited for tunable structures. Although often being ignored, origami structures exhibit additional soft modes beyond rigid folding due to the flexibility of thin sheets that further influence their behaviour. Actual behaviour of origami structures usually involves significant geometric nonlinearity, which amplifies the influence of additional soft modes. To investigate the nonlinear mechanics of origami structures with deformable panels, we present a structural engineering approach for simulating the nonlinear response of non-rigid origami structures. In this paper, we propose a fully nonlinear, displacement-based implicit formulation for performing static/quasi-static analyses of non-rigid origami structures based on 'bar-and-hinge' models. The formulation itself leads to an efficient and robust numerical implementation. Agreement between real models and numerical simulations demonstrates the ability of the proposed approach to capture key features of origami behaviour.

A depth semi-averaged model for coastal dynamics

The present work extends the semi-integrated method proposed by Antuono and Brocchini [“Beyond Boussinesq-type equations: Semi-integrated models for coastal dynamics,” Phys. Fluids 25(1), 016603 (2013)], which comprises a subset of depth-averaged equations (similar to Boussinesq-like models) and a Poisson equation that accounts for vertical dynamics. Here, the subset of depth-averaged equations has been reshaped in a conservative-like form and both the Poisson equation formulations proposed by Antuono and Brocchini [“Beyond Boussinesq-type equations: Semi-integrated models for coastal dynamics,” Phys. Fluids 25(1), 016603 (2013)] are investigated: the former uses the vertical velocity component (formulation A) and the latter a specific depth semi-averaged variable, ϒ (formulation B). Our analyses reveal that formulation A is prone to instabilities as wave nonlinearity increases. On the contrary, formulation B allows an accurate, robust numerical implementation. Test cases derived from the scientific literature on Boussinesq-type models—i.e., solitary and Stokes wave analytical solutions for linear dispersion and nonlinear evolution and experimental data for shoaling properties—are used to assess the proposed solution strategy. It is found that the present method gives reliable predictions of wave propagation in shallow to intermediate waters, in terms of both semi-averaged variables and conservation properties.

Application of a multi-component mean field model to the coarsening behaviour of a nickel-based superalloy

A multi-component mean field model has been applied to predict the particle evolution of the γ′ particles in the nickel based superalloy IN738LC, capturing the transition from an initial multimodal particle distribution towards a unimodal distribution. Experiments have been performed to measure the coarsening behaviour during isothermal heat treatments using quantitative analysis of micrographs. The three dimensional size of the γ′ particles has been approximated for use in simulation. A coupled thermodynamic/mean field modelling framework is presented and applied to describe the particle size evolution. A robust numerical implementation of the model is detailed that makes use of surrogate models to capture the thermodynamics. Different descriptions of the particle growth rate of non-dilute particle systems have been explored. A numerical investigation of the influence of scatter in chemical composition upon the particle size distribution evolution has been carried out. It is shown how the tolerance in chemical composition of a given alloy can impact particle coarsening behaviour. Such predictive capability is of interest in understanding variation in component performance and the refinement of chemical composition tolerances. It has been found that the inclusion of misfit strain within the current model formulation does not have a significant affect upon predicted long term particle coarsening behaviour. Model predictions show good agreement with experimental data. In particular, the model predicts a reduced growth rate of the mean particle size during the transition from bimodal to unimodal distributions.

Radiative transfer equation modeling by streamline diffusion modified continuous Galerkin method.

Optical tomography has a wide range of biomedical applications. Accurate prediction of photon transport in media is critical, as it directly affects the accuracy of the reconstructions. The radiative transfer equation (RTE) is the most accurate deterministic forward model, yet it has not been widely employed in practice due to the challenges in robust and efficient numerical implementations in high dimensions. Herein, we propose a method that combines the discrete ordinate method (DOM) with a streamline diffusion modified continuous Galerkin method to numerically solve RTE. Additionally, a phase function normalization technique was employed to dramatically reduce the instability of the DOM with fewer discrete angular points. To illustrate the accuracy and robustness of our method, the computed solutions to RTE were compared with Monte Carlo (MC) simulations when two types of sources (ideal pencil beam and Gaussian beam) and multiple optical properties were tested. Results show that with standard optical properties of human tissue, photon densities obtained using RTE are, on average, around 5% of those predicted by MC simulations in the entire/deeper region. These results suggest that this implementation of the finite element method-RTE is an accurate forward model for optical tomography in human tissues.

Kinematic assumptions and their consequences on the structure of field equations in continuum dislocation theory

Continuum Dislocation Theory (CDT) relates gradients of plastic deformation in crystals with the presence of geometrically necessary dislocations. Therefore, the dislocation tensor is introduced as an additional thermodynamic state variable which reflects tensorial properties of dislocation ensembles. Moreover, the CDT captures both the strain energy from the macroscopic deformation of the crystal and the elastic energy of the dislocation network, as well as the dissipation of energy due to dislocation motion. The present contribution deals with the geometrically linear CDT. More precise, the focus is on the role of dislocation kinematics for single and multi-slip and its consequences on the field equations. Thereby, the number of active slip systems plays a crucial role since it restricts the degrees of freedom of plastic deformation. Special attention is put on the definition of proper, well-defined invariants of the dislocation tensor in order to avoid any spurious dependence of the resulting field equations on the coordinate system. It is shown how a slip system based approach can be in accordance with the tensor nature of the involved quantities. At first, only dislocation glide in one active slip system of the crystal is allowed. Then, the special case of two orthogonal (interacting) slip systems is considered and the governing field equations are presented. In addition, the structure and symmetry of the backstress tensor is investigated from the viewpoint of thermodynamical consistency. The results will again be used in order to facilitate the set of field equations and to prepare for a robust numerical implementation.

Incorporating root hydraulic redistribution in CLM4.5: Effects on predicted site and global evapotranspiration, soil moisture, and water storage

AbstractWe implemented the Amenu‐Kumar model in the Community Land Model (CLM4.5) to simulate plant Root Hydraulic Redistribution (RHR) and analyzed its influence on CLM hydrology from site to global scales. We evaluated two numerical implementations: the first solved the coupled equations of root and soil water transport concurrently, while the second solved the two equations sequentially. Through sensitivity analysis, we demonstrate that the sequentially coupled implementation (SCI) is numerically incorrect, whereas the tightly coupled implementation (TCI) is numerically robust with numerical time steps varying from 1 to 30 min. At the site‐level, we found the SCI approach resulted in better agreement with measured evapotranspiration (ET) at the AmeriFlux Blodgett Forest site, California, whereas the two approaches resulted in equally poor agreement between predicted and measured ET at the LBA Tapajos KM67 Mature Forest site in Amazon, Brazil. Globally, the SCI approach overestimated annual land ET by as much as 3.5 mm d−1 in some grid cells when compared to the TCI estimates. These comparisons demonstrate that TCI is a more robust numerical implementation of RHR. However, we found, even with TCI, that incorporating RHR resulted in worse agreement with measured soil moisture at both the Blodgett Forest and Tapajos sites and degraded the agreement between simulated terrestrial water storage anomaly and Gravity Recovery and Climate Experiment (GRACE) observations. We find including RHR in CLM4.5 improved ET predictions compared with the FLUXNET‐MTE estimates north of 20° N but led to poorer predictions in the tropics. The biases in ET were robust and significant regardless of the four different pedotransfer functions or of the two meteorological forcing data sets we applied. We also found that the simulated water table was unrealistically sensitive to RHR. Therefore, we contend that further structural and data improvements are warranted to improve the hydrological dynamics in CLM4.5.

A smoothed cap model for variably saturated soils and its robust numerical implementation

Summary This paper deals with the numerical implementation of a cap model for unsaturated soils. It provides a brief review of existing cap model approaches, based on which an improved model formulated in terms of generalised effective stress and matric suction is derived and described in detail. Although the proposed model is a multisurface plasticity model, it can efficiently be implemented using only single-surface projections because of the smoothness of the model, which is obtained by construction. Numerical algorithms are provided for these single-surface stress projections, using a single-equation approach whenever possible. The robustness of the utilised single-equation approaches is enhanced by proposing problem-fitted start-up procedures based on investigations of the nonlinear projection equations. A comparison of the model response with extensive material test data is used to validate the model and to demonstrate the robust application of the approach to silty sands and low to medium plasticity clays. Copyright © 2015 John Wiley & Sons, Ltd.

On the numerical implementation of the higher-order strain gradient-dependent plasticity theory and its non-classical boundary conditions

The higher-order gradient plasticity theory is successful in explaining the size effects encountered at the micron and submicron length scale. Due to the incorporation of spatial gradients of one or more internal variables in these theories and the associated non-classical boundary conditions, special types of elements in the finite element method maybe necessary. This makes the numerical implementation of this higher-order theory not straightforward. In this paper, a robust but straightforward numerical implementation of higher-order gradient-dependent plasticity theories is presented. The novelty of this paper is in (1) the application of the meshless methods, particularly the moving weighted least square method, combined with the finite element method for the numerical computation of plastic strain gradients, and (2) the numerical implementation of different types of higher-order microscopic boundary conditions at internal/external surfaces, interfaces, and moving elastic–plastic boundaries. The proposed numerical implementation algorithms can be easily adapted in the implementation of any form of higher-order gradient-dependent constitutive models. Examples of applying the current numerical approach is demonstrated for capturing mesh-objective shear band formation and size effect and boundary layer formation in thin films on elastic substrates and metal matrix composites with embedded elastic inclusions through the consideration of stiff, intermediate, and soft interfaces.

Phase Field Modeling of Brittle and Ductile Fracture

AbstractThe phase field modeling of brittle fracture was a topic of intense research in the last few years and is now well‐established. We refer to the work [1‐3], where a thermodynamically consistent framework was developed. The main advantage is that the phase‐field‐type diffusive crack approach is a smooth continuum formulation which avoids the modeling of discontinuities and can be implemented in a straightforward manner by multi‐field finite element methods. Therefore complex crack patterns including branching can be resolved easily. In this paper, we extend the recently outlined phase field model of brittle crack propagation [1‐3] towards the analysis of ductile fracture in elastic‐plastic solids. In particular, we propose a formulation that is able to predict the brittle‐to‐ductile failure mode transition under dynamic loading that was first observed in experiments by Kalthoff and Winkler [4]. To this end, we outline a new thermodynamically consistent framework for phase field models of crack propagation in ductile elastic‐plastic solids under dynamic loading, develop an incremental variational principle and consider its robust numerical implementation by a multi‐field finite element method. The performance of the proposed phase field formulation of fracture is demonstrated by means of the numerical simulation of the classical Kalthoff‐Winkler experiment that shows the dynamic failure mode transition. (© 2013 Wiley‐VCH Verlag GmbH & Co. KGaA, Weinheim)

Fluctuating-surface-current formulation of radiative heat transfer: Theory and applications

We describe a novel fluctuating-surface current formulation of radiative heat transfer between bodies of arbitrary shape that exploits efficient and sophisticated techniques from the surface-integral-equation formulation of classical electromagnetic scattering. Unlike previous approaches to non-equilibrium fluctuations that involve scattering matrices---relating "incoming" and "outgoing" waves from each body---our approach is formulated in terms of "unknown" surface currents, laying at the surfaces of the bodies, that need not satisfy any wave equation. We show that our formulation can be applied as a spectral method to obtain fast-converging semi-analytical formulas in high-symmetry geometries using specialized spectral bases that conform to the surfaces of the bodies (e.g. Fourier series for planar bodies or spherical harmonics for spherical bodies), and can also be employed as a numerical method by exploiting the generality of surface meshes/grids to obtain results in more complicated geometries (e.g. interleaved bodies as well as bodies with sharp corners). In particular, our formalism allows direct application of the boundary-element method, a robust and powerful numerical implementation of the surface-integral formulation of classical electromagnetism, which we use to obtain results in new geometries, including the heat transfer between finite slabs, cylinders, and cones.