Chebyshev Spectral Element Analysis for Pore-Pressure of ERDs during Construction Period

Show more

1. Introduction

The finite-element method (FEM) has been the main tool of choice for computational analysis of ERD’s pore-pressure for a long time [2] . FEM remains under continuing development to improve its accuracy, efficiency and reliability with effect. The spectral element method (SEM) combines the high accuracy of spectral analysis with the adoptability of FEM for the complex geometrical configuration of problems. The spectral analysis originates from the Ritz and Galerkin principles to approximate the solutions of differential equations by using orthogonal polynomial of partial differential equation with Sturm-Liouville differential operators as basis expansion of the unknowns. The spectral polynomial poses excellent convergence and accuracy with $O\left({N}^{-1}\right)$ convergence order, where N is the term number of orthogonal series. Despite of excellent features of spectral analysis, it is difficult to handle with the problems in complicated geometrical domains. To overcome this kind of disadvantages of spectral analysis methods, Orszag develops the means of mapping complicated domain into rectangular domain by the coordinate transformation. Even though so, Orszag’s mapping method is yet limited to such simpler problems as the spatial domains decomposed into several polygonal segments. In 1984, Patera (1984) adopted spectral isoparametric transformation, which is similar to the ways in FEM, by spectral polynomials to map a typical irregular element into a standard square spectral-element (SE). This numerical philosophy has been named as spectral element method (SEM) now [1] - [8] .

There have been 86,000 barrages now in China and the ERD account for 95% or more (Lin 2003, Gu 2006) [4] [6] . The high construction rate of ERD is mainly for richly site damming materials such as soil, sand, gravel, rock and other land wastes, and relatively lower costs of foundation treatment than concrete dams. The most of ERD are built by vibro-roller compaction in layered filling embankment-dam. The damming materials, after dam bodies being poured out, are compacted by vibrational compaction equipments to increase the dry density and shear strength of filling soil. The whole settlement of dam is fixed by the pore-pressure and total stress of ERD during construction period. The pore-pressure in the filling material, especially for those with large compressibility and low permeability, of dam body changes greatly for roller compaction under construction.

Biot’s consolidation theory (Ferronato 2010) [2] is an appropriate way to formulate the problem of pore-pressure in constructed dam with the change of total stress as a result of compaction. As well known, the numerical solution to the Biot partial differential equation is a quite difficult affair, for the reason of fully coupled Biot’s model with complicated algebraic systems and nonlinear impediments such as convergence and stabilities of numerical solution of spatial discretization by FEM. Nonlinearity is not a major complication for spectral methods per se. The Chebyshev series without singularities anywhere in the complex plane has super-geometric convergence, i.e. as for the expansion factor is that ${a}_{n}~O\left(\text{exp}\left\{-\left(n/j\right)\text{log}\left(n\right)\right\}\right)$ . It is preferable to use spectral method to integrate Galerkin projection equations in time history in this paper, and it has been verified as unconditionally stable in practice by numerical examples.

In the present paper a SEM in spatial domain, which adopts Chebyshev polynomials as the basis functions of asymptotic expansion of the unknowns, is developed to solve numerically the pore-pressure of ERD compacted by vibration compactor during construction by the governing partial differential equation of Biot consolidation [3] . With regard to non-periodic pore-pressure problems of ERD the Chebyshev spectral methods are common way of asymptotical formulation for the unknowns. The spectral isoparametric transformation is used to handle with the polygonal configuration of ERD, each of spectral elements is mapped into a regular square with 2 unit length of sides just same as in FEM. Unlike ordinary numerical integral methods of spatial discretization equations to time history, we employ Fourier spectral expansions to integrate temporal matrix differential equations so avoiding the predicament of numerical stabilities. A few practical engineering cases of ERD are conducted by the SEM proposed in this paper, through that indicates some interesting and meaningful results of numerical procedure of SEM.

As per numerical examples of compacted ERD one could find out that only one spectral element should be taken for the domains with homogenous materials to get same and even better accuracy than ones of many finite elements, but the spectral polynomials must expand to certain numbers of terms. So, for the problem of ERD with zonal uniform material properties the computational efficiency of SEM seems to be more advantages than one of finite element [5] . In the procedure of domain discretization by Chebyshev spectral expansion, the spectral and finite element have not different in essence, but the biggest difference is that the spectral element can have arbitrary interpolation points only in one spectral element and finite element does not have this feature. Because of this, the SEMs are with the greatest limitation of application in general engineering if have not suitable pre-processor. The numerical calculation indicates that SEM does not need to pay the extra price than FEM, except for different formulation of interpolation functions.

2. General Formulation of Biot’s Consolidation

2.1. Formulation of ERD’s Pore-Pressure during Construction Period

We put forward the model of ERD’s pore-pressure during construction period, which is similar to the classical Biot’s consolidation with a different meaning of the compressibility for a non-completely saturated and compressible solid matrix. The general fluid mass balance can be prescribed by continuity equation (Gaspar, Lisbona 2003, Boyd 2000) [1] [3] :

$\nabla \cdot v+\frac{\partial}{\partial t}\left(\varnothing \beta p+\delta \nabla \cdot d\right)=f$ (1)

where ∇ is Hamilton operator, t time variable, ∅ the medium porosity, β the fluid compressibility, p the fluid pore-pressure, δ the Biot coefficient, d the medium displacement and v is Darcy’s velocity, f a flow source or sink. Equation (1) must be coupled with Darcy’s law (Boyd 2000) [1] :

$v+\frac{k}{\rho g}\nabla p=0$ (2)

where k is the hydraulic conductivity and ρg is the fluid specific weight.

For the pore-pressure evaluation of ERD under construction, we must modify general pore-pressure Equation (1) and Equation (2) corresponding to the different periods and conditions of ERD construction. During the ERD rolling compaction by smooth drum vibrating roller, the porous fluid can be thought as non-compressible and no flow source or sink exists, i.e. β and f equal to zero. So that we can arrive at the combination of Equation (1) and Equation (2), along with the change of total tensor within the bodies of ERD, suitable to rolling compaction ERD as follow (Gu 2006) [4] :

$\frac{\partial p}{\partial t}=\frac{1+\varphi}{{a}_{v}\rho g}\nabla \cdot \left(k\nabla p\right)+\stackrel{\xaf}{B}\frac{\partial \sigma}{\partial t}$ (3)

where ${a}_{v}$ is the core compressibility, σ total tensor, $\stackrel{\xaf}{B}$ coefficient of initial core pressure.

The evaluation of pore-pressure of ERD depends on the nature of damming material and construction rate of ERD (Michael 2006) [7] . If the hydraulic conductivity $k\le {10}^{-7}\text{cm}/\text{s}$ for high ERD or $k\le {10}^{-6}\text{cm}/\text{s}$ for medium ERD, the pore-pressure keeps constant value, i.e. initial pore-pressure ${u}_{0}=\stackrel{\xaf}{B}\sigma $ because of it is difficult for pore-pressure to dissipate in the ERD. The structure sketch of roller compaction ERD is shown as Figure 1.

2.2. One-off Rapid Embankment of ERDs

When the filling ERDs with low height of dam and rapid embankment, we often think the embankment of ERDs as one-off so the governing equation of pore-pressure can be reduced into one-dimension formulation.

$\frac{\partial p}{\partial t}=\frac{1+\varphi}{{a}_{v}\rho g}\frac{{\partial}^{2}p}{\partial {x}^{2}}$ (4)

where coordinate x is in the direction of dam surface. The configuration and coordinate of ERDs can be denoted in Figure 2.

In above the consolidation equation becomes a quite straightforward one-dimensional formulation so we can solve it by the approximate analytic method under non-consideration of the change of medium porosity and other parameters. For most general situations it is necessary to arrive at its solutions by numerical ways.

2.3. Layered Filling of ERDs

For high ERDs the filling procedure of ERDs is in layers, we can consider the embankment of ERDs body as evenly filling at the speed of α. Meanwhile the pore-pressure dissipates in filled soil layer of ERDs body and the pore-pressure of newly filling soil layers increases with total tensor adds with time going. The

Figure 1. Sketch of ERDs of layered roller compaction.

Figure 2. Configuration of ERDs at one-off embankment.

consideration above means that we just take into account the pore-pressure of newly filling soil layers during the embankment of ERDs and ignore the porosity change of soil column caused by the dead weight of overburden. This assumption for low and medium ERDs under construction is suitable and even for high ERDs has acceptable accuracy in damming engineering. So the geometrical sketch of one filling layer of ERDs is depicted in Figure 2.

So the increment of total tensor σ of filling soil layers can be expressed as

$\text{d}\sigma =\alpha \text{d}t$ (5)

in time span dt. By a simple formula transformation $p=\stackrel{\xaf}{p}-\frac{1}{2}\stackrel{\xaf}{B}\alpha {x}^{\text{T}}x$ , rolling

compaction ERD’s Equation (3) can be converted into a formulation similar to Equation (4) as following

$\frac{\partial \stackrel{\xaf}{p}}{\partial t}=\frac{1+\varphi}{{a}_{v}\rho g}\nabla \cdot \left(k\nabla \stackrel{\xaf}{p}\right)$ (6)

We notice that governing Equation (6) of layered filling of ERDs have the same formulation with the governing Equation (4) of one-off rapid embankment at one-dimensional problems, so there aren’t special differences in both cases of the calculation of spectral element methods.

2.4. Two-Dimensional Dewatering Filling of ERDs

In the procedure of rapidly filling of high ERDs both bulk of soil skeleton and effective stresses exist the process of the evolution from initial state to current state. For compressible damming soils, aside from the pore-pressure resulted in by unstable seepage, the compression and total tensor change, which are caused by the rolling process of ERDs, of damming soils generate additional pore-pressure. The pore water is ruled out from the ERDs section. Under the present circumstances, the governing pore-pressure equation of two-dimensional dewatering filling of ERDs takes the general Equation (3), in which first part indicates the effect of natural seepage and second part contributes the compression of pore volume to the pore-pressure. The faster the filling of high ERDs proceeds, the faster the pore-pressure grows up, and the better drainage condition of surface and dam foundation of ERDs are, the faster pore-pressure of ERDs dissipates, vice versa.

3. Chebyshev Spectral Element for Pore-Pressure of Filling ERDs

Variational Formulation and Interpolation Strategy

For generality, we use the partial differential Equation (3) as governing equation to solve numerically the pore-pressure of rolling compaction ERDs whether rapid or layered filling procedure. The variational formulation of spectral element method for pore-pressure of filling ERDs is highly similar to the one of finite element model (Pozrikidis 2005) [9] . However, the major difference is that the spectral element introduces the strategy of interpolation which approximates the unknown by a sequence of high-order and non-overlapping polynomials on the element grid nodes. So the unknown pore-pressure $p\left(x,t\right)$ can be approximated by a set of basic functions, $\phi \left(x\right)={\left[{\phi}_{1},{\phi}_{2},\cdots ,{\phi}_{N}\right]}^{\text{T}}$ , here $x={\left[{x}_{1},{x}_{2},{x}_{3}\right]}^{\text{T}}$ and T denote position coordinates and transposed matrix respectively, and N is the total number of spectral element nodes, Let $u\left(t\right)={\left[{u}_{1},{u}_{2},\cdots ,{u}_{N}\right]}^{\text{T}}$ denotes the pore-pressure of N spectral element nodes at time t, so that the pore-pressure in a typical spectral element at any position x and any time t is

$p\left(x,t\right)=\phi {\left(x\right)}^{T}u\left(t\right)$ (7)

where the basic functions φ(x) relates to Chebyshev polynomials:

${T}_{k}\left(x\right)=\mathrm{cos}\left(k{\mathrm{cos}}^{-1}x\right),\text{\hspace{0.17em}}\text{\hspace{0.17em}}k=0,1,\cdots ,M$ (8)

In above $x\in \left[-1,1\right]$ and M is interpolation number in the direction of x-coordinate. For non-periodic evaluation of pore-pressure of filling ERDs, Chebyshev polynomials are one of the best choices because of their exponential Chebyshev convergence for non-periodic functions. In general cases Chebyshev polynomials are defined in a regular three dimensional domain, ${R}^{3}\in \left[-1,1\right]\times \left[-1,1\right]\times \left[-1,1\right]$ , irregular element domains in practical problems must be mapped into Chebyshev regular domain by the following spectral-isoparametric mapping function:

$x=H{\left(\xi \right)}^{\text{T}}XH\left(\eta \right)$ (9)

where

$H\left(\xi \right)=\left[\begin{array}{ccc}h\left(\xi \right)& 0& 0\\ 0& h\left(\xi \right)& 0\\ 0& 0& h\left(\xi \right)\end{array}\right]\text{,}X=\left[\begin{array}{ccc}{X}_{1}& 0& 0\\ 0& {X}_{2}& 0\\ 0& 0& {X}_{3}\end{array}\right]$ (10)

and the submatrices are

$h\left(\xi \right)=\left[{h}_{1}\left(\xi \right),{h}_{2}\left(\xi \right),\cdots ,{h}_{Nx}\left(\xi \right)\right]\text{,}{X}_{i}=\left[\begin{array}{ccc}{X}_{i}\left({x}_{1,1}\right)& \cdots & {X}_{i}\left({x}_{1,Ny}\right)\\ \vdots & \ddots & \vdots \\ {X}_{i}\left({x}_{Nx,1}\right)& \cdots & {X}_{i}\left({x}_{Nx,Ny}\right)\end{array}\right]$ (11)

Nx of Equation (11) denotes the number of nodes in a spectral element and x_{i,j} of Equation (11) is the coordinates x of ith and jth nodes in spectral-isoparametric coordinate system in the direction ξ and η respectively. And then we discussed the implementation of the general elements with boundary curvature mapping into regular spectral elements above [10] .

$\begin{array}{c}{h}_{i}\left(\xi \right)=\frac{2}{N{C}_{i}}T\left(\xi \right)T\left(\xi \right)\\ =\left[\begin{array}{ccc}{T}_{1}\left({\xi}_{i}\right){T}_{1}\left(\xi \right)/{C}_{1}& {T}_{2}\left({\xi}_{i}\right){T}_{2}\left(\xi \right)/{C}_{2}& \begin{array}{cc}\cdots & {T}_{N}\left({\xi}_{i}\right){T}_{N}\left(\xi \right)/{C}_{N}\end{array}\end{array}\right]\end{array}$

${C}_{i}=\{\begin{array}{l}2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}i=0,N\\ 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}i\ne 0,N\end{array}$ (12)

We mapped an element by isoparametric method similar to in FE from the physical xy plane to the familiar rectangular element with 2 unit length of edge in ξη plane in two dimensions. The spectral element interpolation vectorial functions, $\phi \left(x\left(\xi ,\eta \right)\right)$ , can be derived as follow:

$\phi \left(x\left(\xi ,\eta \right)\right)=H{\left(\xi \right)}^{\text{T}}H\left(\eta \right)$ (13)

It is readily verified that the interpolation vectorial functions,
$\phi \left(x\left(\xi ,\eta \right)\right)$ , satisfy the familiar cardinal interpolation properties
${\phi}_{i}\left({x}_{j}\left(\xi ,\eta \right)\right)={\delta}_{ij}$ , here δ_{ij} is Kronecker’s delta [11] .

For obtaining the approximation of Equation (3) in evaluation of general pore-pressure of rolling compaction ERDs, the Galerkin projections using as weighting functions the global interpolation functions, $\phi \left(x\left(\xi ,\eta \right)\right)$ , are available in such integral form:

$\underset{e}{{\displaystyle \sum}}{\displaystyle \underset{{\Omega}^{e}}{\int}\frac{\partial p}{\partial t}\phi \text{d}\Omega}=\underset{e}{{\displaystyle \sum}}{\displaystyle \underset{{\Omega}^{e}}{\int}\frac{1+\varphi}{{a}_{v}\rho g}\nabla \cdot \left(k\cdot \nabla p\right)\phi \text{d}\Omega}+\underset{e}{{\displaystyle \sum}}{\displaystyle \underset{{\Omega}^{e}}{\int}\stackrel{\xaf}{B}\frac{\partial \sigma}{\partial t}\phi \text{d}\Omega}$ (14)

In which, summation symbol Σ means to sum for whole discrete spectral elements and Ω^{e} denotes spectral-elemental domain. Introducing Chebyshev expansion Equation (7) into Equation (14) and applying Gauss integration rule for above equations, we can obtain following matrix equations considering natural boundary conditions,

$M\frac{\text{d}u}{\text{d}t}=Q-Du+\stackrel{\xaf}{B}\frac{\text{d}\sigma}{\text{d}t}$ (15)

where the global matrices M, Q, and D are formulated such as:

${\left(M\right)}_{ij}={\displaystyle \underset{\Omega}{\int}{\phi}_{i}{\phi}_{j}\text{d}\Omega {\left(D\right)}_{ij}}={\displaystyle \underset{\Omega}{\int}\nabla {\phi}_{i}\cdot \nabla {\phi}_{j}\text{d}\Omega {\left(Q\right)}_{i}}={\displaystyle \underset{\Gamma}{\int}\stackrel{\xaf}{C}{\phi}_{i}n\cdot \nabla u\text{d}l}$ (16)

where Ω is the whole solution domain and Γ the boundary contour, and n the outer normal at boundary flux contour. Noticing that the all of integrations in Equation (16) are in natural coordinates xy plane, we must transfer them into mapping coordinates ξη plane. So the Jacobian matrix of mapping from the physical xy plane to the mapping ξη plane was used in surface integral in two dimensions as

$\text{d}\Omega =\text{Det}\left[J\right]\text{d}\xi \text{d}\eta $ (17)

and line integral

$\text{d}l=\{\begin{array}{l}\sqrt{{\left(\frac{\partial x}{\partial \eta}\right)}^{2}+{\left(\frac{\partial y}{\partial \eta}\right)}^{2}}\cdot \text{d}\eta \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\xi =\text{constant}\\ \sqrt{{\left(\frac{\partial x}{\partial \xi}\right)}^{2}+{\left(\frac{\partial y}{\partial \xi}\right)}^{2}}\cdot \text{d}\xi \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\eta =\text{constant}\end{array}$ (18)

where J is Jacobian matrix defined as

$J={\nabla}_{\xi \eta}\times x,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\nabla}_{\xi \eta}=\left[\begin{array}{cc}\frac{\partial}{\partial \xi}& \frac{\partial}{\partial \eta}\end{array}\right]$ (19)

The integral of Equation (16) over the area of the general quadratic element in the xy plane in two dimensions can be depicted as an integral over the area of the standard rectangle in the mapping ξη plane, and the line integral by the same way [12] .

4. Numerical Evaluation

4.1. The Numerical Evaluation of the Layered Embankment of ERDs

As described in Section 2.2, the layered embankment of ERDs can be considered as one-dimensional consolidation problem, the numerical model illustrated in Figure 1. The evolution of pore-pressure in ERDs is governed by the Equation (6) and following boundary condition at initial time for saturated soil.

$\{\begin{array}{l}p=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{when}\text{\hspace{0.17em}}x=0\\ p=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{when}\text{\hspace{0.17em}}x=2l\\ p=\sigma ,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{when}\text{\hspace{0.17em}}0<x<2l\end{array}$ (20)

where σ means total stress in ERDs. The physical domain [0,2l] is mapped

into a regular domain [−1,1] in spectral element method by $\zeta =-1+\frac{x}{l}$ . So the

governing Equation (6) is rewritten as this:

$\frac{\partial p}{\partial t}=\frac{1+\varphi}{{a}_{v}\rho g{l}^{2}}\frac{{\partial}^{2}p}{\partial {\zeta}^{2}}$ (21)

The pore-pressure function, p(ζ,t), defined over the ith element can be approximated with the jth-degree Chebyshev polynomial, T_{j}(x), at interpolation

points ${x}_{i}=\mathrm{cos}\frac{i\text{\pi}}{N}$ here N the total numbers of interpolation points along the x-direction, as

${p}_{N}\left(\zeta ,t\right)=\underset{j=0}{\overset{N}{{\displaystyle \sum}}}\text{\hspace{0.05em}}{T}_{j}\left(\zeta \right)\stackrel{\xaf}{p}\left({\zeta}_{j},t\right)$ (22)

In there $\stackrel{\xaf}{p}\left({\zeta}_{j},t\right)$ are a priori unknown nodal values at the Chebyshev interpolation points ${\zeta}_{j}$ in the one-dimensional pore-pressure problem. The integral weak formula corresponding to Equation (21) depicts as

$\underset{-1}{\overset{1}{{\displaystyle \int}}}\left(\frac{\partial {p}_{N}}{\partial t}-\frac{1+\varphi}{{a}_{v}\rho g{l}^{2}}\frac{{\partial}^{2}{p}_{N}}{\partial {\zeta}^{2}}\right)\frac{{T}_{j}\left(\zeta \right)}{\sqrt{1-{\zeta}^{2}}}\text{d}\zeta =0,\text{\hspace{0.17em}}j=1,2,\cdots ,N-1$ (23)

As per the orthogonality of Chebyshev polynomials with the weight function,

$\frac{1}{\sqrt{1-{\zeta}^{2}}}$ , in the domain [−1,1], we can readily arrive at one-order ordinary

differential equation on time t and the recurrent formula:

$\begin{array}{l}\frac{\partial \stackrel{\xaf}{p}\left({\zeta}_{j},t\right)}{\partial t}-{b}_{j}^{\left(2\right)}\frac{1+\varphi}{{a}_{v}\rho g{l}^{2}}\stackrel{\xaf}{p}\left({\zeta}_{j},t\right)=0\hfill \\ {A}_{j-2}{b}_{j-2}^{\left(2\right)}+{B}_{j}{b}_{j}^{\left(2\right)}+{C}_{j+2}{b}_{j+2}^{\left(2\right)}=2j\stackrel{\xaf}{p}\left({\zeta}_{j},t\right)\hfill \end{array}\text{\hspace{0.17em}}j=2,\cdots ,N-1$ (24)

where

$\begin{array}{l}{A}_{j-2}=\frac{{\stackrel{\xaf}{c}}_{j-1}{\stackrel{\xaf}{c}}_{j-2}}{2\left(j-1\right)}{B}_{j}=-\left(\frac{{\stackrel{\xaf}{c}}_{j-1}}{2\left(j-1\right)}+\frac{{\stackrel{\xaf}{c}}_{j}}{2\left(j+1\right)}\right){C}_{j}=\frac{1}{2\left(j+1\right)}\hfill \\ {b}_{N+1}^{\left(2\right)}={b}_{N}^{\left(2\right)}=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em} c \xaf j = { 2 \hspace{0.17em} \hspace{0.17em} \hspace{0.17em} \hspace{0.17em} \hspace{0.17em} \hspace{0.17em} \hspace{0.17em} j = 0 1 \hspace{0.17em} \hspace{0.17em} \hspace{0.17em} \hspace{0.17em} 1 \u2264 j \u2264 N \u2212 1}\hfill \end{array}j=2,\cdots ,N-1$ (25)

In the situation of layered embankment of ERDs, the pore-pressure of ERDs is solved by the one-order ODEs (24) in time t and the factors of ODEs (24) are achieved by the recurrent formula (24). In the derivation of formula (24), we use well the orthogonality of Chebyshev polynomial and the periodicity of basic functions. It is a common way in modern numerical method that an unknown is approximated by expansion of a system of orthogonal functions. Theoretically, to a fully smooth a periodic unknown, the solution can easily achieve a spectrum accuracy as long as appropriate expansion functions are chosen. Chebyshev polynomials are a kind of eigen-functions of Sturm-Liouville problem, for any smooth unknown when Chebyshev polynomials are employed to approximate an unknown the expansion function can achieve spectrum accuracy and without limit of boundary condition [13] .

The Equation (24) and (25) must be combined with boundary conditions at terminal points in the domain [−1,1] to solve the unknown functions $\stackrel{\xaf}{p}\left({\zeta}_{j},t\right)$ . Whether spectral Equation (15) or (24), they are all the one-order ordinary differential equation at time t, so we must integrate them with calculus of differences or some other integral methods. For the reason of high precision of spectral method, it is best choice for integrating one-order ODEs aforementioned with the spectral method without any difficulties.

4.2. Time-History Integration Based on Spectral Method

Without loss of generality, we can put Equation (15) and (24) into following general formulation, i.e. a standard problem of evolution of this form (Michael S. 2006) [7] :

$\{\begin{array}{l}\frac{\text{d}}{\text{d}t}{\rm E}\left(t\right)=f\left({\rm E}\left(t\right),t\right)\\ {\rm E}\left(0\right)={{\rm E}}_{I}\end{array}$ (26)

where Ε_{I} is the initial condition of evolution problem and the function Ε(t) depicts a vector with n components, consisting of a group of unknown functions, viz
${\rm E}\left(t\right)=\left[{E}_{1}\left(t\right),{E}_{2}\left(t\right),\cdots ,{E}_{n}\left(t\right)\right]$ , and
$f\left({\rm E}\left(t\right),t\right)$ is a known non-linear vectorial function of unknown Ε(t) and time t. The time variable t is defined in the time-span [0, T]. The time-span [0, T] could be divided into N subdivisions,

so we use Fourier polynomials, ${F}_{k}={C}_{K}{\text{e}}^{-iK\frac{2\text{\pi}t}{T}}\left(K=1,2,\cdots ,N\right)$ where $i=\sqrt{-1}$

is an imaginary unit and C_{K} undetermined coefficients, to expend the temporal function Ε(t)

${\rm E}\left(t\right)=\stackrel{^}{A}F$ (27)

where $\stackrel{^}{A}$ is a second rank tensor and F a vector

$\stackrel{^}{A}=\left[\begin{array}{cccc}{A}_{11}& {A}_{12}& \cdots & {A}_{1N}\\ \vdots & \vdots & \ddots & \vdots \\ {A}_{n1}& {A}_{n2}& \cdots & {A}_{nN}\end{array}\right]$ (28)

$F={\left[{F}_{1},{F}_{2},\cdots ,{F}_{N}\right]}^{\text{T}}$ (29)

and

$\frac{\text{d}{F}_{K}}{\text{d}t}=-iK\frac{2\text{\pi}}{T}{F}_{K}$ (30)

Substituting Equation (4.8)~(4.11) into (4.7) and application of orthogonality of Fourier polynomial, we can get

${A}_{mK}=\frac{i}{K{\delta}_{K}}\underset{0}{\overset{T}{{\displaystyle \int}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\stackrel{\xaf}{{F}_{K}}\cdot f\left({\rm E}\left(t\right),t\right)\text{d}t\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}m=1,2,\cdots ,n;\text{\hspace{0.17em}}K=1,2,\cdots ,N$ (31)

where
$\stackrel{\xaf}{{F}_{K}}$ is a complex conjugate function of Fourier polynomial F_{k}. The initial condition,
$E\left(0\right)={E}_{I}$ , corresponds to the (initial) conditions of time variable function of unknown of specific problem after variable separation.

By the comparison of spectrum time-history integral with finite differential method, a remarkable advantage of the proposed integral method is that the Courant stability condition, just as in some difference methods, can be totally removed, and the time step size is limited only by the numerical accuracy. Except the above statements, one of the reasons which we choose spectrum time-history integral as the solution philosophy of ODE’s in time of Equation (26) is that the formulation of ODE’s in time often are nonlinear, it usually has better convergence properties to integrate time-dependent evolution equation through spectrum method than finite differences by time-marching.

5. Numerical Results

For the several different construction conditions of rolling compaction ERDs, to assess the validity and feasibility of application of spectral element method to ERDs under construction, a few numerical experiments have been performed using above-formulated methods. The physical properties of damming materials are listed in Table 1.

Corresponding to above-stated several working conditions of compaction ERDs, three realistic engineering applications are addressed below.

5.1. Numerical Illustration of Spectral Method for One-Dimensional Rapid and Layered Embankment of Compacted ERDs

The filling materials in a rolling layer can be approximately considered to be homogeneous when rapid embankment of ERDs, so the calculation of pore-pressure of compacted ERDs can be conducted as one-dimensional problem. In application of spectral element method to calculation of pore-pressure of compacted ERDs, to demonstrate the high numerical accuracy and efficiency for the homogeneous domain only one spectral element is taken and the numbers of Chebyshev interpolation points are 6, 10 and 20 respectively for comparison.

5.1.1. Layered Embankment of Compacted ERDs

According to the initial condition given in Equation (20), the maximum error of the solution of Chebyshev spectral element method compares with ones by 2nd central finite difference methods which are demonstrated in Table 2 at end of a rolling layer.

Table 2 clearly demonstrates that the term number of Chebyshv polynomial expansion reaches 6 or larger Chebyshev spectral element methods can provide quite accurate solution for general engineering problems. In order to validate the preliminary conclusion aforementioned and check out effect of spectral element size on term number of spectral expansion, based on the width range 5 ~ 150 m of filling layer of ERDs from top to bottom, we take the same and different number of spectral elements to calculate pore-pressure respectively and compare their differences which are displayed in Table 3. Table 3 depicts the pore-pressures versus position along the filling layer from downstream to upstream surfaces at the ratio of filling loading $\alpha =21\text{\hspace{0.17em}}\text{N}/\text{a}\cdot {\text{cm}}^{2}$ at time 0.75a.

The physical model is formulated in section 4.1 for present example, i.e. the governing equation is stipulated by Equation (21) and the boundary and initial conditions by Equation (20) where σ =14 N/cm^{2}. The width of rolling soil layer

Table 1. Physical properties of damming materials under construction.

Table 2. Max norm error of one-dimensional pore-pressure of single compacted layer.

Table 3. Pore-pressures of a filling layer versus position coordinates.

is 30 m and height 1 m. Table 3 reveals that, when the numbers of spectral elements (N) are 1, 6 and 30, and each SE have the same interpolate nodes, respectively, it doesn’t appear to be much different one from the other for the pore-pressure solution at midpoint of rolling soil layer.

5.1.2. One-off Rapidly Loading in Embankment of Compacted ERDs

In line with the statement of Section 2.2, we can optionally take a filling layer with 1 meter thick, lying certain depth under the crest as a computational example of one-off embankment of ERDs. For the downstream and upstream surfaces of this filling layer are draining profile, the pore-pressures of the drainage surfaces are equal to zero. The initial pore-pressure in the filling layer could be thought as a constant caused by the weight of overlying soil bodies of ERDs.

Because the pore-pressure of roller compacted earth-rockfilled dams submit to similar formulation of governing Equation (4) and (6), and the range of independent variables of the unknown is limited in the regular domain [−1,1], whether in layered embankment or one-off rapidly loading of compacted ERDs, so the item numbers of Chebyshev polynomials are independent of geometric size of spectral elements.

5.2. Biaxial Drainage Filling of ERDs

For more general situations the pore water will be drained off from the rolling surface, the bottom of rolling soil beds and dam slopes in rolling process of earth-rockfilled dams, so the dissipation of pore-pressure of rolling layers could be thought as two-dimensional problems. In order to invalid the spectral element method for the pre-pressure calculation of biaxial drainage filling of ERDs, at present we only investigate the dissipation of the pore pressure in single filled soil layer, the geometric model of spectral element method calculation is shown in Figure 3.

Figure 3. Spectral element model for biaxial drainage filling layer of ERDs.

For a filling soil layer with 1 m height and 30 m width, it is discretized into eight spectral elements along the width direction in spectral element method calculation, from left to right these spectral elements are labeled by SE1, SE2, etc. in turn. The spectral approximation expansion of the unknown just expand to the fourth term, namely internal node numbers of the spectral element are two. Figure 4 depicts the time-history curve of the pore pressure at midpoint (15, 0.5) of filling layer by spectral and finite element method (30 elements taken).

Since there are no exact solutions used as benchmark for this example, so the above comparison is not conceivable to give two algorithms which is better or worse. But generally speaking, increasing the number of finite element is helpful to improve the calculation precision of FE calculation, for this end one can further compare FEM with SEM which is better than other in the calculation precision by the way of increasing the number of FE.

6. Conclusions

Using Chebyshev spectral element method to calculate pore-water pressure of rolling earth-rockfilled dams during construction period is a new attempt, with the aim of taking advantage of high numerical accuracy and possible application of large spectral elements, for homogeneous block material, of spectral element methods. By numerical examples we also confirmed that, due to application of large spectral element for the uniform materials, the geometrical adaptability of spectral element method is also improved by spectral isoparametric transformation. The details of conclusion can be summarized as follows:

Item number of Chebyshev polynomial of an approximation expansion is proportional to node number of spectral element. When the feature sizes of spectral element have the same orders of magnitude and node number N in one spectral element is bigger than 6 not smaller, it is unlikely to help much to augment the number of nodes to improve spectral precision, but it will greatly increase the amount of calculation.

Chebyshev spectral element method uses Chebyshev polynomials as interpolation function and exponential polynomials are applied in general finite element method; both are based on the same Galerkin projection and similar

Figure 4. Comparison of the pore pressure at midpoint (15, 0.5) of filling layer by spectral and finite element method (30 elements taken).

formulae derivation. Therefore the existing finite element codes can be used in spectral element method after a little transformation. However, spectral element method using a large number of internal nodes in one spectral element, and the pre-processor of finite element can’t handle this situation, so it brings about great barrier for spectral element method applied to general engineering problem.

After spatial discretization of governing differential equation, the first-order evolution equation in time is integrated by spectral method, so the time-history integration precision can be greatly improved and avoid the limitation of integral stability, but the expenditure of spectral calculation gets a little more than finite difference methods.

Numerical examples illustrate that Chebyshev spectral element methods have seemingly the same adaptability as finite element methods for the engineering problems of homogeneous material properties in each spectral element when spectral iso-parametric transformation is employed.

Acknowledgements

The research described in this paper was financially supported by the National Natural Science Foundation of China (Grant No. 51179129).

References

[1] Boyd, J.P. (2000) Chebyshev and Fourier Spectral Methods. DOVER Publications, Inc., Mineola, New York.

[2] Ferronato, M., Castelletto, N. and Gambolati, G. (2010) A Fully Coupled 3-D Mixed Finite Element Model of Biot Consolidation. Journal of Computational Physics, 229, 4813-4830. https://doi.org/10.1016/j.jcp.2010.03.018

[3] Gaspar, F.J., Lisbona, F.J. and Vabishchevich, P.N. (2003) A Finite Difference Analysis of Biot’s Consolidation Model. Applied Numerical Mathematics, 44, 487-506. https://doi.org/10.1016/S0168-9274(02)00190-3

[4] Gu, W.C. (2006) Design and Computation of Earth-Rockfilled Dams. Chinese Press of Building Industry, Beijing. (In Chinese)

[5] Kumar, S. (2010) Object-Oriented Finite Element Programming for Engineering Analysis in C++. Journal of Software, 5, 689-696. https://doi.org/10.4304/jsw.5.7.689-696

[6] Lin, Z. (2003) Design of Compaction Earth-Rockfilled Dams. Huanghe Hydraulic Press, Zhengzhou. (In Chinese)

[7] Michael, S. (2006) Computational Engineering—Introduction to Numerical Methods. Springer, Berlin Heidelberg.

[8] Patera, A.T. (1984) A Spectral Element Method for Fluid Dynamics: Laminar Flow in a Channel Expansion. Journal of Computational Physics, 54, 468-488. https://doi.org/10.1016/0021-9991(84)90128-1

[9] Pozrikidis, C. (2005) Introduction to Finite and Spectral Element Methods Using Matlab. Chapman & Hall/CRC, New York.

[10] Azizi, N., Saadatpour, M.M. and Mahzoon, M. (2012) Using Spectral Element Method for Analyzing Continuous Beams and Bridges Subjected to a Moving Load. Applied Mathematical Modelling, 36, 3580-3592. https://doi.org/10.1016/j.apm.2011.10.019

[11] Liu, Z.L., Menouillard, T. and Belytschko, T. (2011) An XFEM/Spectral Element Method for Dynamic Crack Propagation. International Journal of Fracture, 169, 183-198. https://doi.org/10.1007/s10704-011-9593-y

[12] Raju, G.N., Dutt, P., Kumar, N.K. and Upadhyay, C.S. (2014) Spectral Element Method for Elliptic Equations with Periodic Boundary Conditions. 246, 426-439.

[13] Oliveira, S.P. and Azevedo, J.S. (2014) Spectral Element Approximation of Fredholm Integral Eigenvalue Problems. Journal of Computational and Applied Mathematics, 257, 46-56. https://doi.org/10.1016/j.cam.2013.08.016