An adaptive wavelet-based collocation method for solving multiscale problems in continuum mechanics

Computational multiscale methods are highly sophisticated numerical approaches to predict the constitutive response of heterogeneous materials from their underlying microstructures. However, the quality of the prediction intrinsically relies on an accurate representation of the microscale morphology and its individual constituents, which makes these formulations computationally demanding. Against this background, the applicability of an adaptive wavelet-based collocation approach is studied in this contribution. It is shown that the Hill–Mandel energy equivalence condition can naturally be accounted for in the wavelet basis, (discrete) wavelet-based scale-bridging relations are derived, and a wavelet-based mapping algorithm for internal variables is proposed. The characteristic properties of the formulation are then discussed by an in-depth analysis of elementary one-dimensional problems in multiscale mechanics. In particular, the microscale fields and their macroscopic analogues are studied for microstructures that feature material interfaces and material interphases. Analytical solutions are provided to assess the accuracy of the simulation results.


Introduction
Multiscale simulations techniques are used to predict the constitutive response of heterogeneous materials based on a detailed understanding of physical processes that take place at a lower scale.In computational multiscale methods, the lower scale information is incorporated into the model by geometrically resolving individual grains, phases and microscale imperfections such as voids or microcracks in a representative volume element (RVE).In doing so, sophisticated models that have been developed at the level of individual phases may be accounted for.By establishing scale-bridging relations, this detailed representation of the underlying microstructure is then used to derive the effective constitutive response at each macroscale material point, resulting in an hierarchical simulation approach.
First-order computational multiscale methods are meanwhile well-established to study multi-scale multi-physics problems in continuum mechanics [1].Over the past years, intense research efforts focused for instance on purely mechanical [2][3][4][5][6], thermo-mechanical [7][8][9][10][11], and electromechanical problems of electro-active solids [12][13][14] and of electrical conductors [15,16].All of these approaches have in common that they intrinsically rely on an accurate representation of the underlying microstructure, both in terms of its morphology and the constitutive response of individual constituents.For this reason, the repetitive solution of microscale boundary value problems is associated with significant computational costs and memory requirements, in particular in multi-physics applications or when non-linear constitutive relations are to be solved [17].Against this background, various approaches (each with their own merits and limitations) have been proposed to reduce the computational complexity of microscale boundary value problems-i.e. the number of degrees of freedom and the number material point eval-uations.Fast Fourier Transformation (FFT)-based spectral solvers [18][19][20][21] are for instance applied to take benefit of their high computational efficiency.However, convergence issues are encountered for increasing phase contrasts [20], scalability problems occur [1,22,23], and classic FFT-based spectral solvers are associated with a uniform spatial resolution [19] which is not optimal if small-scale features occur only in some (localised) parts of the solution domain.Proper Orthogonal Decomposition-based Reduced Order Models (ROMs) [24] rely on an offline phase in which (deformation) modes are extracted from a selection of snapshots to significantly reduce the problem dimension in the online phase.However, a sufficiently large number of snapshots is required which complicates the application of ROMs when geometric and material non-linearities apply.In an attempt to reduce the numerical costs associated with the evaluation of microscale material models, Empirical Interpolation Methods [25][26][27] or Cubature Methods [28,29] make use of stress modes and a reduced set of quadrature points with optimised weights, respectively.The numerical error that is induced by these approximations has, however, been shown to significantly increase with increasing deformation for non-linear historydependent material behaviour, see [30].
Against this background, this contribution aims to assess the applicability of adaptive wavelet-based approaches to reduce the computational complexity of microscale boundary value problems.Wavelet analysis, often also referred to as the numerical microscope [31], dates back to the seminal work by Morelet [32].It is based on the fundamental concept of multiresolution analysis and has been shown to be particularly useful in the study of physics problems where different resolutions are required in different regions of space [31,33,34], see Remark 1.In particular, by using only a high resolution in regions where it is required, the dimensionality of the problem (and the number of material point evaluations) can be reduced significantly.Moreover, the resolution in different regions can systematically be adapted by adding/removing higher-level wavelets in the expansion of the respective field variable [33,34].For these reasons, wavelet-based solvers for partial differential equations have been in the focus of intense research in the past decades, see for instance [31] for a detailed review.Among the approaches that have been studied are wavelet collocation methods, wavelet Galerkin methods, and wavelet finite element methods.Wavelet collocation methods [35][36][37], in which the primary field is expanded in a wavelet basis, have their merits in dealing with nonlinear problems.Wavelet Galerkin methods [38][39][40], in which trial and test functions are approximated in a wavelet basis, are particularly strong when efficient algorithms for the exact calculation of connection coefficients (i.e.integrals of products of wavelets and their derivatives) can be used.Wavelet finite element methods [41] (as a special case of wavelet Galerkin methods), in which wavelet theory is applied to construct finite element shape functions, make use of the multiresolution property to systematically increase the analysis precision while maintaining the original grid.Focusing on continuum mechanics applications: A wavelet collocation method is used and wavelet analysis is applied in [42] to detect regions of a sample that are damaged or undergo plastic deformation.A wavelet Galerkin approach is adopted in [43] to solve structural optimisation problems for elastically deforming continua.A wavelet finite element method for the efficient solution of continuum mechanics problems, making use of the close relation between hierarchical finite elements and hat-wavelets, is studied in [44] for purely elastic problems and extended to inelastic material behaviour in [45].Moreover, addressing the drawbacks of classic Empirical Interpolation and Cubature Methods, an adaptive wavelet-based reduced order model for the solution of microscale boundary value problems is developed in [46,47].The principal idea of this approach is to project occurring integrands onto a wavelet basis and to significantly reduce the number of material point evaluations by dynamically adapting the set of integration points.Overall there are but a few works available that focus on the development and application of tailored wavelet-based approaches for continuum mechanics problems that include non-linear history-dependent material behaviour.In particular, detailed studies of wavelet-based solution approaches in the context of computational homogenisation-based multiscale mechanics are, to the authors best knowledge, not available despite the multiresolution property of wavelets which seems promising to properly resolve (possibly evolving) microscale features in representative volume elements.
Against this background, the fundamental characteristics of an adaptive wavelet collocation approach applied to typical problems that occur in multiscale mechanics are studied in this contribution.In particular, focus lies on multiphase materials (featuring material interfaces and interphases) that undergo reversible and irreversible deformations.To gain a basic understanding of the capability of wavelet-based approaches in capturing microscale fields and associated macroscale quantities, elementary (but sufficiently complex) one-dimensional problems are studied for which analytical solutions can be furnished.By carrying out an in-depth study of elementary problems, the present contribution thus provides a basis for future developments and contributes to the (long-term) research question: Can adaptive wavelet discretisations of field variables be used to improve the overall computational efficiency of multiscale schemes?In this regard, the following intrinsic properties of the proposed wavelet-based formulation are key: 1. Constitutive relations: The proposed formulation is capable of accounting for inelastic, non-linear material behaviour-i.e.general constitutive models σ (ε, •) relat-ing stresses σ to strains ε and, possibly, to a set of internal variables • can be accounted for.2. Adaptivity: Discretisations based on interpolating wavelet families naturally give rise to grid-refinement criteria that make use of wavelet-coefficients.

State variables:
The dynamic adaption of the calculation grid requires state variables to be mapped between different grids, which can efficiently be done by using the interpolating property of the chosen wavelet family and by carrying out wavelet synthesis and wavelet analysis transformations.4. Scale transition: The vanishing moment property of the chosen lifted-interpolating wavelet family leads to an efficient computational homogenisation scheme relating microscale to macroscale quantities. 5. Hill-Mandel condition: Periodicity conditions on microscale fields that are associated with the Hill-Mandel energy equivalence condition can conveniently be incorporated into the wavelet basis.6. Meshfree methods: Wavelet-based collocation approaches are intrinsically meshfree which facilitates the analysis of multiphase microstructures with complex morphologies.
The article is organised as follow: the fundamentals of function approximation by means of wavelets, with particular focus lying on lifted interpolating wavelets, are briefly summarised in Sect. 2. These serve as a basis to approach the computational multiscale problem posed in Sect.3.More specifically speaking, an adaptive wavelet-based collocation methods is proposed in Sect. 4 and its applicability to multiscale problems that occur in continuum mechanics is studied in detail in Sect. 5 by means of comparison with analytical solutions.The findings are summarised and concluding remarks are drawn in Sect.6.

Remark 1 (Multiscale approaches)
In the continuum mechanics and the wavelet community, the term "multiscale approach" is used with a slightly different meaning.Wavelets operate on different resolution levels and naturally give rise to nested grids in numerical approaches.These resolution levels are frequently referred to as scales, and eventually allow for different resolutions to be adopted in different regions of space.
The classic notion of multiscale methods adopted in the continuum mechanics community is slightly different, with the focus being on so-called coarse graining methods.In these approaches classic material models at the macroscale are replaced by additional boundary value problems on unit cells that take detailed information of the underlying material microstructure into account and scale-bridging relations are established to relate processes at the different length scales considered.
The present contribution focuses on the wavelet-based solution of microscale boundary value problems, i.e. on the proper resolution of microscale features such as grain boundaries and localisation zones.Moreover, the wavelet-based microscale solution approach is additionally embedded into a multiscale approach in the sense of continuum mechanics by establishing closed-form relations for the effective macroscale quantities and, in particular, for the algorithmic consistent macroscale tangent stiffness tensor.

Notation
Tensor notation is adopted throughout this contribution.In particular, with α, β, γ and δ denoting arbitrary first-order tensors and with ⊗ denoting the standard dyadic product, single tensor contractions are understood in the sense gradient and (left-)divergence operators are denoted as ∇ and ∇•, and superscripts • t and • sym indicate the transpose and the symmetric part of a second order tensor, respectively.

Wavelet theory
The fundamentals of function approximation in a onedimensional setting by means of biorthogonal wavelet families are briefly summarised in Sect.2.1 and specified for the particular class of lifted-interpolating wavelets in Sects.2.2 and 2.3.These serve as the basis for the development of the adaptive wavelet-based collocation approach in Sect. 4. A more thorough introduction to the mathematical theory of wavelets is for instance provided in [39,48,49], the interpolating Deslauriers-Dubuc wavelet family is covered in textbooks [33,34], and the suitability of lifted interpolating wavelets for the solution of partial differential equations is discussed in [35].

Fudamentals of biorthogonal wavelet families
Wavelet theories are based on the fundamental concept of multiresolution analysis.A multiresolution analysis of a function space formally consists of a sequence of approximation spaces V j : j ∈ Z and scaling functions φ j k : k ∈ Z spanning V j .Amongst others, the approximation spaces are required to fulfil the nesting property V j ⊂ V j+1 and to contain translates and dilates of functions Taking into account the nesting property, wavelet spaces W j are defined as the complement of V j in V j+1 , i.e.V j+1 = V j ⊕ W j , and wavelets are introduced as functions ψ in a level-1 discretisation.Likewise, the grid point at x = 0.5 mm is associated with wavelet ψ 0 0 and scaling function φ 0 1 , respectively W j .This is exemplified for hat-wavelets in Fig. 1.Likewise, dual scaling functions φ j k and dual wavelets ψ j k of the biorthogonal wavelet families under consideration may be introduced.
With the sequence of nested approximation spaces at hand, a function f (x) may be approximated at different resolution levels j in terms of the scaling function coefficient s j k and wavelet coefficients d j k via ( In this regard, a representation in terms of both scaling functions and wavelets is referred to as a wavelet representation, whereas a representation that is solely based on scaling functions is referred to as a scaling function representation.In the wavelet representation of function f j (x), the wavelet coefficients contain the additional detail information that is required when passing from an approximation of f (x) at resolution level j − 1 to an approximation at resolution level j.Vice versa, by dropping wavelet coefficients smaller than a preset tolerance d , function f j (x) may be represented in a compressed form with only a minor loss in accuracy.In other words, high levels of resolution are only taken into account where needed.The scaling functions φ j k and wavelets ψ i k in (2), that occur at different resolution levels j and that are associated with different translation parameters k, can be expressed as dilates and translates of mother scaling functions and of mother wavelets respectively.Moreover, the nesting property of the sequence of approximation spaces gives rise to the refinement relations for the scaling functions and for the wavelets The coefficients h l and g l that occur in the refinement relations are characteristic for a specific biorthogonal wavelet family and define the latter together with their dual analogues hl and gl .In particular, hl and gl can be interpreted in terms of low-pass and high-pass filters in the forward wavelet transform (wavelet analysis) whereas h l and g l take the interpretation of synthesis filters in the backward wavelet transform (wavelet synthesis) Moreover, by collecting all scaling function coefficients at level j in vector s j and all wavelet coefficients in vector d j it is observed that ( 7) and ( 8) may be expressed in the form of matrix-vector multiplications as follows with F j+1 j and B j j+1 denoting forward and backward transformation matrices.For a prescribed maximal grid level n G one thus obtains the transformation matrices and The transformation matrices B n G 0 and F 0 n G relate level-n G wavelet representations to level-n G scaling function representations.Note that these matrices can be pre-calculated for a given grid.

Lifted interpolating wavelets
In this contribution, lifted interpolating wavelets will be used as basis functions.Their properties and the specific filter coefficients {h l , g l , hl , gl } are controlled by the two parameters N and Ñ as explained in the following.
In order to construct this particular wavelet class, a set of nested dyadic grids G j is introduced, see Fig. 2. At each resolution level j, grid points x j k = 2 − j k with k ∈ Z are taken into account.Regarding a scaling function representation, each grid point x j k is associated with a scaling function φ j k .Likewise, taking into account wavelet representation (1), level zero grid points x 0 k are associated with level zero scaling functions φ 0 k , whereas grid points x j+1 2k+1 for j ≥ 0 are associated with wavelets ψ j k on higher resolution levels.The construction of the lifted interpolating wavelet family as outlined in, e.g., [35] results (for regular periodic grids) in a forward wavelet transform of the form and a backward wavelet transform of the form with w k+l denoting interpolation coefficients from even grid points x j+1 2k+2l to odd grid points x j+1 2k+1 , and with wk+l denoting interpolation coefficients from odd grid points x j+1 2k+2l+1 to even grid points x j+1 2k that occur in the derivation.For different numbers of support points 2 N , respectively Table 1 Interpolation coefficients w l and wl in the construction of (lifted) interpolating wavelets according to (11) and ( 12) observed that ( 11) and ( 12) stipulate filter coefficients h, g, h, g.The parameters N and Ñ thus characterise the lifted interpolating wavelet family, with Ñ controlling the lifting procedure.In particular, it is noteworthy that nonlifted interpolating wavelets are recovered for Ñ = 0 and that scaling function φ is not affected by the lifting procedure.
Focusing in more detail on the properties of the scaling function, it is observed that φ is interpolating, that supp (φ) = [−2N + 1, 2N − 1], holds, and that linear combinations of φ j k exactly reproduce polynomials up to degree 2N − 1, cf. [35].Moreover, N and Ñ can be shown to control the number of vanishing moments of the scaling function and wavelet according to supp(φ) see [35,50,51].The behaviour of the basis functions at the origin and hence the localisation behaviour in Fourier space are thus influenced by the choice of N and Ñ since holds for the Fourier transformed function This observation is important for the stability of the numerical scheme since significant aliasing effects may occur if the wavelet spaces are not well-separated, possibly leading to unstable or inaccurate results, see, e.g., the discussion in [33,35].For different values of N and Ñ , scaling functions and wavelets of the (lifted-)interpolating wavelet family under consideration are exemplarily depicted in Figs. 3, 4, 5 In order to gain a basic understanding of the capability of wavelet-based approaches in capturing microscale fields and associated macroscale quantities, discretisations with basis functions that result from different combinations of the parameters N and Ñ will be in the focus of Sect. 5.

Calculation of derivatives
Due to the interpolating property of the scaling functions, derivatives in physical space can conveniently be approximated in a level-n G scaling function representation by finite difference schemes.More specifically speaking, by interpret- Table 3 Filter coefficients of interpolating wavelets with ing the coefficients of the finite difference scheme in terms of derivative filters, one arrives at For two-and three-point forward difference schemes, the non-zero filter coefficients read two-point forward differences: three-point forward differences: where k denotes the grid spacing on level n G .Furthermore, by introducing filter matrices C n G and G n G 0 , (15) can equivalently be represented in terms of the following matrix-vector multiplications which are frequently used in the remainder of the manuscript to simplify the presentation.

Multiscale modelling: the RVE problem
The governing set of balance equations and constitutive relations is briefly recapitulated in Sect.3.1 before the fundamentals of computational multiscale methods are summarised in Sect.3.2.The solution of the resulting multiscale boundary value problem in a one-dimensional setting by means of a wavelet-based collocation approach will then be addressed Sect. 4.

Balance equations and constitutive relations
The mechanical problems that are in focus of this contribution are governed by the balance equation of linear momentum in its local form with mass density ρ (x), body-force density f (x) and displacement field u (x).The balance equation of angular momentum reduces to the symmetry condition of the Cauchy stress tensor that is accounted for in the specific form of the constitutive relations and in the computational homogenisation scheme to be discussed in Sect.3.2.By additionally assuming a functional dependence of the (volume-specific) Helmholtz free energy density function Ψ on the strain ε (x) = ∇ sym u (x), the temperature θ (x) and possibly several internal variables •, and by considering the Clausius-Planck form of the dissi-pation inequality with s denoting the mass-specific entropy density, thermodynamically consistent constitutive relations are derived.
In particular, (20) stipulates an energy-based definition of stresses and of thermodynamic driving forces * according to

Multiscale problem
In computational multiscale approaches, the constitutive relations at the macroscale are substituted by the homogenised solution of a microscale boundary value problem that enables to incorporate detailed information from the underlying material microstructure, cf.Fig. 7.In the following, subscripts "M" are used to indicate macroscale quantities.Quantities without subscripts refer to the microscale, i.e. no additional subscript is introduced in order to increase the readability of equations.By appealing to a separation of time-and length-scales, the microscale boundary value problem is subjected to the assumptions of quasi-statics and negligible body forces such that (18) reduces to The primary unknown of the microscale boundary value problem is the micro-fluctuation field ω (x) that characterises the deviation from an affine deformation given by the macroscale strain tensor ε M according to and that gives rise to the microscale strain field The macro-micro scale transition thus results from the influence of the macroscale strain tensor ε M on the microscale deformation field.Vice versa, the effective macroscopic stress tensor σ M is calculated as the volume average of its microscopic counterpart which can be shown to be fulfilled for periodic microfluctuation fields.By indicating quantities at opposing parts of the RVE-boundary ∂B with outward unit normal vector n by superscripts • − and • + , the periodicity condition can be cast in the form and implies The solution of the microscale boundary value problem and the treatment of scale-bridging relations when using a wavelet discretisation of field variables is in the focus of Sect. 4.

Adaptive wavelet collocation approach
Based on the fundamentals of computational multiscale methods discussed in Sect.3, an adaptive wavelet-based collocation approach is proposed in this section.At the outset of the developments, we restrict ourselves to the one-dimensional case which simplifies the presentation of main ideas and allows multiple analytical solutions to be derived and to be compared against in Sect. 5.In particular, a wavelet-discretised form of the balance equation is provided in Sect.Since the reduction of the general computational multiscale scheme discussed in Sect. 3 to a one-dimensional setting is straightforward, the governing set of equations is not repeated in this section.

Discretisation
In accordance with Sect.3.2, the balance equation ( 18) is enforced at n B col domain collocation points under the assumptions of negligible body forces and quasi-statics.The Dirichlet-type condition ω = 0 is taken into account at n ∂B ω col boundary collocation points, see Fig. 8a.By additionally introducing an auxiliary field σ (x) that is coupled to the constitutively determined stresses at microscale collocation points χ , one arrives at the discrete residual-type formulation where A denotes the assembly operator.Note that the set of 1 is adaptively controlled, as described in Sect.4.2, and forms a subset of the dyadic grid G n G introduced in Sect.2.1.Moreover, periodicity conditions (27) on the primary fields can naturally be accounted for in the wavelet-based collocation approach by wrapping around subscripts that are out of bounds in all filter-based operations.This is indicated in Fig. 8a by dashed lines that represent the periodic extension of domain B. In addition, boundary collocation points are indicated by ⊗.
In the light of the wavelet representation (29), it is observed that the evaluation of the residual vectors (28) requires the values of stresses σ , micro-fluctuations ω, stress gradients and strains at collocation points.These can be calculated by successive backward wavelet transforms according to (8) and by applying derivative filters according to (15) or by using pre-calculated sparse matrices B n G 0 and G n G 0 .In particular, by collecting the scaling function and wavelet coefficients that occur in (29) in vectors υ • and by denoting submatrices with • and •, the tangent stiffness contributions K * • = dR * /dυ • that are required in a gradientbased iterative solution of the non-linear system of equations (28) take the form

Adaptive grid refinement
In this contribution, a grid refinement algorithm similar to the one proposed in [35] is used.For simplicity, the algorithm will solely be based on the micro-fluctuation field and the same grid will be used for the micro-fluctuation and the micro-stress field.However, there is of course no necessity to do so and different grids or refinement strategies could be used for each field.
To update the calculation grid we proceed as follows: 1. Starting from an empty set C, grid points x 0 k that are associated with level-0 scaling functions φ 0 k are added to the active set of collocation points.8b.Whereas the former grid points are associated with spatial translates of the wavelet centred at x j k , activation of the latter grid points allows higher frequency components of the solution to be accounted for.Likewise, level-1 grid points are activated if | [s ω ] 0 k | > r and only neighbouring grid points on the same resolution level are taken into account if the maximum refinement level is reached.4. In a final step, additional grid points are added to C such that the perfect reconstruction criterion is fulfilled-i.e.such that the forward wavelet transform (11) can be carried out, cf.[35]. 5.If the set of active grid points C is identical with the one of a previous grid iteration step, the grid (and the corresponding numerical solution of the PDE) is accepted and the algorithm proceeds to the next load step.Otherwise, the microscale boundary value problem is solved again on the updated grid and the refinement algorithm is restarted.

Mapping of state variables
In the collocation approach taken, the material model is evaluated at each active grid point and an update in the state variables is carried out.To this end, the state variables from the end of the previous load step are required at each active grid point.However, in the adaptive approach taken the set of active grid points at the end of the previous load step is in general different from the one in the current load step, which requires state variables to be mapped between different calculation grids.
Since the set of active grid points fulfils the perfect reconstruction criterion and since the state variables are available at each active grid point, the state variable profile at the end of each load step can be expanded in the chosen wavelet basis by successive forward wavelet transformations according to (11).Vice versa, by making use of the wavelet representation and by successively applying backward wavelet transformations according to (12), initial values for the state variables at grid points that were active in the previous load step can be retrieved and approximations for the initial values of states variables at grid points that were activated in the present load step can be calculated.

Discrete scale-bridging relations
By invoking (13), the homogenisation scheme (25) reduces in the considered one-dimensional setting to with l denoting the length of the representative volume element and with η 0 denoting the grid spacing at discretisation level 0, see Fig. 8a.Moreover, (31) serves as the basis for the derivation of the algorithmic consistent macroscale tangent stiffness operator which takes the form The evaluation of ( 32) requires sensitivities of level-0 scaling functions coefficients s σ 0 k with respect to prescribed macroscale strains ε M .By considering the total differential of the residual vectors defined in (28) from which the derivatives sought can be extracted.

Validation
The applicability of the proposed wavelet-based computational multiscale formulation is studied in this section by means of comparison with analytical solutions.To this end, a thorough study of elastically deforming microstructures that feature material interfaces and material interphases is carried out in Sect.5.1.The study is extended to elasto-plastic constitutive relations in Sect.5.2 with particular focus on the mapping of state variables between different calculation grids.From a physics point of view, the different material regions could for instance be thought of as plys in a simplified laminate problem that are arranged in a regular periodic pattern.

Reversible material response
This section concerns the validation of the proposed waveletbased computational homogenisation scheme by means of comparison with analytical solutions for a purely reversible material response.To this end, two sample boundary value problems are studied in a one-dimensional setting.In both cases, the material response is assumed to be governed by the Helmholtz energy density function which by evaluation of (21a) results in the linear-elastic constitutive relation More specifically speaking, periodic microstructures consisting of two materials with different elastic stiffness constants E as schematically depicted in Fig. 9 are considered in Sect.5.1.1 and Sect.5.1.2.Sharp transitions at the material interfaces (i.e.discontinuities in the elastic stiffness) are assumed first whereas smooth transition zones are resolved in terms of interphases of finite thickness next.From a numerical point of view, adaptive discretisations based on different interpolating wavelets are considered.In particular, a level-0 grid spacing η 0 = 0.05 mm is used in the discretisation of the periodic unit cells of length l = 1.0 mm and n G = 9 additional discretisation levels are taken into account.

Material interface
The one-dimensional unit cell under consideration is assumed to consist of two materials with different elastic properties according to Fig. 9a and with In particular, E = 210,000 N/mm 2 and α = 0.6 are assumed, and the geometric dimensions (x 1 = 0.2875 mm, x 2 = 0.75 mm, x 3 = 1.0 mm) are chosen such that the discontinuities in the elastic stiffness E occur at level-2 and level-0 grid points, respectively.The analytical solution for  the micro-displacement field reads and the one for the micro-fluctuation field ω ana follows by evaluating (23).Moreover, the micro-stresses σ that occur in (38) can be related to the prescribed macroscale strains ε M by taking into account the total elongation of the representative volume element, i.e.
Focusing on the simulation results of lifted interpolating wavelets with N = 1, Ñ = 1 and using a two-point forward difference scheme, it is observed that although the basis functions are, in principle, able to exactly reproduce the piecewise linear analytical solution, significant discretisation errors occur that are invariant of refinement tolerance r as shown in Fig. 10b.In particular, the first kink in the micro-fluctuation field is predicted at x = 0.3 mm, which corresponds to a level-0 grid point, instead of x = 0.2875 mm, which corresponds to a level-2 grid point, see Fig. 10a.This result can be explained by taking into account the computational grids in Fig. 11a-c which demonstrate that the wavelet coefficient-based refinement strategy fails for the chosen wavelets in the case of a piecewise linear solution profile since higher level grid points are not activated when starting from a level-1 discretisation.This elementary example thus shows that care needs to be taken when using lifted interpolating wavelets with N = 1, Ñ = 1 together with a two-point forward difference scheme in continuum mechanics applications.Note that higher-order finite differences are not trivially applicable in this case due to the limited continuity of the basis functions.
In contrast, taking into account discretisations based on lifted interpolating wavelets with N = 2, Ñ = 1 the computational grid is refined in the neighbourhood of the material interfaces where localised, high-frequency basis functions are required to approximate the solution, cf.Fig. 13 Relative errors in macroscale stresses σ M and in algorithmic consistent tangent stiffness E M = dσ M /dε M for a prescribed macroscale strain state, ε M = 0.005 ticular, it is observed that the maximal allowed discretisation level is activated for moderate values of refinement tolerance r .Moreover, it is noted that the maximal discretisation error does not significantly decrease further with decreasing values of refinement tolerance r , when the maximal discretisation level is reached in simulations using three-point forward differences.The discretisation error can, however, be further decreased by increasing the maximal number of discretisation levels as shown in Fig. 12.
In order to validate the computational homogenisation approach proposed in Sect.4.4, the numerically calculated effective macroscopic stresses and the associated algorithmic consistent tangent stiffness operators are compared against their analytical counterparts.The relative errors are depicted in Fig. 13 as a function of refinement tolerance r , with a significant reduction in the error with decreasing values of r being observed.

Material interphase
In contrast to the boundary value problem discussed in Sect.5.1.1,a smooth transition of material parameters at the material interfaces is assumed in this section according to Fig. 9b.In particular, the interfaces between the two different materials (E = 210,000 N/mm 2 , α = 0.6, x 1 = 0.2 mm, x 2 = 0.3 mm, x 3 = 0.7 mm, x 4 = 0.8 mm, x 5 = 1.0 mm) are resolved as interphases of finite thickness by making use of trigonometric functions, namely, Following the same lines as in Sect.5.1.1,an analytical solution for the micro-displacement field can be derived as a function of the micro-stress field, i.e.
and a relation between micro-stresses and (prescribed) macro-strains can be established Focusing on the simulation results for lifted interpolating wavelets with N = 1, Ñ = 1, first, it is observed that the continuously changing slope in the solution profile causes significant mesh refinement at the material interphases which is controlled by the refinement tolerance r .In particular, the relative error with regard to the analytical solution significantly reduces with decreasing values of r as shown  14a, b, whereby the number of active grid points in Fig. 15a-c significantly increases.In contrast, a notably smaller increase in the number of active grid points with decreasing values of r and a significantly smaller error with regard to the analytical solution are observed when using the same numerical differentiation scheme and setting the wavelet parameters to N = 2 and Ñ = 1, see Fig. 14c, d and Fig. 14d-f.Moreover, by comparing the simulation results for the same wavelet class but different forward difference schemes, a severe influence of the numerical differentiation scheme on the accuracy of the simulation results is revealed in Fig. 14d-f.
The error with regard to the analytical solution is observed to stagnate in Fig. 14d although grid points on level six and higher are not activated.This is in striking contrast to the observations in Sect.5.1.1where grid points at the maximum discretisation level were activated.However, by increasing the maximum number of discretisation levels n G , the maximum error can still be reduced in the present simulations as shown in Fig. 16.In particular, the decrease in maximum error from about 4.2 × 10 −4 to 2.1 × 10 −4 is in good accordance with the O (η)-accuracy of the finite difference scheme that is evaluated at the maximum discretisation level according to (15).Together with the results for the three-point forward difference scheme shown in Fig. 14f, these observations suggest a limitation resulting from the accuracy of the finite difference scheme.
Concluding the analysis of the simulation results, the convergence behaviour of the effective (homogenised) macroscale stresses and tangent stiffness operator is depicted in Fig. 17.In essence, a significant decrease in the relative error with decreasing refinement tolerance r is observed.

Irreversible material response
This section extends the analysis presented in Sect.5.1 to modelling of irreversible material responses that necessitates a mapping of state variables between different grids.As a sample problem, elasto-plastic relations with linear, isotropic hardening are considered.The governing evolution equations are briefly in Sect.5.2.1 and application of the proposed wavelet-based multiscale formulation to plastically deforming microstructures is studied in Sect.5.2.2 with particular focus lying on numerical errors induced by the variable mapping.

Constitutive relations
The elasto-plastic material behaviour is assumed to be governed by the free energy function and the yield function that feature the linear hardening modulus H and the initial yield limit q 0 .The energetic duals to the internal variables ε p and κ, that take the interpretation of plastic strain and accumulated plastic strain, follow from the evaluation of (21b) Moreover, by assuming an associative yield function, (44) stipulates evolution equations for the internal variables in terms of plastic multiplier λ p , namely, that are subjected to the Karush-Kuhn-Tucker conditions The system of evolution equations and constitutive relations, ( 44)- (47), is discretised in time by means of an implicit back-ward Euler scheme and a classic predictor-corrector approach is adopted to resolve the elastic and plastic deformations.Finally, by making use of ( 44)- (47) it is observed that the tangent stiffness tensor at each microscale material point reads for purely elastic loading, unloading or neutral loading whereas it takes the form for plastic loading.

Representative simulation results
In order study the properties of the proposed formulation for plastic processes, a periodic unit cell according to Fig. 9a with x 1 = 0.25 mm, x 2 = 0.75 and x 3 = 1.0 mm is taken into consideration.For simplicity, both material regions reveal the same elastic material properties and hardening moduli.However, their yield limits are assumed to differ by a factor α p .Under these assumptions, the analytical solution for the macroscale stresses σ M for monotonic loading to a maximum prescribed strain ε max M , followed by elastic unloading, reads For ε max M = 1.0%, α p = 0.8 and the specific choice of material parameters E = 210,000 N/mm 2 , H = 15,000 N/mm 2 , q 0 = 525N/mm 2 , (50) results in the piecewise linear macroscale stress-strain relation depicted in Fig. 18a.In particular, it is observed that the material parameters where chosen such that the kinks in the solution profile, which indicate the transition from a purely elastic to an elasto-plastic material response of a material region, occur at the strain-stress pairs 0.2%, 420N/mm 2 and 0.6%, 525N/mm 2 .
The numerical solution that was calculated by using lifted interpolating wavelets with N = 2, Ñ = 1, two-point forward differences, and a maximum number of refinement grids n G = 8 is additionally provided in Fig. 18a where constant load steps ε M = 0.025% were used.Focusing in more detail on the quality of the numerical solution, the relative error with respect to the analytical solution is provided in Fig. 18b for different values of refinement tolerance r .In particular, a characteristic in the error is observed in load step 10 for r = 10 −7 mm and r = 10 −8 mm, in load step 11 for r = 10 −8 mm.These are the load steps where a non-trivial state variable profile is for the first time and a detailed analysis of the numerical error by the state variable mapping indeed suggests that the sudden increase the σ M -error can be related to this mapping.More specifically speaking, the plastic strain ε p at active grid points is provided in Fig. 19a, b for r = 10 −8 mm.Whereas (perfect) piecewise constant profile is observed in load step 9 (-), overshoots in the vicinity of the material interfaces are observed after the mapping of the state variables from step 9 to 10 (--).The error thus induced is not corrected by the state variable update in step 10 (• • • ) since σ is below the initial critical yield stress in the material regions 0.0 mm ≤ x < 0.25 mm and 0.75 mm ≤ x < 1.0 mm.In contrast, focusing on load step 25, i.e. the first load step after the second elasto-plastic transition, the deviation from the piecewise constant state variable profile in load step 24 (-) is corrected by the state variable update in step 25 (• • • ).Accordingly, a significant drop in the σ M -error is observed.
The same observations hold for the calculations that were carried out with r = 10 −7 mm and r = 10 −9 mm, respectively.However, whereas the calculation grid is not adapted for r = 10 −8 mm ( r = 10 −9 mm) after load step 18 (load step 23) and a moderate overhead in microscale calculations is observable in Fig. 20, the significant overhead for r = 10 −7 mm suggests that the refinement algorithm struggles to adapt the calculation grid to the solution profile for "too high" values of r .
To put this into perspective, it is remarked that although the microscale boundary value problem needs to be solved multiple times to dynamically adjust the calculation grid, the number of active grid points that is directly related to the number of degrees of freedom and the number of material model evaluations is significantly reduced.For instance, comparing the maximum number of active grid points for refinement tolerance r = 10 −9 mm (i.e.168) with the number of grid points that would occur in a uniform level 8 discretisation (i.e.20×2 8 = 5120) yields a reduction of about 96.7%.This gain is even expected to increase enormously when considering more general 2D or 3D problems.

Conclusions
In this contribution, a wavelet collocation method for the numerical solution of multiscale problems that occur in continuum mechanics was presented.The formulation is capable of accounting for general non-linear constitutive relations featuring state variables and relies on a wavelet coefficient-based refinement criterion to dynamically adapt the calculation grid.Establishing a basis for future developments, particular focus was put on the accuracy and principal properties of the formulation in capturing microscale fields and effective macroscale quantities derived thereof.To this end, elementary linear and non-linear constitutive relations were considered and fundamental one-dimensional boundary value problems that occur in the simulation of multiphase microstructures were analysed.For each problem, an analytical solution was provided as a reference to assess the quality of the simulation results.
The refinement characteristics and error profiles induced by material interfaces and interphases were studied in a first step for purely elastic material behaviour.In a second step, these studies were extended to elasto-plastic material behaviour as a prototype for the treatment of constitutive relations featuring internal variables.In this regard, emphasis was placed on the additional numerical error induced by the mapping of state variables between different calculation grids when using a wavelet-discretisation.
With regard to the research question posed in the introduction and the principal properties of the wavelet-based formulation studied, the following concluding remarks can be drawn: 1. Constitutive relations: It was that non-linear, history-dependent material behaviour can conveniently be accounted for in the proposed wavelet-based collocation approach.2. Adaptivity: It was shown that the chosen coefficient-based refinement algorithm is capable of dynamically adapting the numerical grid to the solution profile.In the numerical grid was automatically refined close to material and interphases-i.e. in regions where significant changes in the solution profile are expected based on the phasecontrast-and close deformation-induced elasto-plastic transition zones.Hence, the hierarchical character of the wavelet-based approach and the naturally occurring local refinement characteristics allowed for an accurate representation of localised in the microscale fields.This is remarkable compared to, for instance, wellestablished FFT-based spectral approaches for which the treatment of localised deformations and weak discontinuities poses a severe challenge due to the reliance on a regular grid spacing.On the other hand, when considering for instance adaptive FE-based approaches, a close relation between hat-wavelets and hierarchical finite elements is observed, as outlined in the introduction.In this regard, the well-established concept of multiresolution analysis can be considered as a more sophisticated approach to derive hierarchical function spaces and refinement strategies such that wavelet-based approaches excel in terms of grid adaptivity.3. State variables: It was shown that the proposed waveletbased approach naturally gives rise to a mapping algorithm for state variables and does not require intricate patch-recovery techniques that are for instance used in FE-based schemes, cf.[53].4. Scale transition: It was shown that the vanishing moment property of the lifted-interpolating wavelet family considered leads to an efficient integration of the domain and, hence, to an efficient computational homogenisation scheme.Moreover, a closed-form relation for the macroscale algorithmic consistent tangent stiffness tensor could be derived.5. Hill-Mandel condition: It was shown that by accounting for the periodicity of microscale fields in the wavelet syn-thesis and wavelet analysis operations the Hill-Mandel consistency condition could naturally be accounted for.6. Meshfree methods: The proposed wavelet-based approach is intrinsically mesh-free which allows experimental pixel-based data to be efficiently processed as compared to FE-based schemes where the creation of conforming discretisations for complex multiphase microstructures still poses a challenge.However, due to the continuity of the chosen wavelet basis functions, oscillations emerged in the vicinity of weak discontinuities in the solution profile.In this regard, the local adaptivity of the wavelet-based approach is a key property compared to well-established FFT-based spectral approaches.In particular, the oscillations were significantly reduced by activating higher level grid points in the vicinity of phase-boundaries and elasto-plastic transition zones that may evolve in time.This is a remarkable property since, from a physical point of view, these are the locations in the solution domain where a fine mesh is required to accurately represent small-scale features.
The characteristic properties of the proposed waveletbased multiscale formulation were demonstrated in the present contribution by a thorough study of elementary onedimensional problems for which analytical solutions are available.The results are promising since general constitutive relations could be accounted for, since significant reductions in the number of degrees of freedom and the number of material point evaluations were achieved, and since the natural adaptivity of the wavelet-based approach allowed for (evolving) localised microscale features to be accurately resolved.Based on these results, a detailed analysis of the overall efficiency of the proposed formulation when applied to complex evolving microstructures in a multi-dimensional setting is aspired in future contributions.In this regard, it is noted that the proposed one-dimensional approach can naturally be extended to a multi-dimensional setting by adopting, for instance, a tensor-product ansatz.

Fig. 1
Fig.1Hat-wavelets and scaling functions on a periodic domain with period 2.0 mm.The domain is discretised with two grid points on discretisation level-0 and four grid points are used in the corresponding level-1 discretisation.Crosses indicate the spatial positions and the discretisation levels on which grid points occur for the first time.It is

Fig. 2
Fig. 2 Nested (periodic) dyadic grid with period 2 mm.Grid point numbers and associated basis functions are exemplarily provided.Black colour indicates the first occurence of a grid point.

Fig. 7
Fig. 7 Two-scale boundary value problem featuring a polycrystalline microstructure that is subjected to periodic deformations

4 . 1 .Fig. 8
Fig.8 Boundary value problem refinement strategy on a periodic (nested) dyadic grid with period 2.0 mm. a Exemplary wavelet discretisation of a one-dimensional boundary value problem.The domain B with boundary ∂B is indicated in gray-colour, ×-symbols represent domain collocation points, and ⊗-symbols are used for boundary col-

2 .
Grid points x j+1 2k+1 that are associated with wavelet coefficients | [d ω ] j k | > r , j ≥ 0, where r denotes the refinement tolerance, are added to C. 3.For each grid point x j k ∈ C that is associated with a wavelet coefficient | [d ω ] m n | > r , m ≥ 0, neighbouring grid points on the same (i.e.{x j k−2 , x j k+2 }) and on the next higher resolution level (i.e.{x j+1 2k−1 , x j+1 2k+1 }) are added to C, see Fig.

Fig. 10
Fig. 10 Micro-fluctuation field ω and relative error for a prescribed macroscale strain state, ε M = 0.005, and various refinement tolerances r .The relative errors in (d) for r = 10 −9 mm are smaller than 10 −14

Fig. 12
Fig. 12 Dyadic grid and relative error in the micro-fluctuation field for a prescribed macroscale strain state, ε M = 0.005.Lifted interpolating wavelets with N = 2, Ñ = 1 and three-point forward differences

Fig. 14
Fig.14 Micro-fluctuation field ω and relative error for a prescribed macroscale strain state, ε M = 0.005, and various refinement tolerances r

Fig. 16
Fig.16 Dyadic grid and relative error in the micro-fluctuation field for a prescribed macroscale strain state, ε M = 0.005.Lifted interpolating wavelets with N = 2, Ñ = 1 and two-point forward differences

Fig. 17 Fig. 18
Fig.17 Relative errors in macroscale stresses σ M and in algorithmic consistent tangent stiffness E M = dσ M /dε M for a prescribed macroscale strain state, ε M = 0.005

Fig. 19
Fig.19 Evolution of plastic strain ε p at active grid points for lifted interpolating wavelets with N 2, Ñ = 1 and a two-point forward difference scheme.The simulation results are provided for r = 10 −8 mm and a maximum number refinement grids n G = 8 is take into account

( a )Fig. 20
Fig.20 Number of active grid points and number of additional microscale boundary value problems to be solved for lifted interpolating wavelets with N = 2, Ñ = 1, using three-point forward differences and a maximum number of refinement grids n G = 8

Table 2
Filter coefficients of interpolating wavelets with N