To the problem (1), we suppose that u^ε (x, t) can be expressed in the following asymptotic form:

Due to y = x/ε ∈ Y (Y is the unit cell, where Y = (0, 1)ⁿ, n = 2, 3), the chain rule can be denoted as , then substituting Eq. (2) into Eq. (1) and matching terms of the same order of ε , we can immediately obtain

After that, a series of equations are obtained if the coefficients of ε ⁱ (i = − 2, − 1, 0, 1, … ) from both sides of Eq. (3) are set to be equal to each other, which are listed as follows:

Now we solve the asymptotic expansion terms recursively. Firstly, applying the standard asymptotic homogenization theory to Eq. (4), we know that the first term of the asymptotic expansion (2) only depends on macroscopic scale x, which can be expressed as u⁽⁰⁾ = u⁽⁰⁾(x, t).

Secondly, it can be easily recognized that equation (5) is linear with respect to u⁽⁰⁾, so u⁽¹⁾ can be defined as

where N^{α ₁}(x, y, t) is the first-order auxiliary cell function defined in unit cell Y. Substituting u⁽¹⁾ into Eq. (5), we obtain the following equation with attaching boundary condition:

Remark 1 According to Lax– Milgram theorem and assumption (A), it is easy to prove that problem (7) has the unique solution for any fixed time variable t.

Then, inspired by Refs. [7] and [11], we can be define the homogenized material coefficients , , and â _ij(t) as follows:

Further, we can define the homogenized problem as follows:

Remark 2 Under assumptions (A) and (B), we can prove that the homogenized coefficient tensor â _ij(t) is symmetric, positive-definite, and bounded (see Ref. [7]). Besides, it is easy to prove that and are bounded. In addition, the homogenized problem (9) has a unique solution according to the conclusion in Ref. [12].

Finally, from Eqs. (6) and (9), the following equation is obtained:

Thus u⁽²⁾ can be defined as follows:

where N^{α ₁α ₂}(y, t), F(y, t), and G(y, t) are the second-order auxiliary cell functions defined in unit cell Y. Substituting Eq. (11) into Eq. (10), we obtain a series of equations which attach the zero boundary condition

Remark 3 According to the Lax– Milgram theorem and assumptions (A) and (B), it is easy to prove that problems (12)– (14) have a unique solution for any fixed time variable t.

Summing up, we define the FOTS and SOTS of problem (1) as follows:

Notably, the difference from the wave equation in Ref. [7] lies in that there exist second-order auxiliary cell functions F(y, t) and G(y, t) which are used to eliminate the effects of the inertial term and the damping term.

Theorem 1 The hyperbolic– parabolic equations (1) with rapidly oscillating coefficients have an SOTS approximate solution as follows:

where u⁽⁰⁾ is the solution of the homogenized problem (9), N^{α ₁} is the first-order auxiliary cell function defined by Eq. (7), and N^{α ₁α ₂}, F, and G are the second-order auxiliary cell functions defined by Eqs. (12)– (14), respectively.

In this subsection, we give the specific error analysis of FOTS and SOTS in the pointwise sense. Firstly, we denote the error functions of FOTS and SOTS as ω ⁽¹⁾(x, t) = u^{(ε )}(x, t) – u^{(1ε )}(x, t) and ω ⁽²⁾(x, t) = u^{(ε )}(x, t) – u^{(2ε )}(x, t). Suppose that Ω is the union of the entire periodic cells, i.e., , where the index set . Let E_z = ε (z + Y) and ∂ E_z be the boundary of E_z. After a simple computation by substituting ω ⁽¹⁾(x, t) for u^{(ε )}(x, t) in problem (1), we obtain the following residual equations of FOTS for (x, t) ∈ Ω × (0, T) which holds in the distribution sense:

where

Secondly, using the same method, we derive the following residual equation of SOTS which holds similarly in the distribution sense:

where

Now we can give a conclusion about the error analysis in the pointwise sense. From the residual equation (17), one easily finds that the residual of FOTS is O(1) in the pointwise sense due to the term f₀(x, y, t). This is not acceptable in practical engineering calculation. However, it is easy to see from the residual equation (18) that the residual of SOTS is O(ε ) in the pointwise sense. Thus, the SOTS can have the required accuracy for engineering calculations and capture the the micro-scale behavior of composite materials when ε is a sufficiently small constant. This is the main reason and motivation to develop the SOTS.

As we all know, the auxiliary cell functions are defined with the periodic boundary conditions in the classical homogenization methods (see Refs. [22]– [24]). In this case, the auxiliary cell functions have enough regularity on the boundary of the unit cell. However, the normal derivatives of the auxiliary cell functions defined with the Dirichlet boundary conditions only are continuous on the boundary of the unit cell under the geometric symmetry and regularity assumptions of the hyperbolic– parabolic equation’ s coefficients (see Refs. [4]– [10]).

To obtain the regularity of the cell functions, we make the following assumptions similar to Refs. [6]– [10]:

(i) ; ;

(ii) let Δ ₁, … , Δ _n (n = 2, 3) be the middle superplanes of the reference cell Y = (0, 1)ⁿ, then set y = x/ε and assume that a_ii(y, t), ρ (y, t), μ (y, t) are symmetric with respect to Δ ₁, … , Δ _n(n = 2, 3) and a_ij(y, t), i ≠ j are anti-symmetric with respect to Δ ₁, … , Δ _n (n = 2, 3) in Y for any fixed t ∈ [0, T].

Lemma 1 Under assumptions (A)– (C) and (i) and (ii), the normal derivatives σ _Y(N^{α ₁}), σ _Y(N^{α ₁α ₂}), σ _Y(F), and σ _Y(G) can be proved to be continuous on the boundary for any fixed t ∈ [0, T] by using the same method in Refs. [6]– [10], where operator .

Theorem 2 Suppose that Ω ⊂ ℝ ⁿ is the union of the entire periodic cells, i.e., , where the index set . Let u^ε (x, t) and u⁽⁰⁾(x, t) be the weak solutions of problem (1) and the associated homogenized problem (9). The second-order two-scale solution is defined in Theorem 1. Under assumptions (A)– (C) and (i) and (ii), if , u⁽⁰⁾ ∈ L²(0, T; H⁴ (Ω )), ∂ u⁽⁰⁾/∂ t ∈ L²(0, T; H² (Ω )), ∂ ²u⁽⁰⁾/(∂ t²) ∈ L²(0, T; H² (Ω )); then we have the following error estimate of SOTS for the hyperbolic– parabolic equations with rapidly oscillating coefficients:

where C(T) is a constant independent of ε , but dependent on T.

Proof Firstly, from Eqs. (7), (8), and (12)– (14), we have the following result:

This formula will be used in the following proof.

Secondly, we use the residual equation (18) to complete the error estimate of SOTS for the hyperbolic– parabolic equations with rapidly oscillating coefficients. Since ∂ ω ⁽²⁾/∂ t is in L^∞ (0, T; L²(Ω )), we cannot use ∂ ω ⁽²⁾/∂ t as the test function for residual equation (18) directly. To overcome this difficulty, we use the density argument (see Refs. [6], [7], and [24]). In order to simplify the process of proof, we omit this detail. Then, multiplying ∂ ω ⁽²⁾/∂ t on both sides of Eq. (18) and integrating by parts on Ω , we obtain

Using Green’ s formula to simplify the above identity (21), we obtain

where r₂ results from using the Green’ s formula on ∂ E_z (see Refs. [25] and [26]).

Combining Eq. (20) and Lemma 2, we obtain

After that, we can easily derive the following identity by combining Eqs. (22) and (23):

Then, we integrate both sides from 0 to t (0 < t ≤ T) and obtain

After a simple calculation, we obtain the following equation:

Then substituting the initial conditions and the boundary conditions of residual equation (18) into the above identity (26), we have

So far, we have obtained the essential identity for Theorem 3. Starting from here, we will give the final proof. First of all, owing to assumptions (A) and (B) and using the Poincaré – Friedrichs inequality, we have the following inequality by transforming the left side of Eq. (27):

Then, due to and Schwarz’ s inequality,

After that, using the Young’ s inequality ab ≤ (λ a² + λ b²/λ )/2, ∀ λ ∈ ℝ ⁺ with parameter λ and assumptions (A)– (C), and choosing the same parameter λ at ∀ τ ∈ [0, T], we obtain the following inequality by transforming the right side of Eq. (27):

Choosing a sufficiently big λ > 0 fulfilled λ − 2μ ⁰ > 0 and combining Eqs. (28) and (30), we have

Denoting δ ₁ = min(ρ ⁰, C₁) and δ ₂ = max(M, λ − 2μ ⁰), we have

Without loss of generality, we denote C(T) = C(T)/δ ₁ and δ = δ ₂/δ ₁ and obtain the following inequality:

Setting

then we obtain from Eq. (33). It follows from Gronwall’ s inequality (see Ref. [24]) that Θ (t) ≤ C(T)ε ². consequently, we obtain the following result:

Then by using the AM– GM inequality to the left side of inequality (34), the following equality is obtained:

With the arbitrariness of time variable t and squaring root on both sides of inequality (35), we obtain the final convergence result

where C(T) is a constant independent of ε , but dependent on T.

In this section, we give the SOTS’ s convergence theorem of the hyperbolic equation by the degradation of Theorem 2. In general, for ρ ^ε ≠ 0 and μ ^ε = 0, equation (1) is called the hyperbolic wave equation, which describes the wave scattering and diffraction behavior of the composite structure. In the past, the convergence analysis of the hyperbolic equation is done separately. Now, we can obtain the convergence results of the hyperbolic equation from the degradation of Theorem 2.

Theorem 3 Assume that the definitions of Ω , , T_ε , u^ε (x, t), u⁽⁰⁾(x, t) are the same as those in Theorem 2. It is easy to know that and G(y, t) = 0 when ρ ^ε ≠ 0 and μ ^ε = 0. Under assumptions (A)– (C) and (i) and (ii), if , u⁽⁰⁾ ∈ L²(0, T; H⁴ (Ω )), ∂ u⁽⁰⁾/∂ t ∈ L²(0, T; H² (Ω )), ∂ ²u⁽⁰⁾/∂ t² ∈ L²(0, T; H² (Ω )); then we have the following error estimate of SOTS for the hyperbolic wave equation with rapidly oscillating coefficients:

where C(T) is a constant independent of ε , but dependent on T.

Proof Under the situation of ρ ^ε ≠ 0 and μ ^ε = 0, the term in Eq. (22) is equal to 0. From the proof of Theorem 2, it is easy to know that the error estimate of SOTS for the hyperbolic wave equation still fulfils inequality (31)

where the constant μ ⁰ = 0 due to μ ^ε = 0.

Then using the same proof method and denoting δ ₁ = min(ρ ⁰, C₁) and δ ₂ = λ , we can obtain the following inequality:

After that, using Gronwall’ s inequality (see Ref. [6] and [7]) and AM– GM inequality to inequality (42), and mimicking the proof process of Theorem 2, we obtain the final convergence result of the hyperbolic wave equation as follows:

The obtained convergence result of the hyperbolic wave equation is in completely accord with that in Ref. [7].

It is easy to see that auxiliary cell problems (7) and (12)– (14) are all elliptic equations, where the time variable t only plays the role of a parameter. So we can solve these auxiliary cell problems with the standard finite element method. Considering the internal geometric structure complexity of the unit cell, we should use the triangular element in ℝ ² and the tetrahedral element in ℝ ³ to solve the auxiliary cell problems. Without loss of generality, we only discuss the two-dimensional problem in this paper. Let J^h = {K} be a regular family of triangulation of the reference unit cell Y, where h = max_K{h_K}. Define the linear conforming finite element space for solving the auxiliary cell problems as .

As we all know, the homogenized coefficients â _ij(t) are derived from Eq. (8). So â _ij(t) depend on the finite element calculation of the auxiliary cell problem N^{α ₁}(y, t). Suppose that are the corresponding finite element solutions of N^{α ₁}(x, y, t). If we denote as the finite element approximation of the homogenized coefficients â _ij(t), then will be calculated by the formula

Therefore, we need to solve the following modified homogenized equation:

According to our knowledge, the numerical accuracy of the modified homogenization problem (41) directly affects the accuracy of SOTS. To improve the accuracy of the modified homogenization problem (41), we solve the modified homogenized equation with the quadratic element. Let J^h₀ = {e} be a regular family of triangulation of the spatial region Ω , where h₀ = max{h_e}. Define the quadratic finite element space as . Now, we can obtain the semi-discrete scheme for solving problem (51). Firstly, we give the discrete variational formulation for modified homogenization problem (41) as follows:

Secondly, let be a basis of V_h0 and denote . Then, the discrete variational inequality (42) is equivalent to the following linear ordinary differential system when substituting into inequality (42):

where ζ (t) = [ξ ₁(t), ξ ₂(t), … , ξ _M(t)]^T, M(t), C(t), and K(t) are the mass matrix, the damping matrix, and the stiffness matrix, respectively, and Q(t) is the load vector (see Ref. [21]), which can be obtained by the following expressions:

Now, we introduce the Newmark scheme for solving the ordinary differential system (44). As we all know, the Newmark difference scheme is an unconditional stable and convergent difference scheme when scheme parameters δ ≥ 0.5 and α ≥ 0.25 (δ + 0.5)² (see Refs. [21] and [25]) and we choose δ = 0.5 and α = 0.25 in this paper. On the other hand, from Theorem 1, we know that the calculation of SOTS needs the values of acceleration ∂ ²u⁽⁰⁾/∂ t² and velocity ∂ u⁽⁰⁾/∂ t at every time step. The Newmark scheme not only can solve u⁽⁰⁾, but also simultaneously solve ∂ ²u⁽⁰⁾/∂ t² and ∂ u⁽⁰⁾/∂ t. So we use the Newmark scheme to solve the ordinary differential system (43) at any time step. The detailed procedure of the Newmark scheme for system (43) is listed as follows.

Step 1 Construct mass matrix M(t), damping matrix C(t), stiffness matrix K(t), and load vector Q(t).

Step 2 Use the initial and boundary conditions ζ (0) = ζ ⁰ and and linear ordinary differential system to compute the initial .

Step 3 Choose the time step Δ t and parameters δ , α . For each time step, we should know ζ (t), , at first. Then we solve the unknown ζ (t + Δ t) by using the following formula (see Ref. [13]):

where

Finally, we can derive , by the following formulas:

Step 4 Repeat step 3 for all time steps, then we can obtain the value of the ordinary differential system (43) at any time step.

In this subsection, we present a second-order two-scale numerical method to efficiently solve the hyperbolic– parabolic equations with rapidly oscillating coefficients. Firstly, we give the algorithm of the second-order two-scale numerical method. After that, we give the analysis of SOTS’ s computational efficiency and superiority to explain why the second-order two-scale numerical method can greatly reduce the calculation compared with the direct finite element method.

Step I Generate a geometric model of the unit cell and compute auxiliary cell functions , , F_h(y, t), G_h(y, t) on the reference cell Y = (0, 1)ⁿ (n = 2, 3) by solving the auxiliary cell problems (7) and (12)– (14), respectively.

Step II Calculate the homogenized coefficients , , by taking the integral of Eq. (8). Then solve numerically the modified homogenized hyperbolic– parabolic equation (41) on the whole domain Ω × [0, T] with a coarse mesh and a large time step by using the method presented in subsection 4.1.

Step III Solve the spatial variable derivatives ∂ ũ ⁽⁰⁾/∂ x_{α 1}, ∂ ²ũ ⁽⁰⁾/∂ x_{α 1}∂ x_{α 2} of the homogenized solutions ũ ⁽⁰⁾(x, t) by the average technique on relative elements (see Refs. [2] and [6]) and the temporal variable derivatives ∂ ũ ⁽⁰⁾/∂ t, ∂ ²ũ ⁽⁰⁾/∂ t² by using the difference schemes (46) in subsection 4.1. The following formulas are used for evaluating the high order partial derivatives of the spatial and temporal variables:

where M is the node of FEM set J^h₀, J^h₀(M) denotes the set of elements with the same node M, τ (M) is the number of elements in J^h₀(M), denotes the derivative at node M on element e and at time t = t_i;

where m is the number of node freedom on element e, and is the Lagrange’ s type shape function corresponding to node p_j on e;

where ζ _M(t_i) denotes the value of node M at time t = t_i, , are the first-order and the second-order partial derivatives of the temporal variables at time t = t_i.

Step IV Use the auxiliary cell functions in step 1 to assemble SOTS as follows:

Step V For arbitrary point (x, t) ∈ Ω × [0, T], we use the higher-order interpolation method to obtain the high-precision SOTS value as follows:

where

, , F_2h, G_2h likeness,

, F_2h(t), G_2h(t) likeness, and the higher-order interpolation operator and can be found in Refs. [2]– [7] and [26].

The main motivation to develop multiscale methods is to solve the problems which cannot be solved by the conventional numerical methods. So the analysis of the computational efficiency of multiscale methods only can be done when we assume that the multiscale problems can be solved by the conventional numerical methods. Firstly, we give a specific analysis to explain why the SOTS can greatly reduce computational cost compared with the direct FEM. We assume that the whole domain Ω consists of n × n entire periodic cells. If we need m elements to mesh every unit cell in Ω , the total number of elements is n × n × m in the direct FEM calculation. Fortunately, SOTS only needs the grid information of a single cell to solve the auxiliary cell problems. From Section 2, we know that there are two auxiliary cell problems N^{α ₁}(y, t), four N^{α ₁α ₂}(y, t), one F(y, t), and one G(y, t) for the two-dimensional hyperbolic– parabolic equations. Thus, the total number of elements to solve the auxiliary cell problems is 8m. Luckily, this number does not increase even if the cell number of the whole domain Ω increases. On the other hand, we assume that the elements for solving the homogenization problem is c × n × n × m (0 < c ≪ 1). In the general calculation, c is in the range from 1/20 to 1/5 (see Refs. [3]– [7], and [11]). In summary, the total element number of SOTS is (cn² + 8m) at every time step. With the rapid increase of cell number n × n, it is easy to know SOTS can save 19/20– 4/5 of the elements in the two-dimensional hyperbolic– parabolic equations. By a similar analysis, we know that SOTS for the three-dimensional equations still can save 19/20– 4/5 of the elements. So one can understand why SOTS can greatly reduce the computing resources and improve the computational efficiency.

Now, we give a detailed analysis of the SOTS’ s stability. Due to the great difference of material coefficients between different kinds of materials in the composite structure, if we use the finite element method to solve problem (1) directly, the overall damping matrix and stiffness matrix are ill-conditioned matrices. Although they are all symmetric definite matrices, the difference between the maximum eigenvalue and the minimum eigenvalue of the overall damping matrix and the stiffness matrix is too big. So the numerical results will be unstable, which causes a larger error. Fortunately, if we use SOTS to solve problem (1), the overall damping matrix and the stiffness matrix have good properties because the homogenized material coefficients of homogenization problem (41) are functions which vary slowly. Therefore, the numerical algorithms we presented have good stability compared with the direct FEM. At the same time, we can use a large time step to solve homogenization problem (41) because the overall damping matrix and the stiffness matrix are both well-conditioned matrices.