Numerical Difficulties Research Articles

Conservative, bounded, and nonlinear discretization of the Cahn-Hilliard-Navier-Stokes equations

The Cahn-Hilliard equation describes phase separation in a binary mixture, typically modeled with a phase variable that represents the concentration of one phase or the concentration difference between the two phases. Though the system is energetically driven toward solutions within the physically meaningful range of the phase variable, numerical methods often struggle to maintain these bounds, leading to physically invalid quantities and numerical difficulties. In this work, we introduce a novel splitting and discretization for the Cahn-Hilliard equation, coupled with the Navier-Stokes equations, which inherently preserves the bounds of the phase variable. This approach transforms the fourth-order Cahn-Hilliard equation into a second-order Helmholtz equation and a second-order nonlinear equation with implicit energy barriers, which is reformulated and solved with a safeguarded optimization-based solution method. Our scheme ensures the phase variable remains in the valid range, robustly handles large density ratios, conserves mass and momentum, maintains consistency between these quantities, and achieves second-order accuracy. We demonstrate the method's effectiveness through a variety of studies of two-dimensional, two-phase fluid mixtures.

Gradient projection method for enforcing crack irreversibility as box constraints in a robust monolithic phase-field scheme

A phase-field monolithic scheme based on the gradient projection method is developed to model crack propagation in brittle materials under cyclic loading. As a type of active set method, the gradient projection method is particularly attractive to enforce the irreversibility condition imposed on the phase-field variables as bound constraints, or box constraints. This method has the advantages of allowing the rapid change of active constraints during iterations and computing the projected gradient with a negligible cost. The gradient projection method is further combined with the limited-memory BFGS (L-BFGS) method to overcome the convergence difficulties arising from the non-convex energy functional. A compact representation of the BFGS matrix is adopted as the limited-memory feature to avoid the storage of fully dense matrices, making this method practical for large-scale finite element simulations. By locating the generalized Cauchy point on the piecewise linear path formed by the projected gradient, the active set of box constraints can be determined. The variables in the active set, which are at the boundary of the box constraints, are kept fixed to form a subspace minimization problem. A primal approach and a dual approach are presented to solve this subspace minimization problem for the remaining free variables at the generalized Cauchy point. Several two-dimensional (2D) and three-dimensional (3D) examples are provided to demonstrate the capabilities of the proposed monolithic scheme, particularly in enforcing the phase-field irreversibility during crack propagation under cyclic loading. In these numerical examples, the proposed monolithic scheme is combined with an adaptive mesh refinement technique to alleviate the heavy computational cost incurred by the fine mesh resolution required around the crack region. The proposed method is further compared with two other phase-field solving techniques regarding the convergence behavior. To ensure a fair comparison, the same problem settings and implementation techniques are adopted. The proposed monolithic scheme provides a unified framework to overcome the numerical difficulties associated with the non-convex energy functional, effectively enforce the phase-field irreversibility to ensure the thermodynamic consistency, and alleviate the heavy computational cost through adaptive mesh refinement in 2D and 3D phase-field crack simulations.

Numerical Simulation of Free Surface Deformation and Melt Stirring in Induction Melting Using ALE and Level Set Methods.

This study investigates fixed and moving mesh methodologies for modeling liquid metal-free surface deformation during the induction melting process. The numerical method employs robust coupling of magnetic fields with the hydrodynamics of the turbulent stirring of liquid metal. Free surface tracking is implemented using the fixed mesh level set (LS) and the moving mesh arbitrary Lagrangian-Eulerian (ALE) formulation. The model's geometry and operating parameters are designed to replicate a semi-industrial induction melting furnace. Six case studies are analyzed under varying melt masses and coil power levels, with validation performed by comparing experimentally measured free surface profiles and magnetic field distributions. The melt's stirring velocity and recirculation patterns are also examined. The comparative analysis determines an improved performance of the ALE method, convergence, and computational efficiency. Experimental validation confirms that the ALE method reproduces the free surface shape more precisely, avoiding unrealistic topological changes observed in LS simulations. The ALE method faces numerical convergence difficulties for high-power and low-mass filling cases due to mesh element distortion. The proposed ALE-based simulation procedure is a potential numerical optimization tool for enhancing induction melting processes, offering scalable and robust solutions for industrial applications.

A non-interior-point continuation method for the optimal control problem with equilibrium constraints

This study presents a numerical method for the optimal control problem with equilibrium constraints (OCPEC). It is extremely difficult to solve OCPEC owing to the absence of constraint regularity and strictly feasible interior points. To solve OCPEC efficiently, we first relax the discretized OCPEC to recover the constraint regularity and then map its Karush–Kuhn–Tucker (KKT) conditions into a parameterized system of equations. Subsequently, we solve the parameterized system using a novel two-stage solution method called the non-interior-point continuation method. In the first stage, a non-interior-point method is employed to find an initial solution, which solves the parameterized system using Newton’s method and globalizes convergence using a dedicated merit function. In the second stage, a predictor–corrector continuation method is utilized to track the solution trajectory as a function of the parameter, starting at the initial solution. The proposed method regularizes the KKT matrix and does not enforce iterates to remain in the feasible interior, which mitigates the numerical difficulties in solving OCPEC. Convergence properties are analyzed under certain assumptions. Numerical experiments demonstrate that the proposed method can solve OCPEC while demanding remarkably less computation time than the interior-point method.

A pseudo-transient pressure gradient method for solving the incompressible Navier–Stokes equations

The fundamental problem behind the numerical solutions of incompressible viscous flow equations lies in the lack of an update equation for the pressure field. In essence, there is no explicit relation that can be utilized to update the pressure field from a given velocity field. Following the principle of the artificial compressibility method and considering that it is the pressure gradient that drives incompressible flows, this study introduces a new primitive variables method to solve the incompressible Navier–Stokes equations. This method ensures that the continuity condition is rigorously satisfied at every grid point within the flow domain. Its principle lies in modifying the continuity equation to include a pseudo-transient pressure gradient term that vanishes upon convergence, which allows the incompressibility constraint to be satisfied to machine zero. This method eliminates the numerical difficulties associated with pressure-based solvers and improves the convergence of all dependent variables. Numerical solutions for three steady-state validation problems, namely fully developed channel flow, lid-driven cavity flow (2D and 3D), and flow in a constricted channel, are obtained to validate the proposed method. The results show good agreement with solutions available in the literature and demonstrate substantially improved convergence.

A New Nonlinear Time-Domain Flutter Analysis Approach for Distorted Flows

Abstract This paper aims to establish an efficient and accurate computational fluid dynamic (CFD) method for calculating the flutter stability of fan blades in the presence of inlet distortion due to crosswind. Due to the asymmetry of the flow and interaction of fan and distortion, this type of analysis will require a whole assembly computational model. Therefore, reducing the computational cost needed to obtain accurate aerodynamic damping values is paramount. Mode-tracking free-flutter analyses are computationally expensive and troublesome due to the presence of many frequencies in the response signals and, hence, require long physical times to converge the aerodynamic damping. Moreover, the unsteady flow field associated with the distortion can be more significant than the one due to flutter, giving rise to additional numerical difficulties. On the other hand, the energy method approach is commonly performed for a single nodal diameter (ND) at a time, requiring many computations to establish the least stable ND. This paper proposes an alternative approach based on the multi-ND energy method. The underlying hypothesis is that the unsteady pressure scales linearly with the blade displacements, and the scattering of circumferential modes due to the distortion is negligible. The general method presented here is general and can be used to compute the aerodynamic damping for other types of distortion (such as downstream distortions due to pylon and structure).

Fast time-stepping discontinuous Galerkin method for the subdiffusion equation

Abstract The nonlocality of the fractional operator causes numerical difficulties for long time computation for time-fractional evolution equations. This paper develops a high-order fast time-stepping discontinuous Galerkin (DG) finite element method for a time-fractional diffusion equation, which saves storage and computational time. An optimal error estimate of the form $O(N^{-p-1} + h^{m+1} + \varepsilon N^{r\alpha })$ for the proposed time-stepping DG method is rigorously proved, where $N$ denotes the number of time intervals, $p$ is the degree of polynomial approximation on each time subinterval, $h$ is the maximum spatial mesh size, $r\ge 1$, $m$ is the order of finite element space and $\varepsilon>0$ can be arbitrarily small. Numerical simulations verify the theoretical analysis.

Demon Registration for 2D Empirical Wavelet Transforms

The empirical wavelet transform is a fully adaptive time-scale representation that has been widely used in the last decade. Inspired by the empirical mode decomposition, it consists of filter banks based on harmonic mode supports. Recently, it has been generalized to build the filter banks from any generating function using mappings. In practice, the harmonic mode supports can have a low-constrained shape in 2D, leading to numerical difficulties to estimate mappings adapted to the construction of empirical wavelet filters. This work aims to propose an efficient numerical scheme to compute empirical wavelet coefficients using the demons registration algorithm. Results show that the proposed approach is robust, accurate, and continuous wavelet filters permitting reconstruction with a low signal-to-noise ratio. An application for texture segmentation of scanning tunneling microscope images is also presented.

A fully-hybrid finite element formulation for general compressible-incompressible elasticity problems using de Rham compatible spaces

This work proposes a novel fully-hybrid finite element formulation for general elasticity problems, combining $H(\text{div},\Omega)$ conforming vector functions for displacement and $L^2(\Omega)$ discontinuous scalar functions for pressure. This pair is De Rham compatible, which means that within the incompressible regime, the divergence-free constrain will hold strongly at element level. As the $H(\text{div})$ spaces only present continuity of the normal displacements across elements boundaries, the tangential displacement continuity can be weakly imposed performing a hybridization of the shear stresses using $H^1(\partial\Omega_e)$ functions. From past researches conducted at LabMeC, this approach has demonstrated to pose some numerical difficulties as it leads to a saddle point problem with two constraint variables - the pressure and the shear-stresses. In this work, a second hybridization is done by approximating the tangential component of the primal displacement variable using $H^1(\partial\Omega_e)$ space. Two benchmarks are used to verify the developed numerical scheme - the classical Cook's membrane and a tridimensional cantilever beam subjected to an end shear load. Optimum convergence rate is achieved under compressible, quasi-incompressible and full incompressible scenarios.

Critical Assessment of Quadratic Closures for Prediction of Corner Flow Separation

The work proposed in this paper aims at improving the physical understanding and modeling of junction flows by conducting an analysis of several quadratic constitutive relations (QCRs) applied to corner vortices. Reynolds-averaged Navier–Stokes computations are performed on two academic configurations for verification purposes and on a technical application of wing–body juncture featuring corner flows to achieve a better comprehension of the QCR variants’ ability to reproduce turbulent stress combinations relevant to secondary flows generation. This addresses the more theoretical interest of analyzing the onset of corner separation (at an angle of attack α=11 deg) and the role of corner flow vortices. Linear eddy viscosity model and the extended QCR failed to detect these vortical structures; the latter encountered numerical stability difficulties. On the other hand, QCR2000, QCR2020, and QCR(r) showed behavior consistent with a simulation using the Reynolds Stress Model, experiments, and literature results. This better prediction was found to be related to the ability of these models to accurately predict the stress-induced vortex that acts in delaying and reducing the size of the corner separation. Notably, the more recent QCR(r) and QCR2020 have enhanced the estimation of the normal stress combination v′2¯−w′2¯ compared to the more widespread QCR2000, demonstrating their validity as alternative models.

Stress shadow effects in multistage horizontal hydrofracturing of tight reservoirs: a numerical analysis considering perforation cluster spacings and fracturing sequences

AbstractMultistage horizontal fracturing technology of reservoirs has been widely used to enhance tight hydrocarbon resource recovery. Determining the proper perforation cluster spacing is crucial to developing a complex fracture network that connects natural rock fractures and facilitates gas flow in reservoirs. The stress shadow effect that occurs between multiple clusters has a significant effect on the development of the fracture network in reservoirs. The quantification of the stress shadow effect and its influences on the development of fracture networks has not been resolved satisfactorily due to experimental and numerical difficulties with detecting and identifying crack propagation and intersection in deep reservoirs. In this study, we used an adaptive finite element-discrete element method to analyze crack propagation and intersection in tight shale reservoirs by altering perforation spacings and fracture sequences. The effects of five perforation cluster spacings using sequential, simultaneous and parallel fracturing techniques were compared. The non-linear properties of shales, hydromechanical coupling and fluid leak-off effects, proppant transport, dual-rupture criteria of strength and energy, and intersection of hydraulic fractures were taken into account in the simulation. The areas affected by the stress shadow and the optimum perforation cluster spacing in terms of hydraulic fracture propagation lengths, volumes, and microseismic magnitudes were evaluated when adopting different fracturing methods. This study extends the understanding of the impact of the different independent factors that contribute to the stress shadow effect and, therefore, provides new information to help guide wellbore completion design in horizontal fracturing wells.

Preserving large-scale features in simulations of elastic turbulence

Simulations of elastic turbulence, the chaotic flow of highly elastic and inertialess polymer solutions, are plagued by numerical difficulties: the chaotically advected polymer conformation tensor develops extremely large gradients and can lose its positive-definiteness, which triggers numerical instabilities. While efforts to tackle these issues have produced a plethora of specialized techniques – tensor decompositions, artificial diffusion, and shock-capturing advection schemes – we still lack an unambiguous route to accurate and efficient simulations. In this work, we show that even when a simulation is numerically stable, maintaining positive-definiteness and displaying the expected chaotic fluctuations, it can still suffer from errors significant enough to distort the large-scale dynamics and flow structures. We focus on two-dimensional simulations of the Oldroyd-B and FENE-P equations, driven by a large-scale cellular body forcing. We first compare two positivity-preserving decompositions of the conformation tensor: symmetric square root (SSR) and Cholesky with a logarithmic transformation (Cholesky-log). While both simulations yield chaotic flows, only the latter preserves the pattern of the forcing, i.e. its fluctuating vortical cells remain ordered in a lattice. In contrast, the SSR simulation exhibits distorted vortical cells that shrink, expand and reorient constantly. To identify the accurate simulation, we appeal to a hitherto overlooked mathematical bound on the determinant of the conformation tensor, which unequivocally rejects the SSR simulation. Importantly, the accuracy of the Cholesky-log simulation is shown to arise from the logarithmic transformation. We also consider local artificial diffusion, a potential low-cost alternative to high-order advection schemes. Unfortunately, the artificially enhanced diffusive smearing of polymer stress in regions of intense stretching substantially modifies the global dynamics. We then show how the spurious large-scale motions, identified here, contaminate predictions of scalar mixing. Finally, we discuss the effects of spatial resolution, which controls the steepness of gradients in a non-diffusive simulation.

Technology Innovation for Discovering Renal Autoantibodies in Autoimmune Conditions.

Autoimmune glomerulonephritis is a homogeneous area of renal pathology with clinical relevance in terms of its numerical impact and difficulties in its treatment. Systemic lupus erythematosus/lupus nephritis and membranous nephropathy are the two most frequent autoimmune conditions with clinical relevance. They are characterized by glomerular deposition of circulating autoantibodies that recognize glomerular antigens. Technologies for studying renal tissue and circulating antibodies have evolved over the years and have culminated with the direct analysis of antigen-antibody complexes in renal bioptic fragments. Initial studies utilized renal microdissection to obtain glomerular tissue. Obtaining immunoprecipitates after partial proteolysis of renal tissue is a recent evolution that eliminates the need for tissue microdissection. New technologies based on 'super-resolution microscopy' have added the possibility of a direct analysis of the interaction between circulating autoantibodies and their target antigens in glomeruli. Peptide and protein arrays represent the new frontier for identifying new autoantibodies in circulation. Peptide arrays consist of 7.5 million aligned peptides with 16 amino acids each, which cover the whole human proteome; protein arrays utilize, instead, a chip containing structured proteins, with 26.000 overall. An example of the application of the peptide array is the discovery in membranous nephropathy of many new circulating autoantibodies including formin-like-1, a protein of podosomes that is implicated in macrophage movements. Studies that utilize protein arrays are now in progress and will soon be published. The contribution of new technologies is expected to be relevant for extending our knowledge of the mechanisms involved in the pathogenesis of several autoimmune conditions. They may also add significant tools in clinical settings and modify the therapeutic handling of conditions that are not considered to be autoimmune.

Computational toolbox for scattering of focused light from flattened or elongated particles using spheroidal wavefunctions

T-matrix methods, with incident and scattered fields described as sums of multipolar fields, are attractive computational methods for many scattering problems due to their versatility, accuracy, and computational efficiency, especially for repeated calculations. However, numerical difficulties often hamper their use for non-spherical particles with large aspect ratios. Further, even if far-field scattering can be accurately calculated, it can be impossible to accurately calculate near-fields. The use of spheroidal wavefunctions, instead of the commonly-used spherical wavefunctions, can be a useful solution for these problems. However, the mathematical complexity of spheroidal wavefunctions, and the challenging task of accurate numerical calculation of them, are significant barriers to their use. We have developed a computational package of MATLAB routines for performing electromagnetic scattering calculations using spheroidal wavefunctions. These allow the determination of light scattering by non-spherical particles with high aspect ratios that would be inaccessible for double precision calculation using spherical wavefunctions. We demonstrate that our codes can be successfully used for cylinders of aspect ratios from 1/10 and 10, and for metal nanoparticles. We include routines for interoperability with regular T-matrix codes and packages such as our optical tweezers toolbox.

Uncertainty in dynamic econometric input-output models: a Norwegian case study

Input-output models used for macroeconomic impact and policy analysis are often characterised by large data sets and resource-heavy computing. However, the types of results they provide are sensitive to uncertainty in their core assumptions. Various approaches have been proposed to account for the uncertainty in one or more parts of the analysis assumptions. Although it is standard practice to include varied cases in policy analysis, the notion of stochasticity of parameters across periods of time is less widespread. Costly numerical computing and difficulty in interpreting the results complicate this approach. In this work, we adapt an environmentally-extended dynamic econometric input-output model to account for uncertainty of various kinds. We demonstrate the methodology in the case of building renovation. The results provide insight into situations where explicitly considering uncertainty as part of the analysis is useful, as well as others in which no additional information was gained from such treatment.

Evaluating Fracture Energy Predictions Using Phase-Field and Gradient-Enhanced Damage Models for Elastomers

Abstract Recently, the phase-field method has been increasingly used for brittle fractures in soft materials like polymers, elastomers, and biological tissues. When considering finite deformations to account for the highly deformable nature of soft materials, the convergence of the phase-field method becomes challenging, especially in scenarios of unstable crack growth. To overcome these numerical difficulties, several approaches have been introduced, with artificial viscosity being the most widely utilized. This study investigates the energy release rate due to crack propagation in hyperelastic nearly-incompressible materials and compares the phase-field method and a novel gradient-enhanced damage (GED) approach. First, we simulate unstable loading scenarios using the phase-field method, which leads to convergence problems. To address these issues, we introduce artificial viscosity to stabilize the problem and analyze its impact on the energy release rate utilizing a domain J-integral approach giving quantitative measurements during crack propagation. It is observed that the measured energy released rate during crack propagation does not comply with the imposed critical energy release rate, and shows non-monotonic behavior. In the second part of the paper, we introduce a novel stretch-based GED model as an alternative to the phase-field method for modeling crack evolution in elastomers. It is demonstrated that in this method, the energy release rate can be obtained as an output of the simulation rather than as an input which could be useful in the exploration of rate-dependent responses, as one could directly impose chain-level criteria for damage initiation. We show that while this novel approach provides reasonable results for fracture simulations, it still suffers from some numerical issues that strain-based GED formulations are known to be susceptible to.

Reinforcement topology optimization considering the dynamic instability

AbstractThe present study develops a new topology optimization scheme considering the dynamic instability caused by the unsymmetrical properties of system. From a mathematical point of view, the left and right eigenvectors of asymmetric system are observed with the complex eigenvalues. With the dynamic instability, the magnitudes of structural responses are increasing with respect to time and this phenomenon causes many engineering issues. As the dynamic instability is one of the serious problems, the suppression is desired from an engineering point of view. To systematically reduce this dynamic instability, the present study develops a new topology optimization scheme for the reinforcement design. To overcome the numerical difficulties of the mode conversion and the highly nonlinear behavior, this research proposes the summation of the first several complex eigenvalues. To show the issues of the dynamic instability and the validity of the present approach, several topological reinforcement problems are solved.

A stochastic algorithm for the ParaTuck decomposition

This paper introduces a novel stochastic algorithm for the ParaTuck Decomposition (PTD), addressing the challenge of local minima encountered in the traditional alternating least squares (ALS) approach. The proposed method integrates stochastic steps into the ALS framework to avoid the common swamp problems, where numerical difficulties prevent accurate decompositions. Our simulations indicate good convergence properties for PTD, suggesting a potential increase in the efficiency and reliability of this tensor decomposition across various applications.

Effect of Cold Vapor Application on Postoperative Sore Throat, Dry Throat, and Difficulty Swallowing: A Randomized Controlled Trial.

Background: Postoperative sore throat, dry throat, and difficulty swallowing are frequently reported complications following general anesthesia, imposing discomfort and challenges to patients during their recovery. Objective: The aim of this study was to determine the effect of cold vapor application for 15 minutes at 0, 2, and 6 hours after extubation on sore throat, dry throat, and difficulty swallowing. Methods: This was a randomized controlled experimental study. The research was carried out with 64 patients who underwent cholecystectomy in the operating room and general surgery clinic of a state hospital between December 2016 and August 2017. Research data were collected using the American Society of Anesthesiologists classification, Mallampati scoring, patient identification form, intraoperative patient evaluation form, Ramsey Sedation Scale, Numerical Rating Scale, dry throat and difficulty swallowing evaluation form, and cuff pressure gauge. Results: It was found that cold vapor application at the 0th, 2nd, and 6th hours postoperatively had no effect on sore throat and swallowing difficulty. Mean dry throat scores of the patients in the intervention group at the second and sixth hours were significantly lower than those in the control group. At the 24th hour, there was no statistically significant difference between the mean scores of sore throat, swallowing difficulty, and dry throat in the intervention and control groups. Conclusion: According to the results of this study, it was found that cold vapor application in the early postoperative period had no effect on sore throat and swallowing difficulty but reduced throat dryness at the second and sixth hours. Implications for Nursing: In the early postoperative period, cold vapor can be applied to relieve hoarseness due to extubation.

Diffraction by gratings: from the C-method to thestochastic C-method.

The Chandezon, or C-method, is an efficient and versatile numerical method for modeling diffraction problems involving smooth surface relief gratings. For some grating profiles, the C-method is limited when the height-to-period ratio exceeds a factor of three. This is due to the formation of ill-conditioned matrices for inversion. Here, the stochastic C-method (SCM) is introduced as a solution that leverages stochastic differential equations to overcome these numerical difficulties. The SCM is developed by altering the physical model of the grating profile function to include an additive Brownian noise component. The inclusion of noise dramatically expands the applicability of the C-method and enriches the physical model. Numerical experiments show that the SCM achieves a precision on the order of 10-5 for diffracted/transmitted amplitudes on sinusoidal profiles with height-to-period ratios as high as 72. These results are in agreement with those obtained using multi-precision and the rigorous coupled wave analysis (RCWA).