Research Article  Open Access
Axisymmetric Consolidation of Unsaturated Soils by Differential Quadrature Method
Abstract
Axisymmetric consolidation in a sand drain foundation is a common problem in foundation engineering. In unsaturated soils, the excess porewater and poreair pressures simultaneously change during the consolidation procedure; and the solutions are not easy to obtain. The present paper uses the differential quadrature method (DQM) for axisymmetric consolidation of unsaturated soils in a sand drain foundation. The radial seepage of sand drain foundation is considered based on the framework of Fredlund’s onedimensional consolidation theory in unsaturated soils. With the use of Darcy’s law and Fick’s law, the polar governing equations of excess poreair and porewater pressures of axisymmetric consolidation are derived. By using DQM, the two governing equations are transformed into two sets of ordinary differential equations. Then the solutions of excess porewater and poreair pressures can be obtained by RongKutta method. The DQM solution can be used to deal with the case of nonuniform initial poreair and porewater distributions. Finally, case studies are presented to investigate the behavior of axisymmetric consolidation of unsaturated soils. The convergence analysis and average degree of consolidation, the settlements in radial and vertical direction, and the effects of different initial excess pore pressure distributions are presented, and discussed in this paper.
1. Introduction
As a common phenomenon of geotechnical engineering, consolidation is a process of decreasing soil volume when soil is subjected to increased stress. Terzaghi [1] established a classical theory for the analysis of onedimensional (1D) consolidation in saturated soils, which is still widely used in engineering practice. But in real cases, soils related to engineering are usually in a state of unsaturation and are defined by more characteristics than saturated soils, such as the coupled dissipation of porewater and poreair pressures during consolidation. In the case of unsaturated soil, the excess porewater pressure and poreair pressure change simultaneously during the consolidation process.
For the consolidation of unsaturated soil, Biot [2] derived a general theory of consolidation for unsaturated soil which had occluded air bubbles. The theory provides a coupling between the magnitude and progress of settlement. Blight [3] presented a consolidation equation for the air phase, where the air was in a dry, rigid, and unsaturated soil. Based on Biot’s theory of consolidation, Scott [4] derived a consolidation equation for unsaturated soils with occluded air bubbles by changing void ratio and degree of saturation. Subsequently, Barden [5] presented an analysis for 1D consolidation of compacted unsaturated clay. Fredlund and Hasan [6] proposed a 1D consolidation theory for unsaturated soils. This formulation is based on two continuity partial differential equations, one for the water phase and the other for the air phase, which have to be solved simultaneously to give water and air pressures at any time and elevation throughout the soil. According to Fredlund’s consolidation theory [7], with the use of Darcy’s law and Fick’s law, Qin et al. [8] obtained a semianalytical solution for consolidation of unsaturated soils with free drainage well. They obtained the excess poreair and porewater pressures and the soil layer settlement in the Laplace transformed domains by utilizing the Bessel function. More recently, Zhou et al. [9] presented a simple analytical solution to 1D consolidation for unsaturated soil, which can be degenerated into Terzaghi consolidation solution for fully saturated condition.
As the consolidation equations of unsaturated soil are coupled, the analytical solutions cannot be easily obtained. Therefore, numerical solutions such as the finite element method (FEM) and the finite difference method (FDM) are traditionally employed. These two methods approximate the partial derivatives of a function at a grid point, only by using a limited number of function values in the vicinity of that grid point. The accuracy and stability of these methods depend on the size of the grid spacing. Compared with these two numerical methods, differential quadrature method (DQM) is a more efficient numerical method for the rapid solution of linear and nonlinear partial differential equations involving one dimension or multiple dimensions. The fundamental idea of DQM is that the derivatives of each node in a continuous function can be expressed as a weighted linear sum of function values at all grid points. Malik and Civan [10] have made a comprehensive comparison for linear and nonlinear convectiondiffusionreaction problems and have shown that the DQM is superior to FEM and FDM in numerical accuracy as well as computational efficiency. Many studies have successfully employed DQM for solving engineering problems [11, 12]. Examples for the applications of DQM in solving complex problems in geotechnical engineering are presented in these literatures [13–16] among others.
In this paper, DQM is introduced into the analysis of axisymmetric consolidation of unsaturated soils in a sand drain foundation. Following Fredlund’s onedimensional consolidation theory of unsaturated soils, with the use of Darcy’s law and Fick’s law, the polar governing equations of excess poreair pressure and porewater pressure of axisymmetric consolidation are obtained. The behaviors of axisymmetric consolidation of unsaturated soils in the sand drain foundation have been analyzed. The convergence analysis and average degree of consolidation, the effects of settlement in radial and vertical directions, the effects of different initial excess poreair and porewater pressure distributions, and the effects of different boundary conditions are presented and discussed.
2. Mathematical Model and Governing Equations
2.1. Mathematical Model for Axisymmetric Consolidation of Unsaturated Soils
As shown in Figure 1, the unsaturated soil layer is described in a schematic diagram and thickness, . Figures 1(a) and 1(b) are the sectional view and vertical view of the unsaturated soil layer, respectively. The radii of soil layer and sand drain are and . A surcharge is applied on the top surface of the soil layer. and are the permeability coefficients of water and air, respectively. The coordinated origin is selected at the center of the top surface and the coordinate is positive downward. In this paper, the vertical prismatic blocks of soil surrounding the sand drains are simulated by cylindrical blocks, of radius . Square pattern and triangular pattern are two common distributions of sand drains. For instance, if square pattern is adopted as shown in Figure 1(c), then . The basic assumptions are the same as those of Fredlund’s onedimensional consolidation theory for unsaturated soils and the other assumptions are listed as follows.(1)The seepage of water and air phases is considered only in radial direction.(2)The top and bottom surfaces are impermeable and the right boundary surface is impermeable or impeded.(3)The permeability coefficients of water and air phases are assumed to be constant.(4)In the process of consolidation, the strain is small enough and the density of water is assumed to be constant.
(a)
(b)
(c)
2.2. Governing Equations
The net flux of water through the soil layer is computed from the volume of water entering and leaving the soil layer within a period of time, with respect to Darcy’s law in the polar coordinate system: where is the initial total volume of the soil, is the net flux of water per unit volume of the soil, is the excess porewater pressure, is the water unit weight, and is the time variable.
The net flux of water per unit volume of the soil can be obtained by differentiating the water phase constitutive relation with respect to time: where is the coefficient of water volume change with respect to a change in the net normal stress , is the coefficient of water volume change with respect to a change in matric suction , and is the excess poreair pressure.
For constant loading, , substituting (1) into (2), the governing equation for the water phase can be written as where and .
The net flux of air through the soil layer is computed from the volume of air entering and leaving the soil layer within a period of time, with respect to Fick’s law in the polar coordinate system: where is the density of air phase, is the net mass rate of air flow per unit volume of the soil, and is the gravitational acceleration.
The total volume change of an unsaturated soil can be assumed to be small during the consolidation process. The volume of air can be related to the volumemass properties of the soil: , where and are the initial porosity and initial degree of saturation before loading. According to Boyle’s law, assuming there is no initial excess air pressure in the soil before loading, we have [17, 18] where is the universal air constant , is the absolute temperature, is the molecular mass of air, is absolute poreair pressure (i.e., ), and is atmospheric pressure.
The derivative of the air phase constitutive relation with respect to time is equal to the net flux of air per unit volume of the soil: where is the coefficient of air volume change with respect to a change in the net normal stress and is the coefficient of air volume change with respect to a change in matric suction .
Substituting (5) into (6), the governing equation for the air phase under constant loading can be written as where and .
Further, two transformed governing equations of water and air phases from (3) and (7) can be obtained: where , , and .
2.3. Boundary and Initial Conditions
As shown in Figure 1, the wall of sand drain is permeable and lateral boundary surface is impermeable or impeded; the boundary conditions for radial consolidation are when or , ; (9) reflects that the lateral boundary surface is impermeable or impeded, respectively. and are the coefficients of permeability for water and air at lateral boundary, respectively. is the thickness of lateral boundary.
The initial condition can be written as where and are the initial porewater and poreair pressure distributions.
3. Differential Quadrature Formulation
DQM is a numerical solution technique for initial and/or boundary value problems, proposed by Bellman et al. [19] and Bellman and Casti [20]. According to DQM, a partial derivative of a function with respect to a variable can be approximated by a weighted linear sum of the function values at given discrete points. Chen et al. [13] employed DQM to solve onedimensional consolidation problems in multilayered soils.
To show the mathematic detail of DQM, consider a function on the domain and the domain is dispersed as points. Then the general differential quadrature approximation of the function at the discrete point is given by where are the weighting coefficient of derivative, .
This paper adopts a method derived by Quan and Chang [21] which uses Lagrange polynomial to determine the weighting coefficients
The soil layer is dispersed in the radial direction and the number of discrete points is . Then, in order to solve the equations conveniently, the local coordinate is introduced. The relationship between local coordinate and integral coordinate is where is the integral coordinate of soil layer and and are the radius values of and point, respectively.
The differential of (14) can be expressed as
The relationship between partial differential of local coordinate and integral coordinate is shown as follows:
Hence, (16) can be approximated by DQM into where and are the weighting coefficient matrices of the first order and second order derivatives, respectively.
Let , , and , where ; substituting (16) into the governing equations (8), one obtains where , is the radius value of point, and and are the porewater and poreair pressures at the point.
The boundary condition and initial condition are also approximated by DQM; (9)–(11) are transformed into
According to (19) and (20), (18) can be rewritten as where , , , and ;
Therefore, the governing equations of water and air are translated into two sets of ordinary differential equations (ODEs). The solutions of ODEs can be obtained by using RongKutta method. Then we can apply the solutions into the following formulas: the average degree of consolidation, the radial displacement, and the vertical settlement.
In unsaturated soils, the average degree of consolidation can be divided into two parts: the average degree of consolidation with respect to water phase and the average degree of consolidation with respect to air phase . To obtain the average degree of consolidation, two formulations are given by Fredlund and Rahardjo [7]. Consider
According to the two stressstate variable approaches [6, 7], volume strain is represented by the following constitutive equation for soil layer:
By integrating (26) with respect to time from to , we get the expression of volumetric strain :
Volumetric strain consists of vertical strain and horizontal strain. For the axisymmetric consolidation, it is assumed in this paper that onethird of the volume strain is contributed from vertical strain and twothirds are from horizontal strain. So the radial displacement and vertical settlement can be obtained, respectively. Consider
4. An Example and Convergence Analysis
In this section, the convergence analysis, the average degree of consolidation, and the effects of displacements in radial and vertical directions are discussed by using a simple example.
As shown in Figure 1, a vertical loading is applied on the top surface of the soil layer. The lateral boundary is considered as impermeable. Considering a uniform initial excess porewater and poreair pressure distribution: , . Other parameters are , , , , , , , , , and . By applying the DQM solution presented above, the results of axisymmetric consolidation of unsaturated soils in a sand drain foundation were obtained.
Tables 1 and 2 list the solutions of excess porewater and poreair pressures at radius at different time factors and different equally spaced discrete points, respectively. The ratio of permeability coefficient is . From Tables 1 and 2, the accuracy of solutions increases when the number of discrete point becomes big. It is obvious that 9 equally spaced grid points are sufficient to obtain the converged results of excess porewater and poreair pressures. Therefore, the number of discrete points is adopted in most case studies.


Figure 2 shows the curves of average degree of consolidation with respect to water and air phases when . The ratio of permeability coefficient is . From Figure 2, it can be observed that the consolidation of air phase and water phase begins when . When and , the consolidation of air phase and water phase is almost finished, respectively. At the early stages, soil consolidation is mainly caused by the dissipation of excess poreair pressure. But in the later stages, soil consolidation is mainly caused by the dissipation of excess porewater pressure.
Figure 3 shows the radial displacement development with time and Figure 4 describes the settlement of top surface () at different radii with time factor . The number of discrete points is . The ratio of permeability coefficient is . From Figure 4, we can see that the soil settles earlier at the points with smaller radius, that is, closer to the sand drain. Figure 5 shows the settlement at different heights of soil layer at radius with time factor . When the height of soil layer is equal to , the settlement is almost zero.
5. Cases of Different Initial and Boundary Conditions
As the initial conditions do not need to be constant in the present DQM solution, it is easily used to analyze nonuniform initial porewater and poreair distribution problems. In this part, four different initial porewater and poreair distributions and some different boundary conditions are considered. A vertical loading is applied on the top surface of soil layer. Other parameters are , , , , , , , , , and .
5.1. Effect of Different Initial PoreWater/Air Pressure Distributions
In this section, in order to make the curves smooth, the number of discrete points is chosen to be . The lateral boundary is considered as impermeable and the ratio of permeability coefficient with respect to air and water is .
Four different initial porewater/air distributions are taken into account. In Figure 6(a), the case of uniform initial poreair and porewater pressures assumes that and . Figure 6(b) presents the linear increasing initial poreair and porewater distributions. The lateralskewed distributions initial porewater/air pressures is shown in Figures 6(c) and 6(d) and they are represented by the following equations: where is the variable which can control the spread or skewness of the curves in Figures 6(c) and 6(d).
(a)
(b)
(c)
(d)
Under the condition shown in Figure 6(a), the dissipation curves of the excess porewater and poreair pressure at different distances from the sand drain are presented in Figures 7 and 8, respectively. The excess poreair and porewater pressures dissipate faster at the points with smaller radius, that is, closer to the sand drain. As the interval between the two curves becomes smaller, from to , it can be seen that the effect of the distance on the dissipation becomes less obvious.
In Figure 6(b), the initial poreair and porewater pressures are assumed not to be constants and to be linear, increasing when the radius increases. The initial poreair pressure is assumed oneeighth of the initial porewater pressure. In other words, from to , the initial porewater pressure will increase from to and the initial poreair pressure will increase from to . Figures 9 and 10 show the dissipation of excess porewater and poreair pressures at different distances from the sand drain under Figure 6(b) distributions.
Comparing Figures 7 and 9, it is found that the excess porewater pressure dissipates completely at almost the same time. With the change of the initial porewater pressure distribution, the time of dissipation with respect to water is not significantly affected. Comparing Figures 8 and 10, it can be seen that the excess poreair pressure dissipates completely at almost the same time. The initial poreair pressure distribution has little influence on the dissipation time of air phase.
Under the conditions of Figures 6(c) and 6(d), the dissipation curves of the excess porewater and poreair pressures for initial distributions where are shown in Figures 11, 12, 13, and 14. From these figures, whatever the initial porewater and poreair pressures are, the excess porewater and poreair pressures dissipate completely at almost the same time. The plateau period is longer at the points with larger radius.
Here, excess pore pressure isochrones for the above four initial conditions are presented in Figures 15 and 16. Figures 15(a) and 16(a) show the excess porewater and poreair pressure isochrones when the initial excess pore pressure is uniform at all radii. In Figures 15(b) and 16(b), the initial excess pore pressure distributions are linearly increasing. The pore pressure isochrones for lateralskewed initial distributions are shown in Figures 15(c), 16(c), 15(d), and 16(d), respectively. From Figures 15(a) and 16(a), as the initial excess pore pressure is constant, the peak of curves in Figure 15(a) approaches to and the peak of curves in Figure 13(a) approaches to . Comparing Figures 15(a) and 15(b), as the initial excess porewater pressure is linearly increasing, the skewness is more evident during early stages of dissipation in Figure 15(b). But in the later stages, the influence of initial excess porewater pressure becomes less significant. So are the isochrones of the excess poreair pressure in Figures 16(a) and 16(b). Figures 15(c), 16(c), 15(d), and 16(d) show the excess porewater and poreair pressure isochrones for the lateralskewed initial distribution when . It is seen that some intersections appear in the isochrones during the early stages. It means that a redistribution of excess pore pressures occurs toward the region of minimal initial excess pore pressure during the early stages. With the time increasing, the isochrones of excess pore pressures become regular.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
5.2. Effect of Different Boundary Conditions
In this part, the effect of different lateral boundary conditions is studied. Uniform initial excess porewater/air pressure distributions and are adopted. Four different ratios are considered, that is, 0.1, 1, 10, and 100 with . Two different drainage conditions are analyzed: lateral surface is impervious and lateral surface is impeded. For the impeded surface, and in (9) are chosen to be .
Figures 17 and 18 shows the dissipation curves for the water and air pressures at radius under different boundary conditions. First, the excess porewater and poreair pressures dissipate faster in impeded surface than that in impermeable surface. Second, is the same in the cases of four values. It is observed from the results that the bigger is, the faster the excess porewater and poreair pressures dissipate. Third, the dissipation curves of the excess poreair and porewater pressures have almost the same shape under different values of . Comparing Figure 17 with Figure 18, when the excess poreair pressure dissipates almost completely, the dissipation of excess porewater pressure enters into a plateau period. The bigger is, the earlier the excess poreair pressure dissipates completely. Therefore, it is shown that the value of affects the dissipation of porewater pressure during the consolidation process.
6. Conclusion
This paper obtains a general solution to the axisymmetric consolidation of unsaturated soils by using differential quadrature method (DQM) based on Fredlund’s onedimensional consolidation theory for unsaturated soil. With the use of DQM, the two governing equations are transformed into two sets of ordinary differential equations. The solutions to the transformed differential equations are obtained by RongKutta method under different initial and boundary conditions. A case study has been presented for the analysis of unsaturated consolidation in a sand drain foundation. The convergence analysis, the average degree of consolidation for water and air phases, the settlement in radial and vertical directions, effects of different initial excess poreair and porewater pressure distributions, and effects of different boundary conditions have been presented. From the analyses, some conclusions can be drawn: (a) through compiling the programs, it is found that the DQ solution delivers accurate and stable results for unsaturated soil consolidation. Due to the uniform matrix structure of the DQ formulas, it is easy to compile computer programs when considering complicated conditions such as nonuniform porewater/air distribution conditions; (b) soil consolidation is mainly caused by the dissipation of excess poreair pressure in the early stages and by the dissipation of excess porewater pressure in the later stages; (c) the bigger the ratio of the permeability coefficients for water and air is, the faster the excess porewater and poreair pressures dissipate; (d) The excess poreair and porewater pressures dissipate faster at the points closer to the sand drain; (e) the initial distribution has some effects on the early consolidation process of unsaturated soil and boundary condition has a significant effect on the whole consolidation process.
Notation
:  The consolidation coefficient of water 
:  The consolidation coefficient of air 
:  The weighting coefficient 
:  Volumetric strain 
:  The gravitational acceleration 
:  The thickness of soil layer 
:  The permeability coefficients of air 
:  The coefficients of permeability for water at lateral boundary 
:  The permeability coefficients of water 
:  The coefficients of permeability for air at lateral boundary 
:  The coefficient of air volume change with respect to a change in the net normal stress 
:  The coefficient of water volume change with respect to a change in the net normal stress 
:  The coefficient of air volume change with respect to a change in matric suction 
:  The coefficient of water volume change with respect to a change in matric suction 
:  The initial porosity 
:  The discrete point 
:  A vacuum pressure 
:  The density of air phase 
:  The vertical loading 
:  The thickness of lateral boundary 
:  The radius of soil layer 
:  The radius of sand drain 
:  The radius value of th point 
:  The unit weight of water phase 
:  The universal air constant 
:  The distance between two sand drains 
:  The initial degree of saturation 
:  Settlement 
:  The radial displacement 
:  Time 
:  The time factor 
:  The excess poreair pressure 
:  Atmospheric pressure 
:  The initial excess poreair pressure 
:  The excess porewater pressure 
:  The initial excess porewater pressure 
:  The average degree of consolidation with respect to air 
:  The average degree of consolidation with respect to water 
:  The air volume of soil layer 
:  The water volume of soil layer 
:  Initial total volume of the soil. 
Acknowledgments
The author gratefully acknowledges the financial support from the Macau Science and Technology Development Fund (Grant no. FDCT/011/2013/A1) and the University of Macau Research Fund (Grants nos. MYRG189(Y2L3)FST11ZWH and MYRG067(Y2L2)FST12ZWH).
References
 K. Terzaghi, Theoretical Soil Mechanics, Wiley, New York, NY, USA, 1943.
 M. A. Biot, “General theory of threedimensional consolidation,” Journal of Applied Physics, vol. 12, no. 2, pp. 155–164, 1941. View at: Publisher Site  Google Scholar
 G. E. Blight, Strength and consolidation characteristics of compacted soils [Ph.D. thesis], University of London, London, UK, 1961.
 R. F. Scott, Principles of Soil Mechanics, Addison Wesley Publishing Company, Boston, Mass, USA, 1963.
 L. Barden, “Consolidation of compacted and unsaturated clay,” Geotechnique, vol. 15, no. 3, pp. 267–286, 1965. View at: Publisher Site  Google Scholar
 D. G. Fredlund and J. U. Hasan, “Onedimensional consolidation theory: unsaturated soils,” Canadian Geotechnical Journal, vol. 16, no. 3, pp. 521–531, 1979. View at: Google Scholar
 D. G. Fredlund and H. Rahardjo, Soil Mechanics for Unsaturated Soils, John Wiley and Son, 1993.
 A. Qin, D. Sun, L. Yang, and Y. Weng, “A semianalytical solution to consolidation of unsaturated soils with the free drainage well,” Computers and Geotechnics, vol. 37, no. 78, pp. 867–875, 2010. View at: Publisher Site  Google Scholar
 W. H. Zhou, L. S. Zhao, and X. B. Li, “A simple analytical solution to onedimensional consolidation for unsaturated soil,” International Journal for Numerical and Analytical Methods in Geomechanics, 2013. View at: Publisher Site  Google Scholar
 M. Malik and F. Civan, “A comparative study of differential quadrature and cubature methods visàvis some conventional techniques in context of convectiondiffusionreaction problems,” Chemical Engineering Science, vol. 50, no. 3, pp. 531–547, 1995. View at: Google Scholar
 P. Malekzadeh, M. R. Golbahar Haghighi, and M. M. Atashi, “Outofplane free vibration of functionally graded circular curved beams in thermal environment,” Composite Structures, vol. 92, no. 2, pp. 541–552, 2010. View at: Publisher Site  Google Scholar
 P. Malekzadeh, M. R. G. Haghighi, and M. M. Atashi, “Free vibration analysis of elastically supported functionally raded annular plates subjected to thermal environment,” Meccanica, vol. 46, no. 5, pp. 893–913, 2011. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. P. Chen, W. H. Zhou, H. Z. Wang, and Y. M. Chen, “Onedimensional nonlinear consolidation of multilayered soil by differential quadrature method,” Computers and Geotechnics, vol. 32, no. 5, pp. 358–369, 2005. View at: Publisher Site  Google Scholar
 R. P. Chen, W. H. Zhou, H. Z. Wang, and Y. M. Chen, “Study on onedimensional nonlinear consolidation of multilayered soil using differential quadrature method,” Chinese Journal of Computational Mechanics, vol. 22, no. 3, pp. 310–315, 2005. View at: Google Scholar
 H. Z. Wang, R. P. Chen, W. H. Zhou, and Y. M. Chen, “Computation of 1D nonlinear consolidation in doublelayer foundation by using differential quadrature method,” Journal of Hydraulic Engineering, vol. 4, pp. 8–14, 2004. View at: Publisher Site  Google Scholar
 W. H. Zhou and L. S. Zhao, “Onedimensional consolidation of unsaturated soil subjected to timedependent loading under various initial and boundary conditions,” ASCE International Journal of Geomechanics, 2013. View at: Publisher Site  Google Scholar
 E. Conte, “Consolidation analysis for unsaturated soils,” Canadian Geotechnical Journal, vol. 41, no. 4, pp. 599–612, 2004. View at: Publisher Site  Google Scholar
 E. Conte, “Plane strain and axially symmetric consolidation in unsaturated soils,” International Journal of Geomechanics, vol. 6, no. 2, pp. 131–135, 2006. View at: Publisher Site  Google Scholar
 R. Bellman, B. G. Kashef, and J. Casti, “Differential quadrature: a technique for the rapid solution of nonlinear partial differential equations,” Journal of Computational Physics, vol. 10, no. 1, pp. 40–52, 1972. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. Bellman and J. Casti, “Differential quadrature and longterm integration,” Journal of Mathematical Analysis and Applications, vol. 34, no. 2, pp. 235–238, 1971. View at: Publisher Site  Google Scholar  MathSciNet
 J. R. Quan and C. T. Chang, “New insights in solving distributed system equations by the quadrature method. I: analysis,” Computers and Chemical Engineering, vol. 13, no. 7, pp. 779–788, 1989. View at: Google Scholar
Copyright
Copyright © 2013 WanHuan Zhou. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.