An Alternative One-Step Computation Approach for Computing Thermal Stress of Asphalt Mixture: the Laplace Transformation

새로운 아스팔트 혼합물의 저온응력 계산 기법에 대한 고찰: 라플라스 변환

Abstract

Computing low temperature performance of asphalt mixture is one of the important tasks especially for cold regions. It is well known that experimental creep testing work is needed for computation of thermal stress and critical cracking temperature of given asphalt mixture. Thermal stress is conventionally computed through two steps of computation. First, the relaxation modulus is generated thorough the inter-conversion of the experimental creep stiffness data through the application of Hopkins and Hamming's algorithm. Secondly, thermal stress is numerically estimated solving the convolution integral. In this paper, one-step thermal stress computation methodology based on the Laplace transformation is introduced. After the extensive experimental works and comparisons of two different computation approaches, it is found that Laplace transformation application provides reliable computation results compared to the conventional approach: using two step computation with Hopkins and Hamming's algorithm.

겨울철 발생하는 아스팔트 혼합물의 저온균열의 정량적 분석 및 평가를 위해서는 해당 아스팔트 혼합물의 저온응력이 반드시 계산되어야 하며, 이는 현재 대한민국, 미국 북부 및 캐나다 지역에서 포장 유지관리, 설계에 있어서 매우 중요한 사항 중 하나이다. 일반적으로 아스팔트 혼합물의 저온응력은 크리프 시험과 시간중첩이론을 바탕으로 계산되며 전통적으로 두 수학적 단계를 통해 계산된다. 우선 수학적, 수치적 변환과정(홉킨스-해밍 알고리즘)을 통해 크리프-강성응력에서 이완응력이 계산된다. 다음으로 이완응력 지배곡선을 구현한 후 회선적분의 수치해석적 접근을 통해 아스팔트 혼합물의 저온응력이 최종적으로 계산된다. 상기의 과정은 복잡하며, 시간이 오래 걸리는 단점이 있다. 이번 논문에서는 보다 간편한 라플라스 변환을 통해 해당 아스팔트 혼합물의 저온응력을 계산하였으며, 이의 결과를 전통적 계산 기법과 비교, 분석하였다. 결론적으로 새로이 제안된 라플라스 변환 기법은 보다 아스팔트 혼합물의 저온응력을 효과적, 효율적으로 계산할 수 있음이 발견되었다.

1. Introduction

Computing and estimating low temperature performance of asphalt mixture is one of the crucial tasks especially for South Korea, northern U.S. and Canada (Basu, 2002; Velasquez et al., 2009; Moon, 2010; Marasteanu et al., 2009; Moon et. al., 2013; Cannone Falchetto et al., 2018). In many pavement agencies in various countries, the low temperature performance of asphalt pavement layer is measured through the resistance level (and/or value) against longitudinal thermal cracking: the thermal stress and corresponding critical cracking temperature (Moon, 2010; Marasteanu et al., 2009; Moon 2012; Cannone Falchetto et al., 2014a; Cannone Falchetto et al., 2014b; Moon et al., 2014a; Moon et al., 2014b; Cannone Falchetto et al., 2018). Moreover, it is well known that thermal stress and corresponding critical cracking temperature of given asphalt material can be numerically computed (and/or measured) by means of low temperature creep testing (Park and Kim, 1999; Di Benedetto et al., 2009; Moon, 2010; Marasteanu et al., 2009; Moon 2012; Cannone Falchetto et al., 2014a; Cannone Falchetto et al., 2014b; Moon et al., 2014a; Moon et al., 2014b; Cannone Falchetto et al., 2018).

The conventional computation process of thermal stress is consisted of three steps: first, the relaxation modulus: E(t), is generated through the inter-conversion process of the experimental creep compliance data: D(t), by application of Hopkins and Hamming’s algorithm (1967). Secondly, E(t) master curve is generated through various mathematical models and corresponding parameters are computed. Thirdly, thermal stress: σ(T), is numerically computed by solving the E(t) convolution integral (Marasteanu and Anderson, 1999; Basu, 2002; Moon, 2010; Marasteanu et al., 2009; Moon 2012; Cannone Falchetto and Moon, 2015). However, this σ(T) computation approach is relatively complex for many pavement agency employees and researchers to perform. In this paper, an alternative and simple σ(T) computation approach using the Laplace transformation is introduced. As an experimental work, Bending Beam Rheometer (BBR) mixture creep test (Marasteanu et al., 2009; Moon, 2010; Moon, 2012) was performed with 3 different asphalt mixtures. Finally, all the computed σ(T) results by using two different computation approach were compared graphically.

2. Experimental Works

2.1 Asphalt Material Preparations

Three different asphalt mixtures were prepared and tested in the asphalt pavement laboratory at Korea Expressway Corporation (KEC). All the mixtures were compacted with a Superpave Gyratory Compactor (SGC) (Marasteanu et al., 2009) to a target air voids content of 7 % (AASHTO M320-10, 2010) with Nominal Maximum Aggregate Size (NMAS) of 12.5 mm. In addition, all the mixtures were designed with a Styrene Butadiene Styrene (SBS) modified asphalt binder containing Performance Grade (PG) 76-34 (AASHTO M320-10, 2010). The overall design binder content was set to 4.2~4.8 %. Two different amounts of Reclaimed Asphalt Pavement (RAP): 15 % and 25 %, were included in the mix design. All the specimens were long term aged based on current AASHTO specification mentioned elsewhere (AASHTO R028-12, 2012). Schematic information of the materials used in this paper is shown in Table 1.

Table 1. Asphalt Mixtures

Note: NMAS=Nominal Maximum Aggregate Size, VMA = Voids between Mineral Aggregates, V_be = Volume of effective asphalt binder, Air = Air void, Asphalt = Portion of asphalt binder, For mixtures 1 to 3, maximum filler content is around 1.4~1.6%, VMA=V_be+Air

2.2 Bending Beam Rheometer (BBR) Mixture Creep Test

To compute σ(T), Bending Beam Rheometer (BBR) mixture creep experimental works (Marasteanu et al., 2009) were performed on thin asphalt binder mixture beams (length : 102.05 mm width : 12.70.5 mm height : 6.250.5 mm) (see Fig. 1) for 1,000 testing seconds with applying 4,000~6,000 mN of constant load amount.

Fig. 1. BBR Binder Mixture Test Specimen (a) and Testing Equipment Setup (b)

Similar to the BBR low temperature creep testing on asphalt binder currently mentioned as testing specification in AASHTO(AASHTO T313-06, 2006), the creep stiffness: S(t) can be computed based on Bernoulli-Euler beam theory in BBR creep mixture test as (Moon, 2010; Moon, 2012):

\(S(t)=\frac{\sigma}{\varepsilon(t)}=\frac{P \cdot l^{3}}{4 \cdot b \cdot h^{3} \cdot \delta(t)}=\frac{1}{D(t)}\)       (1)

In Eq. (1), D(t), σ, ε(t), P, l, b, h, δ(t) and t mean creep compliance, maximum bending stress, bending strain, applied load, length, width, height, deflection in the mid span and testing time, respectively. To generate relaxation modulus: E(t), master curves and to compute thermal stress: σ(T), the BBR mixture creep tests were performed at two different temperatures: low PG + 10°C and low (PG +10°C)+12°C (i.e. -24: reference temperature and -12°C: higher temperature). Three replicated were used at each testing condition therefore, total six mixtures were used for one set of prepared asphalt mixture in this study.

3. Two Different Approaches for Computing Thermal Stress of Asphalt Mixture

3.1 Traditional Approach using Hopkins & Hamming’s Algorithm (1967) and CAM Model (1999)

It is well known that the values of D(t) and E(t), for a viscoelastic material are correlated through the convolution integral as (Findley et al., 1976; Ferry, 1980; Ebrahimi et al., 2014):

\(\int_{0}^{t} E(t-\tau) \cdot D(t) d \tau=t \text { or } \int_{0}^{t} E(t) \cdot D(t-\tau) d \tau=t\)       (2)

Based on Eq. (2), E(t) was generated numerically by applying Hopkins and Hamming’s algorithm (1967) in this paper as following computation steps:

1) Select time interval as:

\(t_{0}=0, t_{1}=1, t_{2}=2, \cdots t_{1,000}=1,000\)       (3)

2) Express the integral form of D(t) and compute f(t) as (see Eq. (4)):

\(\begin{array}{l} f(t)=\int_{0}^{t} D(t) d t \Rightarrow f\left(t_{n+1}\right)= \\ f\left(t_{n}\right)+\frac{1}{2} \cdot\left(D\left(t_{n+1}\right)+D\left(t_{n}\right)\right) \cdot\left(t_{n+1}-t_{n}\right) \end{array}\)       (4)

3) Modify Eq. (2) based on Eq. (4) then express in terms of time: t, as (see Eq. (5)):

\(\begin{array}{l} t_{n+1}=\int_{0}^{t_{n}+1} E\left(t^{\prime}\right) \cdot D\left(t_{n+1}-t^{\prime}\right) d t^{\prime}= \\ \sum_{i=0}^{n} \int_{t_{i}}^{t_{i}+1} E\left(t^{\prime}\right) \cdot D\left(t_{n+1}-t^{\prime}\right) d t^{\prime} \end{array}\)       (5)

4) Approximate each element of the integral in Eq. (5) and rewrite Eq. (6) as (see Eq. (8)):

\(\begin{array}{l} \int_{t_{i}}^{t_{i}+1} E\left(t^{\prime}\right) \cdot D\left(t_{n+1}-t^{\prime}\right) d t^{\prime}= \\ -E\left(t_{i+1 / 2}\right) \cdot\left[f\left(t_{n+1}-t_{t+1}\right)-f\left(t_{n+1}-t_{i}\right)\right] \end{array}\)       (6)

where interval time can be expressed as:

\(t_{i+\frac{1}{2}}=\frac{1}{2} \cdot\left(t_{i+1}+t_{i}\right)\)       (7)

\(\begin{aligned} &t_{n+1}=-\sum_{i=0}^{n-1} E\left(t_{i+1 / 2}\right) \cdot\left[f\left(t_{n+1}-t_{i+1}\right)-f\left(t_{n+1}-t_{i}\right)\right]\\ &+E\left(t_{n+1 / 2}\right) \cdot f\left(t_{n+1}-t_{n}\right) \end{aligned}\)       (8)

5) Derive E(t) from Eq. (8) by solving for E(t_n+1/2) :

\(\begin{aligned} &E\left(t_{n+1 / 2}\right)=\\ &\frac{t_{n+1}-\sum_{i=0}^{n-1} E\left(t_{i+1 / 2}\right) \cdot\left[f\left(t_{n+1}-t_{i}\right)-f\left(t_{n+1}-t_{i+1}\right)\right]}{f\left(t_{n+1}-t_{n}\right)}=\\ &\frac{t_{n+1}-\sum_{i=0}^{n-1} E\left(t_{i+1 / 2}\right) \cdot\left[f\left(t_{i+1}\right)-f\left(t_{i}\right)\right]}{f\left(t_{n+1}\right)-f\left(t_{n}\right)} \end{aligned}\)       (9)

where, \(f\left(t_{0}\right)=0, E\left(t_{0}\right)=0, E\left(t_{1}\right)=\frac{t_{1}}{f\left(t_{1}\right)}\)       (10)

After computation step (5), E(t) master curve was generated based on Christensen, Anderson and Marasteanu (CAM) model (Marasteanu and Anderson, 1999) on BBR experimental data (see Eqs. (11) and (12): Shift factor) as can be seen in Eq. (11):

\(E(t)=E_{g} \cdot\left[1+\left(\frac{t}{t_{c}}\right)^{v}\right]^{-\frac{w}{v}}\)       (11)

\(a_{t}=10^{C_{1}+C_{2} \cdot T}\)       (12)

In Eq. (11), E_g means glassy modulus (i.e. approximately 30~35 GPa, Moon, 2010; Moon et al., 2013) and t_c, v, w, C₁, C₂ represent fitting parameters. Finally, thermal stress: σ(T), was calculated numerically by solving the following convolution integral:

\(\begin{array}{l} \sigma(\xi)=\int_{-\infty}^{\xi} \frac{d \varepsilon\left(\xi^{\prime}\right)}{d \xi^{\prime}} \cdot E\left(\xi-\xi^{\prime}\right) d \xi^{\prime}= \\ \int_{-\infty}^{t} \frac{d(\alpha \cdot \Delta T)}{d t^{\prime}} \cdot E\left(\xi(t)-\xi^{\prime}(t)\right) d t^{\prime} \end{array}\)       (13)

\(\varepsilon(t)=\alpha \cdot \Delta T\)       (14)

In Eq. (14), α means the coefficient of thermal expansion or contraction; in this study it is assumed as 30.2810^-6 (Moon, 2010; Moon et al., 2013). The σ(T) was computed at temperature ranged from 22 to -40°C (temperature interval was 0.5°C) with consideration of 2°C/hour of material cooling rate. Moreover, 24 Gauss points integration (points used for Gauss-Legendre numerical integration) was considered to solve Eq. (12) numerically (Moon, 2010; Moon et al., 2013, Moon et al., 2014b; Moon et al., 2014b). Detailed about the computation process is mentioned elsewhere (Moon, 2010).

3.2 Alternative Approach using Laplace Transformation

As an alternative computation approach for computing σ(T), Laplace transformation was used in this paper. A brief summary of the computation procedure is presented as following steps:

1) Determine the shift factor, a_T, and hence, parameters C₁ and C₂ in Eq. (12) based on the BBR test results (i.e. LogD(t) versus Log(t) curve). Then re-express Eq. (12) as:

\(\begin{array}{l} a_{T}=10^{C_{1}+C_{2} \cdot T}=10^{C_{1}+C_{2} \cdot\left(T_{i}-C_{0} \cdot t\right)} \\ =10^{\left(C_{1}+C_{2} \cdot T_{i}\right)-C_{2} \cdot C_{0} \cdot t}=10^{C_{3}+C_{4} \cdot t}=A_{0} \cdot 10^{C_{4} \cdot t} \end{array}\)       (15)

2) Construct the reduce time based master curve of D(ξ)based on BBR experimental data at two temperature conditions: low PG+10°C and low(PG+10)+12°C, as (see Eq. (16)):

\(D(\xi)=a \cdot \xi^{b}+c \cdot \xi^{d}\)       (16)

where: a, b, c and d = Fitting parameters, and

\(\xi=\frac{t}{a_{T}}=\) Reduced time (at reference temperature T_i = 22°C).

3) Based on Eq. (13), the relationship between strain and thermal stress can be expressed as:

\(\varepsilon_{t}=\int_{0}^{\xi} D\left(\xi-\xi^{\prime}\right) \frac{\partial \sigma}{\partial \xi^{\prime}} d \xi^{\prime}+\int_{0}^{\xi} \alpha\left(\xi-\xi^{\prime}\right) \frac{\partial(\Delta T)}{\partial \xi^{\prime}} d \xi^{\prime}=0\)       (17)

By applying the Laplace transformation, Eq. (17) can be rewritten as:

\(s \cdot \bar{D}(s) \cdot \bar{\sigma}(s)+s \cdot \bar{\alpha}(s) \cdot \Delta \bar{T}(s)=0\)       (18)

\(\bar{\sigma}(s)=-\frac{\bar{\alpha}(s) \cdot \Delta \bar{T}(s)}{\bar{D}(s)}\)       (19)

In Eqs. (18) and (19), parameters α and T can be expressed as:

\(\alpha(t)=\alpha_{0} \Rightarrow \bar{\alpha}(s)=\frac{\alpha_{0}}{s}\)       (20)

\(\Delta T=-C_{0} \cdot t \Rightarrow \Delta \bar{T}(s)=-\frac{C_{0}}{s^{2}}\)       (21)

4) Transform (i.e. inverse transformation) Eq. (19) from Laplace domain to reduced time domain (ξ) using Stehfest algorithm (D.Amore et al., 2018; Stehfest 1970; Villinger H., 1985) then the results of numerically generated σ(ξ) can then be fitted with the following power function (see Eq. (22)):

\(\sigma(\xi)=a+b \cdot \xi^{c}\)       (22)

Finally, thermal stress can be computed by converting the thermal stress function from the reduced time domain, σ(ξ) (Eq. (22)) to the actual time domain, σ(T), with application of Eq. (23):

\(\begin{aligned} &\xi=\int_{0}^{t} \frac{d t^{\prime}}{a_{T}\left(T\left(t^{\prime}\right)\right]}=\int_{0}^{t} \frac{d t^{\prime}}{A_{0} \cdot 10^{C_{i} \cdot t^{\prime}}}=\frac{1}{A_{0}} \int_{0}^{t} 10^{-C_{4} \cdot t^{\prime}} d t^{\prime}\\ &=\frac{1}{A_{0}} \cdot\left[-\frac{10^{-C_{t} t^{r}}}{C_{4} \cdot \ln 10}\right]_{0}^{t}=\frac{1}{A_{0}} \cdot\left[\frac{1}{C_{4} \cdot \ln 10} \cdot\left(1-10^{-C_{i} t}\right)\right]\\ &=\frac{1}{A_{0} \cdot C_{4} \cdot \ln 10} \cdot\left[1-10^{-C_{i} t}\right]=A_{1} \cdot\left[1-10^{-C_{i} t}\right] \end{aligned}\)       (23)

\(A_{0}=10^{C_{3}}=10^{C_{1}+C_{2} \cdot T_{i}} \text { and } C_{4}=-C_{2} \cdot C_{0}\)       (24)

Based on Eqs. (20), (21), (22) and (23), thermal stress was computed at temperature ranged from 22 to -40°C with 0.5°C of temperature interval similar to the previous section. Then the computed results of σ(T) from section 3.1 and 3.2 were visually compared.

4. Data Analysis (Comparison of Thermal Stress Results)

From Eqs. (2) to (23), values of σ(T) were computed with using two different computation approaches: traditional Hopkins and Hamming’s algorithm and Laplace transformation. All the results are presented in Table 2 and Figs. 2 and 3.

Table 2. Comparison Results of Computed Thermal Stress [MPa] (Averaged Results)

Note: M1(H): Mixture 1, Hopkins and Hamming’s algorithm, M1(L): Mixture 1, Laplace transformation

M1(D): Mixture 1, Differences = ABS[(M(H)-M(L))/M(H)], %

Fig. 2. Thermal Stress Comparison Results

Fig. 3. Thermal Stress Comparison Results

From the results in Table 2 and Figs. 2 to 3, it can be found that dramatically higher σ(T°C) results were computed with addition of RAP (e.g. 15 % and 25 %). This data generation trend can support the previous findings that addition of RAP provides stiffer characteristics at low temperature compared to the conventional asphalt mixture.

It also needs to mentioned that Laplace transformation provided reasonable and reliable σ(T°C) computation results compared to the conventional σ(T°C) computation method: using Hopkins & Hamming’s algorithm and CAM model (1999). Below the 0°C degree the σ(T°C) result differences between two computation approaches decreased significantly down to 0.8~2.1 % (see Fig. 3). In addition, almost identical σ(T°C) computation results were observed visually between Laplace transformation and Hopkins & Hamming’s algorithm in all different types of asphalt mixture. Based on the findings it can be said the one-step Laplace transformation approach can be a successful alternative for computing low temperature cracking resistant performance of given asphalt materials. In addition, more experimental works and mathematical parametric studies are needed to further narrow the computation results between two different approaches.

5. Summary and Conclusion

In this paper, two different computation approaches: Hopkins & Hamming’s algorithm and Laplace transformation, for obtaining thermal stress of three different asphalt mixtures were investigated. As an experimental work, newly developed and modified BBR mixture creep test was performed with various types of asphalt mixture including different portion of RAP. Based on the graphical and numerical comparison some crucial findings, results and limitations can be derived.

(1) It can be said that applying one-step Laplace transformation approach, can successfully be an alternative method for computing thermal stress compared to the conventional and sophisticated computation approach. Based on this approach, it is expected that more university and pavement agency employees are able to compute thermal stress results of asphalt material in more convenient manner than before.

(2) However, only three limited asphalt mixtures and BBR mixture creep testing were considered in this study; therefore, additional extensive experimental work with wider range of asphalt materials and various testing approaches (i.e. IDT creep, IDT strength test) is highly recommended to further verify the findings acquired in this study.