Transient deformation of a poroelastic channel bed
P. C. Hsieh
a,y, W. P. Shih
band L. H. Huang
ba
Department of Soil and Water Conservation, National Chung-Hsing University, Taichung 40227, Taiwan, ROC
bDepartment of Civil Engineering, National Taiwan University, Taipei 106, Taiwan, ROC
SUMMARY
The coupled transient response of a poroelastic bed form due to stream flow and non-linear water waves is investigated numerically. The theory of potential flow is applied to channel flow while Biot’s theory of poroelasticity (J. Appl. Phys. 1962; 33(4):1482) is adopted to deal with the deformable porous bed. A boundary-fitted co-ordinate system is used to calculate the variation in the bed form. The result of a simple periodic wave form over a soft poroelastic bed agrees well with the analytical solution of Hsieh et al. (J. Eng. Mech., ASCE 2000; 126(10):1064). However, due to the rapidly damping second dilatational wave inside the soft poroelastic bed, the solution for transient bed form near the interface is not easy to compute accurately. In order to overcome this difficulty, a simplified numerical model based on the boundary layer correction concept proposed by Hsieh et al. (2000) is established, which neglects Darcy’s terms. The transient deformation of an irregular poroelastic bed that includes a trench and a downward step at the channel bed is simulated successfully. Copyright # 2002 John Wiley & Sons, Ltd.
KEY WORDS: Transient deformation; soft poroelastic bed; boundary layer correction concept; boundary-fitted co-ordinate transformation
1. INTRODUCTION
The problem of flow past a porous deformable bed is always an important topic to hydraulic and ocean engineers. In addition to the complexity of interrelated mechanics between homogeneous water flow and porous deformable bed, the interaction between the skeleton and the pore-saturated fluid of the porous bed is even more complicated. For these subjects, some researches based on experiments [1–3] or simplified theories, which may apply empirical mass transport formulas [4–7] to avoid the analysis of complicated deformable porous media mechanics, were published to explain the physical behaviour. In addition, Darwin [8] conducted experiments on sand-ripple caused by oscillatory bed motion and concluded that the eddies induced by a series of vortex acting on the sandy bed would probably render unstable ripple into stable state. Bagnold [9] explained that the instability of dune in the desert was due to the impact
Received 12 June 2001
yCorrespondence to: P. C. Hsieh, Department of Soil and Water Conservation, National Chung-Hsing University,
Taichung 40227, Taiwan, R.O.C. E-mail: [email protected]
Contract/grant sponsor: National Science Council of the Republic of China; contract/grant numbers: NSC 86-2621-E002-016 and NSC 88-2611-E-002-028.
of sand by wind. Exner [10] established a differential erosion equation for two-dimensional flow, which indicated that the change in bed elevation was due to dilatational variation of bottom velocity. Anderson [11] applied Exner’s [10] erosion equation to explain the mechanism of ripple formation. However, the mechanism of dune and antidune formation is somewhat different from that of ripple. Kennedy [12] applied an empirical sediment transport formula to govern the continuity of the porous bed. Also, he used instability analysis of potential theory to obtain his famous discovery about the dune and antidune formations in alluvial channels. Unfortunately, owing to the constraint of instability analysis, Kennedy [12] could only find the dominant wavelength of stable bed forms instead of the whole bed form.
In fact, fluid within porous material interacting with a deforming solid skeleton is a more complicated two-phase flow problem for a realistic analysis. Biot [13] developed a poroelasticity theory to describe elastic waves in a fluid saturated porous solid. Mei and Foda [14] proposed a boundary layer correction to simplify the analysis; however, their approach was based on physical intuition but without systematic analysis. By treating the bed as a poroelastic material, Huang and Song [15] solved the problem of oscillatory linear water waves interacting with a deformable bed by using three decoupled Helmoltz equations derived by Huang and Chwang [16]. In their solution, five non-dimensional parameters were derived. Chen et al. [17] also applied Huang and Chwang’s [16] approach along with the conventional Stokes expansion of
deep water waves based on the perturbation parameter e1¼ k0a ðk0and a are the wave number
and the amplitude of incoming water wave, respectively) to investigate the dynamic response of permeable bed material to non-linear water waves. They found that the conventional Stokes expansion is only valid for hard poroelastic bed material but invalid for soft materials even though the Ursell parameter is small. (Here ‘hard’ and ‘soft’ are defined according to a parameter P; a Mach number of the second longitudinal wave inside the poroelastic material, in their paper.) Huang and Chiang [18] simulated stable bed forms under a constant current and linear oscillatory water waves by an approach similar to that of Huang and Song [15]. In Huang and Chiang [18], various bed forms such as dune, antidune and flat bed under linear oscillatory water wave accompanied with a constant current were obtained. The ambiguous lagged distance d; introduced in Kenney [12], which has bothered researchers in river mechanics for many years was clarified. Huang et al. [19] and Hsieh et al. [20] used a two-parameter ðe1; e2¼ k0=k2; k2 is
the wave number of the second longitudinal wave inside the poroelastic material) perturbation expansion to discuss the dynamic response of soft poroelastic bed to linear and non-linear oscillatory water waves by a boundary layer correction approach.
The aforementioned studies are mostly focused on the stable bed form problem, however, the transient response of a constant current flowing over a semi-infinite soft poroelastic deforming bed along with non-linear water waves (Figure 1) remains to be analysed. In the present study, a finite difference numerical model based on the poroelastic theory is proposed to discuss the transient bed form.
2. FORMULATION
This study analyses the problem of a non-linear water wave along with a uniform channel flow over a poroelastic deformable bed, as defined schematically in Figure 1. Region 1 is homogeneous water, considered as potential flow, while region 2 is a semi-infinite porous medium saturated with water governed by Biot’s theory of poroelasticity [21]. The co-ordinates
of region 1 range from y ¼ gðx; tÞ to y ¼ h þ zðx; tÞ; and region 2 from y ¼ gðx; tÞ to y ! 1: The symbols z and g represent the displacements of waves from the mean free surface ðy ¼ hÞ and mean bed interface ðy ¼ 0Þ; respectively.
2.1. Boundary value problem
Assuming that the homogeneous channel flow is a potential flow, the flow velocity can be written as a given constant current ðU Þ in the x direction (i.e. U
%
ex) plus a perturbed velocity
% u1:
Thus the perturbed velocity potential f1 satisfies
rf1¼
%
u1 ð1Þ
The equations of continuity and momentum in terms of velocity potential f1 become
r2f 1¼ 0 ð2Þ @f1 @t þ 1 2 2U @f1 @x þ @f1 @x 2 þ @f1 @y 2 " # þP1 rf ¼ 0 ð3Þ
where P1 is the perturbed pressure, and rf is the water density.
Referring to Biot’s theory of poroelasticity [13], the linear momentum equations of solid skeleton and fluid for the porous bed may be written as
r %% s ¼ r11@ 2 % d @t2 þ r12 @2 % D @t2 þ b @ % d @t @ % D @t ð4Þ r %% S ¼ r12@ 2 % d @t2 þ r22 @2 % D @t2 b @ % d @t @ % D @t ð5Þ with %% s ¼ %% t ð1 n0ÞP2 %% I ð6Þ %% t ¼ 2G %% e þ lðr % dÞ %% I ð7Þ
%% e ¼1 2½rd þ ðr% dÞ% t ð8Þ %% S ¼ n0P2 %% I ð9Þ r11 ¼ ð1 n0Þrsþ ra ð10Þ r12¼ ra ð11Þ r22¼ n0rfþ ra ð12Þ b ¼mn20=kp ð13Þ where %%
s is the solid stress tensor, %%
t the effective stress tensor of the solid, %%
e the strain tensor of solid,
%%
S the normal stress tensor of the fluid, % d and
%
D are the solid and fluid displacement
vectors, respectively, P2 the perturbed pressure of the fluid inside the porous medium; rs the
solid density; ra the mass coupling effect (neglected in this study), n0 the porosity, m the fluid
viscosity, kp the specific permeability, G and l are the Lame’s constants of elasticity and
%% I the identity matrix.
Combining the continuity equations of the solid and the fluid with the state equation of the fluid, and after linearization of the porosity (see Reference [22]), we find
@P2 @t ¼ K n0 ð1 n0Þr @ % d @t þ n0r @ % D @t ð14Þ for the perturbed pressure. In Equation (14), K is the bulk modulus of compressibility of fluid inside the porous bed.
According to Huang and Chwang [16], the displacement vectors of skeleton and fluid could be separated into dilatational and rotational parts
% d ¼ rjdþ r % Hd ð15Þ % D ¼ rjDþ r % HD ð16Þ
The first terms on the right-hand sides of Equations (15) and (16) represent the dilatational components (two dilatational waves) of the skeleton and fluid displacements, respectively; while the second terms represent rotational components of the skeleton and fluid displacements. The former ones are irrotational while the latter ones are solenoidal. Referring to Morse and Feshbach [23], % Hd and % HDcan be expressed as % Hd ¼ od % eZ and % HD¼ oD % eZ ð17Þ where %
ezis the unit vector normal to the x–y plane.
Substituting Equations (15)–(17) into Equations (4) and (5), we obtain ð2G þ lÞr2jd ð1 n0ÞP2¼ r11 @2j d @t2 þ r12 @2j D @t2 þ b @jd @t @jD @t ð18Þ n0P2 ¼ r12 @2j d @t2 þ r22 @2j D @t2 b @jd @t @jD @t ð19Þ
Gr2od ¼ r11 @2o d @t2 þ r12 @2o D @t2 þ b @od @t @oD @t ð20Þ 0 ¼ r12@ 2o d @t2 þ r22 @2o D @t2 b @od @t @oD @t ð21Þ Substituting Equations (15)–(17) into Equation (14), we get
@P2 @t ¼ K n0 @ @t½ð1 n0Þr 2j dþ n0r2jD ð22Þ
Equations (18)–(22) are the governing equations.
There are three boundaries that require boundary conditions. They are (i) the free surface ½y ¼ h þ zðx; tÞ; (ii) the channel bed ½y ¼ gðx; tÞ; and (iii) the deep far field of the porous bed ½y ! 1:
On the free surface, the kinematic boundary condition is @f1 @y ¼ @z @tþ @z @x U þ @f1 @x ð23Þ and the dynamic boundary condition is
@f1 @t þ 1 2 2U @f1 @x þ @f1 @x 2 þ @f1 @y 2 " # þ gz ¼ 0 ð24Þ
On the porous bed, the continuity of pressure implies that
P1¼ P2 ð25Þ
while the continuity of fluid flux gives @f1 @y @g @x @f1 @x ¼ U @g @xþ @ @t ð1 n0Þ @jd @y @od @x þ n0 @jD @y @oD @x ð1 n0Þ @g @x @jd @y @od @x n0 @g @x @jD @y @oD @x : ð26Þ
The continuity of effective stress of the solid, i.e. n2
%% t ¼ % 0; gives Gð2jd;xyþ od;yy od;xxÞ @g
@x½ð2G þ lÞjd;xxþ 2God;yxþ ljd;yy ¼ 0 ð27Þ
ð2G þ lÞjd;yyþ 2God;xyþ ljd;xx
@g
@xGð2jd;xyþ od;yy od;xxÞ ¼ 0 ð28Þ
At the far field of the porous bed, i.e. y ! 1; the boundary conditions consist of vanishing displacement vectors, i.e.
lim
and no perturbation is propagated to the upstream and downstream boundaries, i.e. lim x!1 @ @nfjd; jD; odg ¼ 0 lim x!1 @f1
@n ikf1¼ 0; periodic wave problem
@f1
@n ¼ 0; transient wave problem
8 > > < > > : ð30Þ
In addition, the kinematic boundary condition at the porous bed interface results in @g @t¼ @ @t @jd @y @od @x @g @x @jd @x þ @od @y ð31Þ which allows to determine the surface elevation of the porous bed, g:
The governing equations (2), (18)–(21) with perturbed pressure in porous bed, Equation (22), surface elevations z and g; Equations (24) and (31), together with other boundary conditions, Equations (23), (25)–(30), and a properly ‘warm-up’ initial condition form a well-posed problem.
2.2. Non-dimensionalization through boundary layer correction
When the poroelastic bed is composed of hard material, the inertial terms of the governing equations (first two terms on the right-hand side of Equations (18)–(21)) are unimportant and they can be neglected (see Reference [24]). The deformation is also negligible. However, when the poroelastic material is soft, the inertial terms are important, and the second dilatational wave travelling through the porous bed has a much shorter wavelength than that of the wave travelling through homogeneous water (see Reference [17]). Thus there exists a small boundary layer near the interface, inside the soft poroelastic bed. Because the effect of the boundary layer will render errors of the partial derivative in the vertical direction for the second longitudinal wave, the traditional estimation of the deformable bed by Stokes expansion based
on e1 will fail [17]. Hsieh et al. [20] proposes a boundary layer correction approach based on
multiple length scales to overcome this difficulty. Since we are interested in bed form deformation in this study, the following discussion will concentrate on soft poroelastic material only.
The coefficients of the governing equations, Equations (18)–(21), range too far to be solved directly by numerical computation so that non-dimensionalization is needed to reduce the difference of the magnitude of each term. We define
x ¼ Lxn ; y ¼ Lyn ð32Þ f1¼ L ffiffiffiffiffiffi gL p fn 1; t1¼ ffiffiffiffiffiffiffiffi L=g p tn 1 ð33Þ jd¼ L2j n d; od¼ L2o n d ð34Þ jD¼ L2jn D; oD¼ L2o n D ð35Þ P2¼ K n0 Pn 2; t2¼ ffiffiffiffiffiffiffiffi L=g p tn 2 ð36Þ
In which, all symbols with an asterisk stand for the dimensionless variables. The order of tn
2 is
found to be Oð1Þ: Equations (2) and (18)–(22) can be rewritten as
#aa1f1;xx*n þ 2#aa2f1;xZ*n þ#aa3f1;ZZ*n þ#aa4f1;x*nþ#aa5f1;Z*n ¼ 0 ð37Þ
ð2G þ lÞ K r 2 *j n d ð1 n0ÞP n 2 ¼ r11gL K @2jn d @t* 2 þ r12gL K @2jn D @t* 2 þ b ffiffiffi g L r L2 K @jn d @tn @jn D @tn ð38Þ n0P n 2 ¼ r12gL K @2jn d @t* 2 þ r22gL K @2jn D @t* 2 b ffiffiffi g L r L2 K @jn d @tn @jn D @tn ð39Þ G Kr 2 *o n d ¼ r11gL K @2on d @t* 2 þ r12gL K @2on D @t* 2 þ b ffiffiffi g L r L2 K @on d @tn @on D @tn ð40Þ 0 ¼r12gL K @2on d @t* 2 þ r22gL K @2on D @t* 2 b ffiffiffi g L r L2 K @on d @tn @on D @tn ð41Þ @Pn 2 @tn ¼ 1 n0 @ @tn½ð1 n0Þr 2 *j n dþ n0r2*j n D ð42Þ
in which, the coefficients #aa1;#aa2;#aa3;#aa4 and #aa5 are referred in Appendix A.
Adding Equation (38) to Equation (39) to eliminate Darcy’s terms, we get ð2G þ lÞ K r 2 *j n d P n 2 ¼ ðr11þ r12ÞgL K @2jn d @t* 2 þ ðr12þ r22ÞgL K @2jn D @t* 2 ð43Þ
and similarly adding Equation (40) to Equation (41), we have
G Kr 2 *o n d¼ ðr11þ r12ÞgL K @2on d @t* 2 þ ðr12þ r22ÞgL K @2on D @t* 2 ð44Þ
The huge coefficient b of Darcy’s terms in Equation (39) is not suitable for numerical computation. However, for soft poroelastic material the magnitude of the second longitudinal wave is rather small and the displacement potentials of solid and fluid are very close, i.e. jn d ffi j n Dand o n d ffi o n
D; by the boundary layer correction proposed by Hsieh et al. [20]. Since
Darcy’s terms are mainly due to the second longitudinal wave, we find that when the porous material is soft, Darcy’s terms are negligible. We therefore can simplify Equations (39) and (41), respectively, to n0 K bL2 ffiffiffi L g s Pn 2 ¼ r12 b ffiffiffi g L r @2jn d @t* 2 þ r22 b ffiffiffi g L r @2jn D @t* 2 ð45Þ 0 ¼r12 b ffiffiffi g L r @2on d @t* 2 þ r22 b ffiffiffi g L r @2on D @t* 2 ð46Þ
Note that Equations (39), (41), (43), (44) and (42) are the full non-dimensional governing equations for the soft porous bed, while Equations (43)–(46) and (42) are the simplified non-dimensional governing equations for the soft porous bed. The non-non-dimensional boundary
conditions are as follows:
1. Kinematic free surface boundary condition @fn 1 @yn¼ 1 g @2fn 1 @t* 2 þ 2 U n þ@f n 1 @xn @2fn 1 @xn @tnþ U n þ@f n 1 @xn 2 @2fn 1 @x* 2 " þ @f n 1 @yn @2fn 1 @yn @tnþ @fn 1 @yn U n þ@f n 1 @xn @2fn 1 @xn @yn : ð47Þ
2. Dynamic free surface boundary condition @fn 1 @tn þ 1 2 2U n@f n 1 @xn þ @fn 1 @xn 2 þ @f n 1 @yn 2 " # þ zn ¼ 0 ð48Þ
3. Continuity of flux at water/porous bed interface @fn 1 @yn @gn @xn @fn 1 @xn ¼ U n@g n @xnþ @ @tn ð1 n0Þ @jn d @yn @on d @xn þ n0 @jn D @yn @on D @xn ð1 n0Þ @gn @xn @jn d @yn @on d @xn n0 @gn @xn @jn D @yn @on D @xn : ð49Þ
4. Continuity of effective stresses (i.e. txy¼ 0 and tyy ¼ 0) at water/porous bed interface
G Kð2j n d;xnynþ o n d;ynyn o n d;xnxnÞ @g n @xn ð2G þ lÞ K j n d;xnxnþ 2G Ko n d;ynxnþ l Kj n d;ynyn ¼ 0 ð50Þ ð2G þ lÞ K j n d;ynyn 2G Ko n d;xnynþ l Kj n d;xnxn @g n @xn G Kð2j n d;xnynþ o n d;ynyn o n d;xnxnÞ ¼ 0 ð51Þ
5. Continuity of pressure at water/porous bed interface
Pn 2K þrfgL @fn 1 @tn þ 1 2 2U n@f n 1 @xnþ @fn 1 @xn 2 þ @f n 1 @yn 2 " # ( ) ¼ 0 ð52Þ
6. Far field of porous bed
lim yn!1ff n 1; j n d; j n D; o n dg ¼ 0 ð53Þ
7. Upstream/downstream boundary lim xn!1 @ @xnfj n d; j n D; o n dg ¼ 0 lim xn!1 @fn 1 @xn ik0f n
1 ¼ 0; periodic wave problem
@fn
1
@xn ¼ 0; transient wave problem
8 > > < > > : ð54Þ
In addition, the non-dimensional kinematic boundary condition at the bed surface is @gn @tn¼ @ @tn @jn d @yn @on d @xn @gn @xn @jn d @xnþ @on d @yn ð55Þ
3. COMPUTATION AND DISCUSSION
To estimate the variation of channel bed accurately and effectively, the proposed numerical computation adopts the boundary-fitted co-ordinate system of Thompson et al. [25] to transform the physical domain into the computational domain by substitution of the Dirichlet boundary conditions. Taking the Laplace equations
xxxþ xyy¼ 0 ð56Þ
Zxxþ Zyy¼ 0 ð57Þ
as the grid generator, the orthogonality of grids on the boundaries can be maintained. Besides the orthogonal boundary-fitted co-ordinate transformations, the new equations are expressed in semi-differential form using the finite difference method. Time and space are discretized separately. The time scale is discretized by the first-order backward difference scheme while the space domain is discretized using the central difference method. The computation is divided into channel flow region and porous bed region. These two regions are carried out separately and sequentially, and they are linked by the boundary conditions at the interface. At the moving boundary, the interpolation adopts a cubic spline to keep the shape of porous bed and re-generate the grids to make sure the orthogonality condition is maintained at each time step. The expressions of the independent variables by the finite difference method are written as:
Fin;j¼5 2F n1 iþ1;j 4 2F n1 iþ2;jþ 1 2F n1 iþ3;j ð58Þ Fin;j¼5 2F n1 i1;j 4 2F n1 i2;jþ 1 2F n1 i3;j ð59Þ
where F is the independent variable for the channel flow or the porous bed, the superscript n refers to the time step and the subscripts i; j are the grid co-ordinates.
In order to verify our formulation, this numerical model without simplification (Equations (37), (39), (41)–(44), along with the boundary conditions equations (47), (49)–(54) and the surface elevation equations (48), (55)) is applied to simulate the problem of a simple periodic wave progression in combination with a constant current. Taking the current
Figure 2. Refined grid near the interface of the computational domain.
U ¼0:3 m=s; water depth h ¼ 3 m; time to warm up t ¼ 60 s; time step Dt ¼ 0:1 s for channel flow, but 0:01 s for porous bed, and the thickness of channel bed ¼ 3 m; the computed results of free surface and bed form are compared with the analytic solution proposed by Hsieh et al. [20]. Figure 2 shows that the meshes are taken as 10 by 10 grids for region 1, but they are much finer close to the interface with region 2. Figure 3 shows good agreement for the free surface and bed form with the solution of Hsieh et al. [20]. Furthermore, the results for the displacement potentials of both the dilatational and rotational waves at selected sections (indicated by the downward arrows in Figure 4) also agree very well with the Hsieh et al. [20] analytical solutions.
Figure 4. Comparison of computed displacement potentials of dilatational wave ðjdÞ and rotational wave ðodÞ with exact solution.
Figure 5. Schematic diagram of a trench at the channel bed.
Figure 6. Transient variation of free surface and a trench at bed form, full formulation ðDt ¼ 1 s; U ¼ 0:3 m=sÞ:
For the first application, the transient response of a small trench at the channel bed, an important problem of river bed dredging, is simulated. Figure 5 is the schematic diagram of a small trench that is analysed. Figure 6 shows the transient transport of the trench by adopting the original complete formulation (Equations (37), (39), (41), (43), (44) with boundary conditions (47)–(54)). Figure 6 shows that the fluctuation of water surface is very large and continuously increases while the bed form remains unchanged. From the results in Figure 6, we find that although small meshes are used near the interface and the time step is small, the numerical model without simplification is still very unstable and will diverge after a very short period of time, say 10 s: This is because that the dilatational waves for water and solid inside the soft porous bed are almost the same, hence Darcy’s terms multiplied by a very large coefficient b (the third terms on the right-hand side of Equations (39) and (41)) are impossible to estimate accurately even after non-dimensionalization. Referring to Hsieh et al. [20], it is noticeable that when the porous bed is soft, the wavelength of the second dilatational wave near the interface is much shorter than that in the homogeneous water, which makes convergence problematic. Thus, the traditional approach of refining meshes near the interface is unsuccessful. However,
Figure 7. Transient variation of free surface and a trench at bed form, simplified formulation ðDt ¼ 60 s; U ¼ 0:3 m=sÞ:
Figure 8. Schematic diagram of channel bed with a downward step.
Figure 9. Transient variation of free surface and a downward-step channel bed, simplified formulation ðDt ¼ 60 s; U ¼ 0:3 m=sÞ:
the magnitude of the second longitudinal wave is rather small and requires only the so-called boundary layer correction (see Reference [20]). Since Darcy’s terms are mainly due to the second longitudinal wave, i.e. the relative velocity of the irrotational part, we find that when the porous material is soft, Darcy’s terms are negligible. After neglecting Darcy’s terms, Equations (39) and (41) turn out to be Equations (45) and (46). This simplified model is stable, larger time step is allowable, and the meshes near the interface no longer need to be refined.
Indeed, this becomes apparent when the same problem as in Figure 6 is analysed using the simplified procedure. (The finite difference scheme is described in Appendix A.) The results are shown in Figure 7. The time step in Figure 7 is 60 s instead of 1 s in Figure 6. After about 12 min; the channel bed begins to deform.
A second application involves the channel bed with a downward step under a uniform flow, as shown in Figure 8. This is an important problem for bridge pier protection applications when river beds are dredged. The transient variations of the free surface and the channel bed are shown in Figure 9, which indicates the deformation of the head-cutting of channel bed together with the upstream going water wave.
5. CONCLUSION
The transient response of poroelastic bed form due to stream flow and non-linear water wave are investigated numerically. The theory of potential flow is applied to channel flow while Biot’s theory of poroelasticity [21] is adopted to deal with the deformable porous bed. The numerical model without simplification can simulate the periodic wave/current problem with a simple or regular geometry only. But it is easy to diverge within very short time for the transient problem with irregular geometry. The transient problem with irregular geometry always needs small grids near the bed/water interface as well as very small time steps. Even then, convergence is not assured. The present work establishes a simplified numerical model based on the boundary layer concept proposed by Hsieh et al. [20] by neglecting Darcy’s terms for a soft poroelastic bed. This makes the study of transient deformations of a channel bed with a trench or a downward step feasible.
ACKNOWLEDGEMENTS
This study is sponsored by the National Science Council of the Republic of China under grants NSC 86-2621-E002-016 and NSC 88-2611-E-002-028.
APPENDIX A: FINITE DIFFERENCE FORMULATIONS Governing equations:
Region 1:
where #aa1¼ x2xnþ x 2 yn ðA2Þ #aa2¼ xxnZxnþ xynZyn ðA3Þ #aa3¼ Z2xnþ Z 2 yn ðA4Þ #aa4¼ xxnxnþ xynyn ðA5Þ #aa5¼ Zxnxnþ Zynyn ðA6Þ Region 2: 2G þ l K ð#aa1f *n
d;xxþ 2#aa2fd*;xZn þ#aa3fd*;ZZn þ#aa4fd*;xnþ#aa5fd*;ZnÞ P2*n
¼ O1 f*n d 2f* n1 d þ f* n2 d Dt*2 þ O2 f*n D 2f* n1 D þ f* n2 D Dt*2 ðA7Þ O3P2*n¼ O4 f*n d 2f* n1 d þ f* n2 d Dt*2 þ O5 f*n D 2fD*n1þ fD*n2 Dt*2 ðA8Þ G Kð#aa1o *n d;xxþ 2#aa2o *n d;xZþ#aa3o *n d;ZZþ#aa4od*;xn þ#aa5o *n d;ZÞ ¼ O1 o*n d 2od*n1þ od*n2 Dt*2 þ O2 o*n D 2oD*n1þ oD*n2 Dt*2 ðA9Þ 0 ¼ O4 o*n d 2od*n1þ od*n2 Dt*2 þ O5 o*n D 2oD*n1þ oD*n2 Dt*2 ðA10Þ P*n 2 2P2*n1þ P2n2 Dt*2 ¼ 1 n0 1 Dtnfð1 n0Þ½ð#aa1f *n
d;xxþ 2#aa2fd*;xZn þ#aa3fd*;ZZn þ#aa4fd*;xn þ#aa5fd*;ZnÞ
ð#aa1fd*;xxn1þ 2#aa2fd*;xZn1þ#aa3fd*;ZZn1þ#aa4fd*;xn1þ#aa5fd*;Zn1Þ
þ n0½ð#aa1fD*;xxn þ 2#aa2fD*n;xZþ#aa3fD*n;ZZþ#aa4fD*;xn þ#aa5fD*;ZnÞ
ð#aa1fD*;xxn1þ 2#aa2fD*n;xZ1þ#aa3fD*n;ZZ1þ#aa4fD*;xn1þ#aa5fD*;Zn1Þg ðA11Þ
where O1¼ ðr11þ r12ÞgL K ðA12Þ O2¼ ðr12þ r22ÞgL K ðA13Þ O3¼ n0 K bL2 ffiffiffi L g s ðA14Þ
O4¼ r12 b ffiffiffi g L r ðA15Þ O5¼ r22 b ffiffiffi g L r ðA16Þ
Boundary conditions at the free surface: (a) Kinematic boundary conditions:
ðf*n 1;xxynþ f*n 1;ZZynÞ ¼1 g f*n 1 2f*n 1 1 þ f*n 2 1 Dt*2 þ 2ðUn þ f*n1 1;x xxnþ f*n1 1;Z ZxnÞ f*n 1;xxxnþ f*n 1;Z f1;x*n1 f *n1 1;Z Zx Dtn ! þ ðf*n 1;xxx 2 xnþ 2f*n 1;xZxxnZxnþ f*n 1;ZZZ 2 xnþ f*n 1;xxxnxnþ f*n 1;ZZxnxnÞ ðUn þ f*n1 1;Z xxnþ f*n1 1;Z ZxnÞ2þ ðf*n1 1;x xynf*n1 1;Z ZynÞ f*n1 1;x xynþ f*n1 1;Z Zyn f*n2 1;x xyn f*n2 1;Z Zyn Dtn ! þ ðf*n1 1;x xynþ f*n1 1;Z ZynÞðU n þ f*n1 1;x xxnþ f*n1 1;Z ZxnÞ ½f*n1 1;xxxxnxynþ f*n1 1;xZðxynZxnþ xxnZynÞ þ f*n1 1;ZZZxnZynþ f1;x*n1xxnynþ f*n1 1;Z Zxnyng ðA17Þ
(b) Dynamic boundary condition: f*n 1 f* n1 1 Dtn þ 1 2½2U n ðf*n 1;xxxnþ f*n 1;ZZxnÞ þ ðf*n 1;xxxnþ f*n 1;ZZxnÞ2 þ ðf*n 1;xxynþ f*n 1;ZZynÞ 2 þ zn¼ 0 ðA18Þ
Boundary conditions at the porous bed interface: (a) Continuity of flux:
ðf*n 1;xxynþ f*n 1;ZZynÞ y*n1 x x*n1 x ðf*n 1;xxxnþ f*n 1;ZZxnÞ ¼ Uny *n1 x x*n1 x þ 1 Dtn ð1 n0Þ 1 y*n1 x x*n1 x ! " ðj*n d;xxynþ jd*n ;ZZyn od*;xnxxn od*n ;ZZxnÞ (
þ n0 1 y*n1 x x*n1 x ! ðj*n D;xxynþ jD*n ;ZZyn oD*n ;xxxn oD*n ;ZZxnÞ : ð1 n0Þ 1 y*n2 x x*n2 x ! ðj*n1 d;x xynþ j*n1 d;Z Zyn od*n1 ;x xxn o*n1 d;Z ZxnÞ " þ n0 1 y*n2 x x*n2 x ! ðj*n1 D;x xynþ jD*n1 ;Z Zyn oD*n1 ;x xxn oD*n1 ;Z ZxnÞ :g ðA19Þ
(b) Continuity of effective stress:
G Kf2½j *n d;xxxxnxynþ j*n d;xZðxxnZynþ xynZxnÞ þ j*n d;ZZZxnZynþ jd*;xnxxnynþ jd*n ;ZZxnyn þ ½o*n d;xxx 2 ynþ 2od*;xZn xynZynþ od*n ;ZZZ 2 ynþ od*;xnxynynþ od*n ;ZZynyn þ ½o*n d;xxx 2 xnþ 2od*;xZn xxnZxnþ od*n ;ZZZ 2 xnþ od*;xnxxnxnþ od*n ;ZZxnxng y *n1 x x*n1 x 2G þ l K ½j *n d;xxx 2 xnþ 2jd*n ;xZxxnZxnþ jd*;ZZn Z2xnþ jd*n ;xxxnxnþ j*n d;ZZxnxn þ 2G K½o *n d;xxxxnxx yn þ o *n d;xZðxxnZynþ xynZxnÞ þ o*n d;ZZZxnZynþ od*n ;xxxnynþ o*n d;ZZxnyn þ l K½j *n d;xxx 2 ynþ 2jd*;xZn xynZynþ jd*;ZZn Z2ynþ jd*;xnxynynþ j*n d;ZZynyn : ¼ 0 ðA20Þ 2G þ l K ½j *n d;xxx 2 ynþ 2jd*;xZn xynZynþ jd*n ;ZZZ 2 ynþ jd*;xnxynynþ jd*n ;ZZynyn 2G K½o *n d;xxxxnxynþ o*n d;xZðxxnZynþ xynZxnÞ þ o*n d;ZZZxnZynþ od*n ;xxxnynþ o*n d;ZZxnyn þ l K½j *n d;xxx 2 xnþ 2jd*n ;xZxxnZxnþ jd*;ZZn Z2xnþ jd*n ;xxxnxnþ jd*n ;ZZxnxn : y *n1 x x*n1 x G Kf2½j *n d;xxxxnxynþ jd*n ;xZðxxnZynþ xynZxnÞ þ j*n d;ZZZxnZynþ j*n d;xxxnynþ jd*n ;ZZxnyn þ ½o*n d;xxx 2 ynþ 2od*;xZn xynZynþ od*n ;ZZZ 2 ynþ od*;xnxynynþ o*n d;ZZynyn þ ½o*n d;xxx 2 xnþ 2od*;xZn xxnZxnþ od*n ;ZZZ 2 xnþ od*;xnxxnxnþ od*n ;ZZxnxng ¼ 0 ðA21Þ
(c) Continuity of pressure: P*n 2 K þrfgL f*n 1 f *n1 1 Dtn þ 1 2½2U n ðf*n 1;xxxnþ f*n 1;ZZxnÞ þ ðf*n 1;xxxnþ f*n 1;ZZxnÞ2þ ðf*n 1;xxynþ f*n 1;ZZynÞ2 : ¼ 0 ðA22Þ
(d) Far field of porous bed: fn 1¼ j n d ¼ j n D¼ o n d¼ o n D¼ 0 ðA23Þ
(e) Upstream boundary:
f*n 1i;j¼ 5 2f *n1 1iþ1;j 4 2f *n1 1iþ2;jþ 1 2f *n1 1iþ3;j ðA24Þ 5j*n diþ1;j 4j* n diþ2;jþ j* n diþ3;j¼ 0 ðA25Þ 5j*n
Diþ1;j 4jDiþ*n 2;jþ jDiþ*n 3;j¼ 0 ðA26Þ
5o*n diþ1;j 4o* n diþ2;jþ o* n diþ3;j¼ 0 ðA27Þ (f) Downstream boundary: f*n 1i;j¼ 5 2f *n1 1i1;j 4 2f *n1 1i2;jþ 1 2f *n1 1i3;j ðA28Þ 5j*n
di1;j 4jdi*n2;jþ jdi*n3;j¼ 0 ðA29Þ
5j*n Di1;j 4j* n Di2;jþ j* n Di3;j¼ 0 ðA30Þ 5o*n
di1;j 4odi*n2;jþ odi*n3;j¼ 0 ðA31Þ
After the solutions are obtained, the variation of bed form can be found as g*n ¼ g*n1þ ðj*n d;xxynþ jd*n ;ZZyn o*n d;xxxn od*n ;ZZxnÞ h y *n1 x x*n1 x ðj*n d;xxxnþ jd*n ;ZZxnþ od*;xnxynþ od*n ;ZZynÞ ðj*n1 d;x xynþ j*n1 d;Z Zyn od*n1 ;x xxn o*n1 d;Z ZxnÞ h y *n2 x x*n2 x ðj*n1 d;x xxnþ j*n1 d;Z Zxnþ od*n1 ;x xynþ o*n1 d;Z ZynÞ : ðA32Þ REFERENCES
1. Graf WH. Hydraulics of Sediment Transport. Table 11.1 and Figure 11.2. McGraw-Hill: New York, 1971; 277. 2. Shibayama T. Sediment transport mechanism and two-dimensional beach transformation due to waves. D. Eng.
Dissertation, University of Tokyo, Tokyo, 1984.
3. Horikawa K (ed.). Nearshore Dynamics and Coastal Processes. Figure 3.2. University of Tokyo Press: Tokyo, 1988; 172.
5. Sleath JFA. A numerical study of the influence of bottom roughness on mass transport by water waves. Proceedings of the International Conference on Numerical Methods Fluid Dynamics, Southampton, 1973.
6. Longuet-Higgins MS. Oscillating flow over steep sand ripples. Journal of Fluid Mechanics 1981; 107:1–35. 7. Vittori G, Blondeaux P. Sand ripples under sea waves}Part 3. Brick-pattern ripple formation. Journal of Fluid
Mechanics 1992; 239:23–45.
8. Darwin GH. On the formation of ripple-marks. Proceedings of the Royal Society, London, 1884; 36:18–43. 9. Bagnold RA. The movement of the desert sand. Proceedings of the Royal Society, London Series A 1936;
892(157):594–620.
10. Exner FM. .UUber die wechselwirkung zwischen wasser und geschiebe in flussen. Sitzungsberichte der Akademie der Wissenschaften. Wein, Heft 3–4, 1925.
11. Anderson AG. The characteristics of sediment waves formed by flow in open channels. Proceedings of the Third Midwestern Conference in Fluid Mechanics, Minneapolis, March 1953; 379–395.
12. Kennedy JF. The mechanics of dunes and antidunes in erodible-bed channels. Journal of Fluid Mechanics 1963; 16:521–544.
13. Biot MA. Theory of propagation elastic waves in a fluid saturated porous solid. I. Low-frequency range. Journal of American Statistical Association 1956; 28:168–178.
14. Mei CC, Foda MA. Wave-induced responses in a fluid filled poroelastic solid with a free surface-a boundary layer theory. Geophysical Journal of the Royal Astronomical Society 1981; 66:597–631.
15. Huang LH, Song CH. Dynamic response of poroelastic bed to water waves. Journal of Hydraulic Engineering, ASCE 1993; 119:1003–1020.
16. Huang LH, Chwang AT. Trapping and absorption of sound waves. II: a sphere covered with a porous layer. Wave Motion 1990; 12:401–414.
17. Chen TW, Huang LH, Song CH. Dynamic response of poroelastic bed to nonlinear water waves. Journal of Engineering Mechanics, ASCE 1997; 123(10):1041–1049.
18. Huang LH, Chiang YL. Bed form formation in alluvial channels. Proceedings of the Seventh International Symposium of River Sedimentation, Hong Kong, 1998.
19. Huang LH, Hsieh PC, Wang TW. Linear oscillatory water wave propagating over soft poroelastic bed. Biot Conference on Poromechanics, Louvain-la-Neuve, Belgium, 1998.
20. Hsieh PC, Huang LH, Wang TW. Dynamic response of soft poroelastic bed to nonlinear waver wave}boundary layer correction approach. Journal of Engineering Mechanics, ASCE 2000; 126(10):1064–1073.
21. Biot MA. Mechanics of deformation and acoustic propagation in porous media. Journal of Applied Physics 1962; 33(4):1482–1498.
22. Verruijt A. In Elastic Storage of Aquifers}Flow through Porous Media, De DeWiest RJM (ed.). Academic Press: New York, 1969.
23. Morse PM, Feshbach H. Method of Theoretical Physics. McGraw-Hill: New York, 1978.
24. Madson OS. Wave-induced pore pressure and effective stress in a porous bed. Geotechnology 1978; 28:377–393. 25. Thompson JF, Thames FC, Mastin CW. Automatical numerical generation of body-fitted curvilinear coordinate
systems for fields containing any number of arbitrary two-dimensional bodies. Journal of Computational Physics 1974; 15:299–317.
