Numerical Flux Research Articles

DOT-type schemes for hybrid hyperbolic problems arising from free-surface, mobile-bed, shallow-flow models

Free-surface, mobile-bed, shallow-flow models may present Hybrid hyperbolic systems of partial differential equations characterised by conservative and non-conservative fluxes that can only be expressed in primitive variables. This paper presents the effort we made to derive DOT-type schemes (Osher-type schemes derived by Dumbser and Toro, 2011 [1]) for these kinds of systems formulated for the one-dimensional case. Firstly, for a Hybrid system, we managed to write a quasi-linear form characterised by the presence of a matrix, expressed as a function of the primitive variables, that multiplies the spatial derivative of the conserved variables. Next, we derived the first numerical flux by adapting the approach of Leibinger et al., 2016 [2] to this quasi-linear form. We called this result DOTHCP flux. To achieve a faster algorithm, instead of using an integration path in the space of conserved variables, as in the previous case, we employed a path in the space of primitive variables. We called this second formulation DOTHPP flux. Subsequently, we managed to account for certain physical constraints arising from the generalised Rankine-Hugoniot relations in the expression of one term of the previous flux formulation, thus obtaining the DOTHZR flux. Finally, we showed that these methods can also be applied to Combined systems characterised by conservative and non-conservative fluxes expressed in conserved variables. Several tests show the characteristics and good performances of the proposed methods when applied to Riemann problems of Hybrid and Combined systems deriving from free-surface models. Finally, thanks to the general formulation of the proposed DOT-type fluxes, these can also be applied to Hybrid and Combined hyperbolic systems deriving from different physical problems.

A fully well-balanced hydrodynamic reconstruction

Abstract The present work focuses on the numerical approximation of the weak solutions of the shallow water model over a non-flat topography. In particular, we pay close attention to steady solutions with nonzero velocity. The goal of this work is to derive a scheme that exactly preserves these stationary solutions, as well as the commonly preserved lake at rest steady solution. These moving steady states are solution to a nonlinear equation. We emphasize that the method proposed here never requires solving this nonlinear equation; instead, a suitable linearization is derived. To address this issue, we propose an extension of the well-known hydrostatic reconstruction. By appropriately defining the reconstructed states at the interfaces, any numerical flux function, combined with a relevant source term discretization, produces a well-balanced scheme that preserves both moving and non-moving steady solutions. This eliminates the need to construct specific numerical fluxes. Additionally, we prove that the resulting scheme is consistent with the homogeneous system on flat topographies, and that it reduces to the hydrostatic reconstruction when the velocity vanishes. To increase the accuracy of the simulations, we propose a well-balanced high-order procedure, which still does not require solving any nonlinear equation. Several numerical experiments demonstrate the effectiveness of the numerical scheme.

A geometrically and thermodynamically compatible finite volume scheme for continuum mechanics on unstructured polygonal meshes

We present a novel Finite Volume (FV) scheme on unstructured polygonal meshes that is provably compliant with the Second Law of Thermodynamics and the Geometric Conservation Law (GCL) at the same time. The governing equations are provided by a subset of the class of symmetric and hyperbolic thermodynamically compatible (SHTC) models introduced by Godunov in 1961. Specifically, our numerical method discretizes the equations for the conservation of momentum, total energy, distortion tensor and thermal impulse vector, hence accounting in one single unified mathematical formalism for a wide range of physical phenomena in continuum mechanics, spanning from ideal and viscous fluids to hyperelastic solids. By means of two conservative corrections directly embedded in the definition of the numerical fluxes, the new schemes are proven to satisfy two extra conservation laws, namely an entropy balance law and a geometric equation that links the distortion tensor to the density evolution. As such, the classical mass conservation equation can be discarded. Firstly, the GCL is derived at the continuous level, and subsequently it is satisfied by introducing the new concepts of general potential and generalized Gibbs relation. The new potential is nothing but the determinant of the distortion tensor, and the associated Gibbs relation is derived by introducing a set of dual or thermodynamic variables such that the GCL is retrieved by dot multiplying the original system with the new dual variables. Once compatibility of the GCL is ensured, thermodynamic compatibility is tackled in the same manner, thus achieving the satisfaction of a local cell entropy inequality. The two corrections are orthogonal, meaning that they can coexist simultaneously without interfering with each other. The compatibility of the new FV schemes holds true at the semi-discrete level, and time integration of the governing PDE is carried out relying on Runge-Kutta schemes. A large suite of test cases demonstrates the structure preserving properties of the schemes at the discrete level as well.

A finite volume scheme for solving an alternative morphodynamic model with porosity

AbstractThis work proposes an extended numerical scheme based on the finite volume method that is intended to solve a new morphodynamic model with porosity. The St-Venant equations and the Exner equation are coupled in this model, accounting for the variation of the porosity function . In the morphodynamic system, the Jacobian matrix's sign is taken into account when reconstructing the numerical fluxes using an improved Roe's approach to assure correct numerical computations.

WENO scheme on characteristics for the equilibrium dispersive model of chromatography with generalized Langmuir isotherms

Column chromatography is a laboratory and industrial technique used to separate different substances mixed in a solution. Mathematically, it can be modeled using non-linear partial differential equations whose main ingredients are the adsorption isotherms, which are non-linear functions modeling the affinity between the different substances in the solution and the solid stationary phase filling the column. The goal of this work is twofold. Firstly, we aim to extend the techniques of Donat, Guerrero and Mulet (2018) [3] to other adsorption isotherms. In particular, we propose a family of generalized Langmuir-type isotherms and prove that the correspondence between the concentrations of solutes in the liquid phase (the primitive variables) and the conserved variables is well defined and admits a global smooth inverse that can be computed numerically. Secondly, to establish the well-posedness of the mathematical model, we study the eigenstructure of the Jacobian of the mentioned correspondence and use this characteristic information to get oscillation-free sharp interfaces on the numerical approximate solutions. To do so, we determine the structure of the Jacobian matrix of the system and use it to deduce its eigenstructure. We combine the use of characteristic-based numerical fluxes with a second-order implicit-explicit scheme proposed in the cited reference and perform some numerical experiments with Tóth's adsorption isotherms to demonstrate that the characteristic-based schemes produce accurate numerical solutions with no oscillations, even when steep gradients appear in the solutions.

Analysis of the boundary conditions for the ultraweak-local discontinuous Galerkin method of time-dependent linear fourth-order problems

In this paper, we study the ultraweak-local discontinuous Galerkin (UWLDG) method for time-dependent linear fourth-order problems with four types of boundary conditions. In one dimension and two dimensions, stability and optimal error estimates of order k + 1 k+1 are derived for the UWLDG scheme with polynomials of degree at most k k ( k ≥ 1 k\ge 1 ) for solving initial-boundary value problems. The main difficulties are the design of suitable penalty terms at the boundary for numerical fluxes and the construction of projections. More precisely, in two dimensions with the Dirichlet boundary condition, an elaborate projection of the exact boundary condition is proposed as the boundary flux, which, in combination with some proper penalty terms, leads to the stability and optimal error estimates. For other three types of boundary conditions, optimal error estimates can also be proved for fluxes without any penalty terms when special projections are designed to match different boundary conditions. Numerical experiments are presented to confirm the sharpness of theoretical results.

A fragile points method with a numerical-flux-based interface debonding model to simulate the delamination migration in composite laminates

Predicting damage in composite laminates is essential for a better utilization of the material, but it is challenging especially when the inter- and intra-ply damages are coupled. In this paper, a Fragile Points Method (FPM) with a numerical-flux-based interface debonding model is proposed to simulate the laminates’ damage. The FPM is a novel meshless method using the discontinuous Galerkin weak formulations and the point-based discontinuous polynomial trial functions. In this work, an orthotropic interface debonding model has been developed based on the numerical fluxes to simulate the mixed-mode inter- and intra-ply damage. For any interior interface in an FPM model, the interface damage is triggered when the damage criterion is satisfied and the interface is separated when the fracture energy is reached. Due to the discontinuous characteristics of the FPM, the inter- and intra-ply damage modes are introduced explicitly, and thus the coupling between them can be naturally captured. In this paper, the proposed FPM is verified by four benchmark examples for Mode I, Mode II and Mixed Mode delamination. Then the proposed method is applied to simulate the delamination migration, used as an example for the inter- and intra-ply damage. The FPM results show that this method is able to provide reliable prediction on the load–displacement curves as well as the damage and crack patterns, and the mechanism governing the delamination migration is also discussed based on the FPM results.

Improved well-balanced AWENO schemes with hydrostatic reconstruction for the Euler equations under gravitational fields

The Euler equations under gravitational fields allow the hydrostatic equilibrium states, which requires that the numerical scheme of the system should also have this characteristic. In our previous work, a well-balanced finite difference conservative AWENO scheme has been constructed to preserve the isothermal equilibrium state accurately (Fu et al., Appl. Numer. Math., 180:1-15, 2022). However, the scheme fails to maintain more complicated hydrostatic equilibrium states, such as the isentropic equilibrium state. In this study, the improved well-balanced AWENO schemes with hydrostatic reconstruction are designed to maintain the general equilibrium states numerically. The key technologies to achieve the well-balanced property are rewriting the source term into an appropriate form based on the equilibrium equation, acting on the conserved variables with an affine-invariant WENO or linearized WENO interpolation operator, employing hydrostatic reconstruction methods to create auxiliary conserved variables, and modifying the numerical flux with auxiliary conserved variables and higher-order conserved terms based on pressure. The proposed well-balanced schemes maintain the equilibrium solution near the machine precision and capture small perturbations of the equilibrium state without spurious oscillations even at the coarse mesh. Furthermore, they can be applicable to the general hydrostatic equilibrium more than the isothermal equilibrium state. Several numerical experiments are tested to show that the developed well-balanced schemes perform well in maintaining equilibrium states and capturing small perturbations.

Penalty‐free discontinuous Galerkin method

AbstractIn this article, we present a new high‐order discontinuous Galerkin (DG) method, in which neither a penalty parameter nor a stabilization parameter is needed. We refer to this method as penalty‐free DG. In this method, the trial and test functions belong to the broken Sobolev space, in which the functions are in general discontinuous on the mesh skeleton and do not meet the Dirichlet boundary conditions. However, a subset can be distinguished in this space, where the functions are continuous and satisfy the Dirichlet boundary conditions, and this subset is called admissible. The trial solution is chosen to lie in an augmented admissible subset, in which a small violation of the continuity condition is permitted. This subset is constructed by applying special augmented constraints to the linear combination of finite element basis functions. In this approach, all the advantages of the DG method are retained without the necessity of using stability parameters or numerical fluxes. Several benchmark problems in two dimensions (Poisson equation, linear elasticity, hyperelasticity, and biharmonic equation) on polygonal (triangles, quadrilateral, and weakly convex polygons) meshes as well as a three‐dimensional Poisson problem on hexahedral meshes are considered. Numerical results are presented that affirm the sound accuracy and optimal convergence of the method in the norm and the energy seminorm.

An efficient second-order cell-centered Lagrangian discontinuous Galerkin method for two-dimensional elastic-plastic flows

An efficient second-order cell-centered Lagrangian discontinuous Galerkin (DG) method for solving two-dimensional (2D) elastic-plastic flows with the hypo-elastic constitutive model and von Mises yield condition is presented. First, starting from the governing equations of conserved quantities in the Euler framework, the integral weak formulation of them in the Lagrangian framework is derived. Next, the DG method is used for spatial discretization of both the weak formulation of conserved quantities and the evolution equation of deviatoric stress tensor. The Taylor basis functions defined in the reference coordinates provide the piecewise polynomial expansion of the variables, including the conserved quantities and the deviatoric stress tensor. The vertex velocities and Cauchy stress tensor on the edges are computed using a nodal solver equipped with a variant of Li's new Harten-Lax-van Leer-contact approximate Riemann solver [Li et al., “An HLLC-type approximate Riemann solver for two-dimensional elastic-perfectly plastic model,” J. Comput. Phys. 448, 110675 (2022)], in which the longitudinal wave velocity in the plastic state is modified. Then the vertex velocities and Cauchy stress tensor on the edges are used to compute numerical fluxes. A second-order total variation diminishing Runge–Kutta scheme is used for time discretization of both the governing equations of conserved quantities and the evolution equation of deviatoric stress tensor. After solving the evolution equation of deviatoric stress tensor, a radial return algorithm is performed at the Gauss points of each element according to the von Mises yield condition. And then the coefficients of the DG expansion for the deviatoric stress tensor on each element are modified by a least squares procedure using the deviatoric stress tensors at these Gauss points. To achieve second-order accuracy, the least squares procedure is used for piecewise linear reconstruction of conserved quantities and the deviatoric stress tensor, and the Barth–Jespersen limiter is used to suppress the nonphysical numerical oscillation near the discontinuities. After that, the coefficients of the DG expansion are modified through L2 projection using the reconstructed polynomials. Finally, a second-order cell-centered Lagrangian DG scheme is established. Several tests demonstrate that the new scheme achieves second-order accuracy with good robustness, and that the DG method of updating the deviatoric stress tensor has comparable accuracy and much higher efficiency with mesh refinement compared with previous works.

Positivity-preserving and entropy-bounded discontinuous Galerkin method for the chemically reacting, compressible Euler equations. Part II: The multidimensional case

In this second part of our two-part paper, we extend to multiple spatial dimensions the one-dimensional, fully conservative, positivity-preserving, and entropy-bounded discontinuous Galerkin scheme developed in the first part for the chemically reacting Euler equations. Our primary objective is to enable robust and accurate solutions to complex reacting-flow problems using the high-order discontinuous Galerkin method without requiring extremely high resolution. Variable thermodynamics and detailed chemistry are considered. Our multidimensional framework can be regarded as a further generalization of similar positivity-preserving and/or entropy-bounded discontinuous Galerkin schemes in the literature. In particular, the proposed formulation is compatible with curved elements of arbitrary shape, a variety of numerical flux functions, general quadrature rules with positive weights, and mixtures of thermally perfect gases. Preservation of pressure equilibrium between adjacent elements, especially crucial in simulations of multicomponent flows, is discussed. Complex detonation waves in two and three dimensions are accurately computed using high-order polynomials. Enforcement of an entropy bound, as opposed to solely the positivity property, is found to significantly improve stability. Mass, total energy, and atomic elements are shown to be discretely conserved.

High order conservative Lagrangian schemes for two-dimensional radiation hydrodynamics equations in the equilibrium-diffusion limit

Radiation hydrodynamics equations (RHE) refer to the study of how interactions between radiation and matter influence thermodynamic states and dynamic flow, which has been widely applied to high temperature hydrodynamics, such as inertial confinement fusion (ICF) and astrophysical gaseous stars. Solving RHE accurately and robustly even under the equilibrium diffusion approximation is a challenging task. To address this, we develop two types of high order conservative Lagrangian schemes for RHE in the equilibrium-diffusion limit for the two dimensional case on the Lagrangian moving mesh. Based on the multi-resolution WENO reconstruction for the spatial discretization and strong stability preserving Runge-Kutta (SSP-RK) time discretization, we first develop an explicit Lagrangian scheme with the HLLC numerical flux to achieve high order accuracy in space and time. We also discuss the positivity-preserving property of the high order explicit Lagrangian scheme. To overcome the severe time step restriction arising from the nonlinear radiation diffusion term in the explicit scheme, we further present a high order explicit-implicit-null (EIN) Lagrangian scheme. By adding a sufficiently large linear diffusion term on both sides of the scheme, we treat the complicated nonlinear parts explicitly and efficiently, and treat the added linear diffusion term on the right-hand side implicitly with a relaxed time step restriction. According to our numerical experiments, these two types of Lagrangian schemes are high order accurate, conservative and can capture the interfaces automatically. Additionally, the explicit scheme is found to be non-oscillatory and can preserve positivity while maintaining the original high order accuracy.

An Entropy Stable Discontinuous Galerkin Method for the Two-Layer Shallow Water Equations on Curvilinear Meshes

We present an entropy stable nodal discontinuous Galerkin spectral element method (DGSEM) for the two-layer shallow water equations on two dimensional curvilinear meshes. We mimic the continuous entropy analysis on the semi-discrete level with the DGSEM constructed on Legendre–Gauss–Lobatto (LGL) nodes. The use of LGL nodes endows the collocated nodal DGSEM with the summation-by-parts property that is key in the discrete analysis. The approximation exploits an equivalent flux differencing formulation for the volume contributions, which generate an entropy conservative split-form of the governing equations. A specific combination of a numerical surface flux and discretization of the nonconservative terms is then applied to obtain a high-order path-conservative scheme that is entropy conservative. Furthermore, we find that this combination yields an analogous discretization for the pressure and nonconservative terms such that the numerical method is well-balanced for discontinuous bathymetry on curvilinear domains. Dissipation is added at the interfaces to create an entropy stable approximation that satisfies the second law of thermodynamics in the discrete case, while maintaining the well-balanced property. We conclude with verification of the theoretical findings through numerical tests and demonstrate results about convergence, entropy stability and well-balancedness of the scheme.

High‐order gas kinetic flux solver for viscous compressible flow simulations

AbstractAlthough the gas kinetic schemes (GKS) have emerged as one of the powerful tools for simulating compressible flows, they exhibit several shortcomings. Since the local solution of continuous Boltzmann equation with the Maxwellian distribution function is used to calculate the numerical fluxes at the cell interface, the flux expression in GKS is usually more complicated. In this paper, a high‐order simplified gas kinetic flux solver (GKFS) is presented for simulating two‐dimensional compressible flows. Circular function‐based GKFS (C‐GKFS), which simplifies the Maxwellian distribution function into the circular function, combined with an improved weighted essentially non‐oscillatory (WENO‐Z) scheme is applied to capture more details of the flow fields with fewer grids. As a result, a simple high‐order accurate C‐GKFS is obtained, which improves the computing efficiency and reduce its complexity to facilitate the practical application of engineering. A series of benchmark‐test problems are simulated and good agreement can be obtained compared with the references, which demonstrate that the high‐order C‐GKFS can achieve the desired accuracy.

Numerical simulation of water jet load induced by spherical bubble collapse under underwater explosion

In this study, a numerical method combining the HLLC solver and FV-WENO scheme is proposed to accurately simulate the behavior of multiphase compressible fluids and investigate the formation process of water jets. The HLLC solver approximates the conservation equations of the fluid by using numerical fluxes from the approximate Riemann problem, while the FV-WENO scheme improves the accuracy of the numerical solution by employing the weighted essentially non-oscillatory (WENO) reconstruction method. The combination of these methods allows for a more accurate simulation of the behavior of multiphase compressible fluids, leading to an accurate depiction of the formation process of water jets. By simulating the interaction between underwater shock waves from explosions and single or multiple spherical bubbles, we study the generation, development, and collapse of water jets during the bubble collapse process. By revealing the characteristics of the load during the formation of water jets and the interaction mechanisms between bubbles, our study provides support for further research on water jet phenomena. The formation process of water jets is of great significance for understanding the propagation of underwater shock waves and structural damage. Through this study, we delve into the mechanism of water jet formation, providing a new perspective on the propagation of underwater shock waves and structural damage. Our research results demonstrate that the application of the HLLC solver and FVWENO scheme in water jet simulation effectively improves the accuracy and precision of the numerical solution.

Grad's 13 moments-based gas kinetic flux solver with triangle unstructured meshes for simulating continuum and rarefied fluid flows

A numerical framework based on the gas kinetic flux solver (GKFS) with unstructured meshes of triangle cells has been proposed and assessed in this work. Like conventional GKFS, the present scheme only updates the cell-average macroscopic conservative variables in time, thereby saving computational costs. At each cell interface, the numerical normal flux is computed with the GKFS scheme proposed by Liu et al. [J. Comput. Phys. 415, 109548 (2020)], in which the distribution functions at the surrounding points are reconstructed with the Grad's 13 moments (G13). Additionally, a solid wall boundary condition treatment has been proposed in the present G13-GKFS framework, which allows one to model from purely specular to perfectly thermalized solid walls. The present framework has been validated with the steady cylindrical Couette flow, the lid-driven cavity flow, the unsteady Rayleigh flow, and the rarefied flow around the NACA0012 airfoil test-cases, where good agreements are found with references for a quite wide range of flow regimes, from continuum to transitional flow regimes.

Application of higher-order FV-WENO scheme to the interaction between shock wave and bubble

The high-order finite volume-WENO (Weighted Essentially Non-Oscillatory) scheme combines the finite volume method with the WENO method, allowing for high-order accuracy and accurate simulation of complex physical phenomena. It has advantages in handling shock waves and bubble interactions. The idea of this method is to discretize the physical equations into a set of conservation equations and use the WENO method to calculate the numerical fluxes on each finite volume cell. This approach is effective in dealing with problems such as shock wave propagation, bubble deformation, and evolution. In fluid dynamics simulation and research, the high-order FV-WENO scheme has significant application prospects. It can provide accurate numerical solutions and simulate complex physical phenomena, making it widely applicable in scientific research and engineering. In this study, we simulated the interactions between shock waves and single or double bubbles, obtaining the complete process of bubble collapse with clear bubble interfaces and strong program stability.

A variant of improved discrete velocity method for efficient simulation of flows in entire Knudsen number regimes

In this paper, a variant of the improved discrete velocity method (VIDVM) is proposed for flows in the whole Knudsen number regimes. This method retains the advantage of the improved discrete velocity method (IDVM), which calculates numerical fluxes through a self-adaptive strategy by combining the microscopic reconstruction and the macroscopic reconstruction. Like the IDVM, the microscopic reconstruction for VIDVM is also based on the collisionless Boltzmann solver. However, different from IDVM, the macroscopic reconstruction for VIDVM is established on the Euler solver instead of the Navier–Stokes solver. Considering that the Euler solver merely computes the inviscid fluxes while the Navier–Stokes solver additionally calculates the viscous fluxes, the present method could be more efficient than IDVM. To validate the accuracy and efficiency of the present scheme, some benchmark cases from the continuum regime to the free molecular regime are conducted. Results reveal that the present scheme can predict the flow as well as IDVM, but the present solver is more efficient than IDVM.

Matric Flux Potential in Time Independent Infiltration Problems from a Single Triangular and a Trapezoidal Irrigation Channel

In this paper, steady infiltration problems into a homogeneous soil from a single triangular and trapezoidal irrigation channel are considered. The governing equation is Richard's equation that represents the movement of water in unsaturated soil. It is a non-linear equation and can be solved by linearizing to become a modified Helmholtz equation. Dual Reciprocity Boundary Element Methods (DRBEM) are used in this study to numerically solve the modified Helmholtz equation. Therefore, by using a provided solution, the numerical Matric Flux Potential (MFP) is calculated. This method was applied to the homogeneous soil problem of stationer infiltration from triangular and trapezoidal single irrigation. Both numerical solutions were compared. The result show that the MFP value from the triangular irrigation is higher than the trapezoidal irrigation. This indicates that content water from the triangular irrigation channel is higher than the trapezoidal irrigation channel.

Numerical model for determining the effective heat capacity of macroencapsulated PCM for building applications

This paper presents a finite difference model of macroencapsulated PCM panels coupled with the genetic algorithm for the determination of effective heat capacity of whole panels via inverse method. This provides an accurate characterization of the thermal properties of macroencapsulated PCMs for building envelope applications. A novel definition of the effective heat capacity is proposed based on the superimposition of two Gaussian curves, applicable to any PCM whose phase transition is characterized by a single peak. Three PCMs were tested, subjected to temperature variation rates typically experienced in building envelopes: 0.5 °C/h and 1 °C/h. Surface temperature and heat flux were measured and used in the inverse method procedure. The developed model is accurate, as numerical results greatly agree with the experiments: the root mean square difference between the experimental and numerical heat fluxes ranged between 0.543 and 1.246 W/m2. Significant differences in the effective heat capacity were found between the whole macrocapsule and small quantities of PCM (specified in the datasheets). The effective heat capacity specified in the datasheets is sensibly greater than that of the whole macrocapsules determined through the inverse method: the specific heat in the solid phase was up to 107.39 % higher in the datasheet values, the specific heat in the liquid phase up to 184.04 %, and the peak effective heat capacity, between 18.28 % and 164.11 %. The same happened to the enthalpy: datasheet values were 61.24 % – 175.55 % greater than inverse method results. This proves that latent heat is overestimated if small quantities of PCM are analyzed, and not the whole panels. The scale effect was assessed by comparing two capsules with the same material, but with different quantities of PCM: 0.5 kg and 1 kg. A greater mass of PCM over the total mass of the capsule implies a different relationship between the effective heat capacity and temperature, with higher peak heat capacity. The capsule with 1 kg of PCM showed a peak effective heat capacity up to 30.65 % greater than that of the panel with 0.5 kg of PCM. Thus, adequate modeling in building applications requires characterization of whole macroencapsulated PCMs. The determination of the relationship between temperature and effective heat capacity of macroencapsulated PCMs presented in this work could easily be incorporated into other simulation software, facilitating the assessment of adaptive envelopes with PCM macrocapsules.