Xem mẫu
- T0 g H e − mx
2π 2π
ξ= − n cos (t − nx / C 0 ) + r sin (t − nx / C 0 ) +
2π D s T0 T0
1 πH n(1 + 2s ) − 2 mx 4π
(t − nx / C 0 ) −
e cos
2 βDT0 s T0
CHƯƠNG 2 – NHỮNG PHƯƠNG PHÁP SỐ TRỊ
4π
q (1 + 2l ) − px
(t − qx / C 0 ) +
e cos
TÍNH THỦY TRIỀU
l T0
1 H s − 0,5 −2 mx 4π
(t − nx / C 0 ) −
e sin
Trong chương 1 đã xét những khái niệm cơ bản về hiện tượng thủy
4 D sn T0
triều trong đại dương và những lý thuyết giải thích sự hình thành những
4π
0,5(l − 0,5) − px đặc điểm cơ bản, chung nhất của hiện tượng triều xảy ra trong biển thực.
(t − qx / C 0 ) .
e sin (1.134)
Tuy nhiên như đã nhận xét, những thuyết này chưa thể cung cấp những
lq T0
công thức, những phương pháp để tính toán những đặc trưng thủy triều
với độ chính xác cần thiết trong thực tiễn. Điều này chủ yếu do ở biển và
đại dương thực các sóng thủy triều lan truyền trong những điều kiện vật
lý, điều kiện hình học đường bờ và địa hình đáy biển phức tạp hơn nhiều
so với những sơ đồ đã xét bằng giải tích. Do đó, trong chương này, chúng
ta sẽ xem xét những phương pháp thủy động số trị để giải các phương
trình chuyển động thủy triều nhằm tính tới được những điều kiện gần
thực ở biển.
2.1. PHƯƠNG PHÁP DEFANT
Xét chuyển động thủy triều trong kênh nửa kín. Giả sử kênh rất hẹp,
có thể bỏ qua ảnh hưởng của lực Coriolis. Ma sát ở đáy và thành kênh
không có. Chuyển động ngang của các hạt nước không đổi trong mặt
phẳng vuông góc với phương truyền thủy triều, tức trong thiết diện ngang
kênh. Tốc độ thành phần ngang u có thể là một hàm số chỉ theo hướng
x và thời gian t [3].
39
- Bây giờ chúng ta sẽ nhận các phương trình thuận tiện cho việc tích
phân bằng số. Đặt gốc toạ độ lên mặt phẳng đáy, trục x hướng dọc theo trong đó ζ và ξ tuần tự là các biên độ của dao động mực nước và
kênh, trục z thẳng đứng hướng lên trên.
quãng đường dịch chuyển ngang của hạt nước trong chuyển động triều.
Phương trình chuyển động theo hướng trục x (1.19) và phương
Ký hiệu diện tích mặt cắt ngang kênh là S , chiều rộng kênh là b và
trình liên tục (1.28) sẽ có dạng đơn giản sau đây: D = S /b .
∂ζ
∂u
= −g Khi đó các phương trình (2.4) và (2.5) sẽ dẫn đến dạng các phương
; (2.1)
∂t ∂x trình cho biên độ các dao động [6]:
∂ζ ∂u 2
= −D . dζ 1 2π
(2.2)
ξ;
=
∂t ∂x (2.6)
dx g T
Nếu sử dụng đại lượng di chuyển ngang ξ của hạt nước liên hệ với
d [ S ( x)ξ )
tốc độ u theo định nghĩa = ζ b( x ) . (2.7)
dx
t
ξ = udt , (2.3) Dùng điều kiện triệt tiêu chuyển động ngang ở đầu kín của kênh
( x = 0 ) làm điều kiện biên theo x :
0
thì phương trình chuyển động (2.1) được viết lại thành ξ =0 (2.8)
x =0
∂ξ ∂ζ
2
= −g , (2.4) và cho trước dao động thẳng đứng của mực nước ở cửa mở của kênh
2
∂t ∂x ( x = ):
và phương trình liên tục (2.2) thành
ζ =ζ . (2.9)
x =
∂ζ ∂ 2ξ
= −D . (2.5) Như vậy hệ phương trình (2.6), (2.7) và các điều kiện biên (2.8) và
∂t ∂x∂t
(2.9) hoàn toàn xác định trường dao động triều trong kênh.
Giả sử dao động thủy triều của mực nước và di chuyển ngang là các
Bây giờ ta chia kênh ra làm nhiều đoạn bằng một loạt các thiết diện
hàm điều hoà thời gian dạng:
thẳng đứng vuông góc với trục dọc kênh (hình 2.1). Khoảng cách giữa
2π
ζ = ζ cos hai thiết diện liền nhau bằng Δx . Ký hiệu Δζ là số gia biên độ mực
t
T
nước qua khoảng Δx . Từ phương trình (2.6) sẽ nhận được
2π
ξ = ξ cos t 4π 2
T Δζ = ξ Δx . (2.10)
gT 2
40
- triều, tính theo công thức
ζ j + ζ j −1
q j = q j −1 + Rj ,
2
với q j = 0 ở đầu kín của kênh j = 0 theo điều kiện (2.8).
Sternec và Defant khi mới xây dựng phương pháp này, năm 1915-
1919, đã dùng nó để tính thủy triều cho Đại Tây Dương, biển Ađriatic,
Địa Trung Hải và nhiều biển khác. Kết quả tương đối thoả mãn khi tính
Hình 2.1. Sơ đồ kênh trong phương pháp tích phân từng bước dao động trung bình theo thiết diện ngang của kênh. Tuy nhiên phương
pháp vừa trình bày không tính đến lực Coriolis, nên không thể áp dụng
Dịch chuyển ngang ξ được tìm nhờ phương trình (2.7). Tích phân
đối với những biển không có dạng kênh hẹp. Ngày nay sơ đồ tính toán
phương trình này theo từ 0 đến x và dùng điều kiện biên (2.8) ta được
trên với những cải tiến nhất định có thể sử dụng để tính sự truyền triều
x
b trong các vùng cửa sông, các sông. Về sau này Hansen (năm 1949, 1952)
ξ = − ζ dx . (2.11)
S và sau nữa là Polukarov (năm 1956, 1957, 1960) [10] đã đưa ra những
0
mô hình số trị đầy đủ hơn, tránh được những thiếu sót của phương pháp
Bây giờ ta tích phân hệ phương trình (2.10), (2.11) được thực hiện
Defant. Chúng ta sẽ tìm hiểu các phương pháp này qua việc xem xét mô
bằng phương pháp số “từng bước về phía trước”. Đối với trường hợp
hình số trị của Hansen ở mục tiếp theo.
sóng triều là sóng đứng, các công thức (2.10), (2.11) chuyển thành dạng:
ζ j + ζ j −1 2.2. PHƯƠNG PHÁP HANSEN
ζ j = ζ j −1 + a ; (2.12)
2
ξ j −1
2.2.1. Các phương trình và điều kiện biên
1
ξj = − q j −1 + ζ j −1 + a , (2.13)
aR 4
2
Hansen đã xuất phát từ hệ phương trình chuyển động sóng dài có kể
S j 1 − j
đến ma sát rối thẳng đứng, trong đó các ứng suất ma sát rối tại đáy được
4S j
xấp xỉ bằng quy luật tuyến tính (xem [6]). Trong trường hợp này hệ các
trong đó phương trình chuyển động và phương trình liên tục có dạng (xem các
4π 2 phương trình (1.31), (1.32) và (1.28))
Δx ,
a=
∂ζ
∂u
gT 2
− fv = − g − ru ; (2.14)
∂t ∂x
R j − diện tích mặt kênh giữa hai thiết diện; q − lưu lượng của dòng
41
- ∂ζ
∂v Trước hết nhân hai phương trình chuyển động với D . Sau đó lấy
+ fu = − g − rv ; (2.15)
đạo hàm phương trình (2.18) theo x , lấy đạo hàm phương trình (2.19)
∂t ∂y
theo y rồi cộng hai phương trình lại (thực hiện toán tử phân kỳ ngang),
∂ζ ∂uD ∂vD
+ + = 0. (2.16) nhận được:
∂t ∂x ∂y
∂u D ∂v D ∂v D ∂u D ∂D ∂ζ ∂D ∂ζ
δ − f ∂x + ∂y = − gD∇ ζ − g ∂x ∂x + ∂y ∂y .
2
+
Khi hệ số ma sát được cho trước thì các phương trình (2.14)−(2.16)
∂x ∂y
liên hệ ba hàm số cần tìm: các thành phần tốc độ u , v và độ cao ζ của
(2.21)
mặt biển so với mực trung bình.
Lấy đạo hàm phương trình (2.18) theo y , phương trình (2.19) theo
Cũng như trong mục trước, các đại lượng u, v và ζ biến thiên với
x rồi lấy phương trình thứ hai trừ đi phương trình thứ nhất (thực hiện
thời gian theo quy luật điều hoà đơn giản, viết dưới dạng phức như sau
toán tử xoáy), nhận được
u u
−i σ t ∂v D ∂u D ∂u D ∂v D ∂D ∂ζ ∂D ∂ζ
v = v ⋅ e , δ
(2.17) − + + = − g ∂x ∂y − ∂y ∂x .
f
∂y
∂x ∂x ∂y
ζ ζ
(2.22)
trong đó σ − tốc độ góc của dao động triều; u , v , ζ − những biên độ
Trong các biểu thức nhận được u D và v D là những thành phần
phức của các hàm tương ứng.
dòng toàn phần của triều lưu. Bây giờ nếu loại xoáy vận chuyển toàn
Thế (2.17) vào hệ các phương trình (2.14)−(2.16) và giản ước thừa phần ra khỏi hai phương trình vừa nhận được (bằng cách nhân phương
− iσ t
số e trình thứ nhất với δ , phương trình thứ hai với f rồi cộng hai phương
ta được hệ phương trình viết cho các biên độ
∂ζ trình lại), ta có
δu − fv = − g ; (2.18)
∂x ∂u D ∂v D
(δ 2 + f 2 ) = − gδD∇ ζ − gδ I ( D, ζ ) − gf J ( D, ζ )
2
+
∂ζ ∂x ∂y
δv + fu = − g ; (2.19)
∂y
(2.23)
∂u D ∂v D với các ký hiệu
− iσζ = 0 .
+ (2.20)
∂x ∂y
∂D ∂ζ ∂D ∂ζ
I ( D, ζ ) = +
ở đây δ = r − iσ .
∂x ∂x ∂y ∂y
Bây giờ ta biến đổi các phương trình này để nhận được một phương
trình cho một ẩn là hàm ζ .
42
- ∂D ∂ζ ∂D ∂ζ còn ở biên lỏng G2 biết trước giá trị mực nước
J ( D, ζ ) = −
∂x ∂y ∂y ∂x ζ ( x, y ) G = ϕ ( x , y ) , (2.26)
2
Dùng phương trình (2.20) để loại biểu thức phân kỳ dòng toàn phần ở đây α và β − các góc giữa pháp tuyến trong của bờ với các trục x và
ra khỏi phương trình (2.23), giả thiết (δ 2 + f 2 ) khác không, cuối cùng y (hình 2.2). Bài toán này gọi là bài toán biên hỗn hợp.
ta nhận được phương trình vi phân mô tả dao động mặt biển
Tính đơn trị của nghiệm cũng tồn tại cả trong trường hợp khi các giá
I ( D, ζ ) iσ
f trị của hàm ζ được biết trước trên khắp vòng biên vùng biển nghiên
∇ζ + J ( D, ζ ) + (δ 2 + f 2 ) ζ = 0 .
2
+ (2.24)
δD gDδ
D cứu:
Phương trình (2.24) là phương trình vi phân đạo hàm riêng dạng ζ ( x, y ) G = ψ ( x, y ) , (2.27)
elliptic với các hệ số phức của hàm phức ζ .
(tức bài toán toán biên loại một) [6].
Sự khác nhau giữa bài toán biên loại một và bài toán biên hỗn hợp là
ở chỗ trong bài toán biên loại một các giá trị của hàm mực nước được cho
trước trên toàn đường biên, khi giải phương trình (2.24) cho hàm mực
nước ta chỉ cần tính giá trị của hàm này cho những điểm bên trong của
Hình 2.2. Biên cứng
miền tính. Với bài toán biên loại hỗn hợp cần ít thông tin đầu vào hơn vì
điều kiện biên (2.25) thực chất là điều kiện lý thuyết thuần tuý, không
yêu cầu dữ liệu thực. Song với bài toán này khi giải phương trình (2.24)
ta cần tính hàm mực nước cho cả các điểm bên trong miền tính và các
điểm trên biên cứng và do đó về phương diện kỹ thuật giải số bài toán
Hansen (1952) đã chứng minh rằng đối với trường hợp vùng nghiên
này sẽ khó khăn hơn.
cứu có hệ số ma sát không bằng không, nghiệm của phương trình (2.24)
Nhiệm vụ tiếp theo là tìm các biểu thức tính biên độ tốc độ dòng
khi cho trước điều kiện biên hỗn hợp sẽ xác định đơn trị. Vì vậy nếu vùng
triều theo mực nước. Muốn vậy sử dụng các phương trình (2.18) và
biển giới hạn bởi đường biên kín G , một phần G1 của nó là đường bờ,
(2.19). Nhân phương trình (2.18) với δ , nhân phương trình (2.19) với f
phần còn lại G2 là biên lỏng, thì hàm ζ được xác định đơn trị trong
rồi cộng hai kết quả với nhau ta sẽ được biểu thức của u và trừ hai kết
khắp vùng biển khi ở biên cứng G1 cho trước điều kiện không chảy
quả cho nhau ta sẽ được biểu thức của v :
xuyên qua biên
(u cos α + v cos β ) G1 = 0 , (2.25)
43
- ∂ζ ∂ζ ∂D ∂ζ ∂D ∂ζ
g
I ( D, ζ ) =
δ + ≅
u =− +f
2
δ + f ∂x ∂x ∂y ∂y ∂y
2
∂y
(2.28)
[ ]
1
∂ζ ∂ζ
g ( Dl +1, k − Dl −1, k )(ζ l +1, k − ζ l −1, k ) + ( Dl , k +1 − Dl , k −1 )(ζ l , k +1 − ζ l , k −1 )
f ∂x + δ ∂y
v= 2
4h 2
δ + f2
1
I P ( D, ζ )
≡
Nếu bên trong vùng nghiên cứu và trên các biên của nó đã tính được
4h 2
hoặc cho trước các giá trị hàm ζ , thì theo các biểu thức (2.28) dễ dàng
∂D ∂ζ ∂D ∂ζ
J ( D, ζ ) = − ≅
tính được u và v . ∂x ∂y ∂y ∂x
[ ]
1
( Dl +1, k − Dl −1, k )(ζ l , k +1 − ζ l , k −1 ) − ( Dl , k +1 − Dl , k −1 )(ζ l +1, k − ζ l −1, k )
2.2.2. Sơ đồ sai phân hữu hạn giải các phương trình
4h 2
Vùng biển được chia bằng mạng lưới đều (hình 2.3). Đối với bài 1
J P ( D, ζ )
≡
toán loại một, theo các điều kiện biên (2.27) ta xác định các giá trị hàm 4h 2
ζ ở dãy nút ngoài của vùng lưới G ' :
iσh 2 2
(δ + f 2 ) ; h − bước lưới;
μ − thông số không thứ nguyên, bằng
ζ ( x, y ) G ' = ψ ′( x, y ) . (2.29) gDδ
các chỉ số l , k xác định vị trí của từng nút bên trong vùng lưới.
Ở các nút trong của lưới phương trình vi phân (2.24) được được thay
bằng tương tự sai phân hữu hạn của nó Nếu l biến thiên từ 0 đến N , và k từ 0 đến M , thì lưới sẽ chứa
( N − 1)(M − 1) nút trong. Giá trị của hàm ζ ở mỗi nút trong là những
1
f
) + μζ = 0 ,
∇2 ζ + I P ( D , ζ ) + δ J P ( D, ζ (2.30)
P
4D giá trị cần tìm. Vậy nếu viết phương trình (2.30) cho từng điểm trong
l , k của lưới thì ta sẽ có một hệ phương trình đại số gồm
trong đó ∇ 2 ζ , I P ( D, ζ ), J P ( D, ζ ) − tuần tự là các tương tự sai phân
P
( N − 1)(M − 1) phương trình chứa đúng ( N − 1)(M − 1) ẩn số. Như vậy
của các toán tử Laplacian, I ( D, ζ ) và J ( D, ζ ) nhận được bằng phép
giải bài toán biên loại một dẫn đến giải hệ ( N − 1)( M − 1) phương trình
xấp xỉ sai phân hữu hạn trung tâm:
đại số tuyến tính.
∇ζ
2
1
∇ 2ζ ≅ (ζ l +1, k + ζ l , k +1 + ζ l −1, k + ζ l , k −1 − 4ζ l , k ) ≡
P
2 2
h h
44
- Cách đơn giản nhất để xấp xỉ sai phân các phương trình này cho
những điểm biên là thay các đạo hàm bằng sai phân hữu hạn một chiều.
Kết hợp những phương trình sai phân vừa nhận được cho các điểm
nút biên với các phương trình sai phân cho những nút bên trong lưới ta sẽ
được một hệ phương trình đại số tuyến tính trong đó số phương trình
bằng số điểm nút bên trong cộng với số nút ở biên cứng.
Hệ phương trình đại số tuyến tính nhận được sẽ có nghiệm đơn trị
chỉ trong trường hợp định thức các hệ số của hệ khác không [6]. Nếu điều
kiện này thoả mãn thì có thể giải hệ phương trình đại số tuyến tính bằng
phương pháp thế hoặc phương pháp ma trận. Cũng có thể giải hệ đó bằng
phương pháp lặp, nhưng mỗi lần giải phải kiểm tra tính hội tụ của
phương pháp. Hansen cho biết rằng khi định thức có trị số nhỏ tính hội tụ
của nghiệm bài toán bị phá vỡ.
2.2.3. Nhận xét về phương pháp Hansen qua thực tế tính thủy
triều
Hình 2.3. Sơ đồ vùng tính và lưới sai phân trong phương pháp Hanxen
Khi giải bài toán biên loại hỗn hợp hàm ζ ở từng nút trong của Những công trình tính thủy triều ở Đại Tây Dương (Hansen, 1949;
Boris, 1961) và Thái Bình Dương (Bogđanov, Kim, Magaric, 1964) và ở
vùng lưới cũng cần phải thoả mãn phương trình sai phân (2.30). Tuy
các biển khác như Bắc Hải (Hansen, 1952), Hoàng Hải (Boris, 1958),
nhiên, khác với trường hợp đã xét ở trên, các giá trị của hàm ở các nút
biển Nauy và biển Grinlen (Nhecrasov, 1962, 1965)... xác nhận rằng
trên vòng biên cứng bây giờ lại phải xác định dựa theo điều kiện biên
phương pháp Hansen không những cho bức tranh chung, mà cả những nét
(2.25).
chi tiết trong sự phân bố thủy triều trên các biển này. Các nhà khoa học
Kết hợp các phương trình (2.28) và điều kiện biên (2.25) có thể nhận
Việt Nam như Nguyễn Ngọc Thụy (1969) [18], Đặng Công Minh (1975)
được các phương trình tính ζ cho những điểm trên biên cứng như sau:
[14] cũng đã sử dụng phương pháp Hansen để nghiên cứu đặc điểm
∂ζ ∂ζ truyền sóng thủy triều ở biển Đông.
δ +f =0 cho biên kinh tuyến
∂x ∂y Tuy nhiên phương pháp này có những thiếu sót sau:
∂ζ ∂ζ a) Không thể tính thủy triều cho những biển sâu nằm gần vùng vĩ
+δ =0
f cho biên vĩ tuyến.
∂x ∂y tuyến "tới hạn", nơi tốc độ góc của phân triều cần tính xấp xỉ bằng thông
45
- số Coriolis; b) Cách đánh giá ứng suất ma sát đáy trong mô hình rất thô. ma sát trên mặt tự do (điều kiện (1.25)) và cho áp suất khí quyển trên mặt
tự do P0 = const . Bây giờ nếu tính tới hiệu ứng ma sát do gió tác động
Hansen khi tính toán đã cho hệ số ma sát tỉ lệ với độ sâu biển và tốc độ
triều lưu. Nhưng bản thân tốc độ dòng triều là đại lượng chưa biết cần tìm lên mặt nước
trong khi giải bài toán và trong thực tế hệ số ma sát phải xem như đã
Tx Ty
∂u ∂v
=− =−
k và k (2.31)
được biết trước (theo kết quả đo triều lưu cực đại).
ρ ρ
∂z ∂z
Những nghiên cứu lý thuyết và thử nghiệm (Kagan, 1968) [6] chỉ ra
và khi tích phân phương trình thủy tĩnh chú ý tới sự biến đổi của áp suất
rằng ma sát rối thẳng đứng chỉ góp phần ảnh hưởng tới sự phân bố thẳng
khí quyển theo các phương ngang (xem phương trình (1.18)), thì hệ
đứng theo độ sâu của tốc độ dòng triều ở lớp biên gần đáy biển. Trong
phương trình chuyển động sóng dài sẽ được bổ sung bằng các số hạng
toàn bề dày còn lại của biển với độ sâu lớn có thể bỏ qua lực ma sát rối.
chứa ứng suất gió và građien khí áp như sau:
Điều này làm cho phương pháp Hansen không áp dụng được cho các
1 ∂P a Tx
∂u ∂u ∂u r
vùng vĩ độ "tới hạn". Một trong những cách khắc phục nhược điểm này là
u 2 + v2 u
+u + v − fv = − + −
ρ ∂x ρ ( D + ζ ) D + ζ
đề xuất của Nhecrasov và Kagan (1965, 1966) đưa thành phần ma sát rối ∂t ∂x ∂y
ngang vào các phương trình chuyển động [6].
1 ∂P a Ty
∂v ∂v ∂v r
u2 + v2 v
+ u + v + fu = − + −
Trong các mô hình tính thủy triều hiện đại người ta có thể tính tới cả ρ ∂y ρ ( D + ζ ) D + ζ
∂t ∂x ∂y
những số hạng phi tuyến trong các phương trình chuyển động, sử dụng
∂ζ ∂ ( D + ζ )u ∂ ( D + ζ )v
những phương trình đầy đủ dưới dạng (1.31)−(1.32), tính toán thủy triều =− − (2.32)
∂t ∂x ∂y
có kể tới sự tương tác của nó với những dao động mực nước tổng cộng,
ngoài dao động thủy triều có thể tính tới những dao động nguồn gốc do Trong các phương trình trên bây giờ ta dùng ký hiệu T x , T y − ứng
gió, nước dâng, ảnh hưởng của các dòng nước lục địa... suất gió lên mặt nước tuần tự theo các trục x và y , P a − áp suất khí
quyển trên mặt biển. Khi cho trước điều kiện biên ở cửa biển là dao động
2.3. MÔ HÌNH DAO ĐỘNG MỰC NƯỚC TỔNG CỘNG TRONG
thủy triều, thì hệ này sẽ mô tả sự lan truyền thủy triều từ đại dương vào
BIỂN VEN
thủy vực đang xét dưới ảnh hưởng của trường gió và trường khí áp, tức
Trong mô hình này chuyển động của nước trong thủy vực cũng tuân có thể khảo sát được hiệu ứng tổng cộng của thủy triều và các quá trình
theo hệ phương trình chuyển động sóng dài trong nước nông và phương khí quyển .
trình cân bằng thể tích nước (gọi là hệ phương trình sóng dài trong nước Khi đó điều kiện tại biên lỏng (phía biển) là cho trước dao động thực
nông) nhưng có tính tới khá đầy đủ các lực gây dao động mực nước. Như tổng cộng của mực nước
đã thấy, khi xây dựng các phương trình chuyển động thủy triều (1.31) và ζ = ζ ( x, y , t ) , (2.33)
(1.32) ở mục 1.5 chương 1, chúng ta đã cho điều kiện triệt tiêu ứng suất
hoặc cho biến thiên mực nước bằng phương trình độ cao mực nước thủy
46
- triều (xem chương 3) nếu chỉ khảo sát dao động thủy triều: Trong thực hành tích phân số hệ phương trình trên máy tính có nhiều
cách khác nhau để hiện thực các thủ tục sai phân hoá các phương trình và
n
ζ t = f i H i cos[qi t + (V0 + u ) i − g i ] . (2.34) điều kiện biên vừa nhận xét. Dưới đây là thí dụ các công thức sai phân
i =1
tổng quát đơn giản viết cho trường hợp bỏ qua các số hạng phi tuyến
Các điều kiện tại biên cứng (bờ biển) vẫn tương tự như trong trường không gian trong các phương trình chuyển động của (2.32):
hợp bài toán Hansen.
Δt ~ Δt ~
~ ~
ζ i', j = ζ i , j − ( Di , j ui , j − Di , j −1ui , j −1 ) − ( Li , j vi , j − Li−1, j vi−1, j ) ;
Khi muốn tính tới hiệu ứng của dòng nước sông thì tại các điểm biên
Δx Δx
gắn với cửa sông cho trước lưu lượng sông hoặc tốc độ dòng chảy sông. x a a
Δt Ti , j Pi , j +1 − Pi , j
gΔt '
~
Tại thời điểm ban đầu t = 0 cho các trường mực nước và vận tốc (ζ i , j +1 − ζ i', j ) +
ui , j + fΔtK i , j − −
ρ Di , j ρΔx
Δx
bằng không. ui' , j = ;
rΔt ~
1 + ~ (ui2, j + K i2, j )1 / 2
Giải hệ phương trình với các điều kiện biên sẽ tìm được dòng chảy
Di , j
và độ cao mực nước tổng cộng tại mỗi điểm của vùng biển theo thời gian.
y a a
Δt Ti , j Pi +1, j − Pi , j
Cần nhận xét rằng hệ phương trình (2.32) ngoài những bổ sung đã gΔt '
~'
(ζ i +1, j − ζ i , j ) +
'
+ fΔtSi , j − ~−
vi , j
nêu trên đây, nó còn tính tới hiệu ứng phi tuyến khá đầy đủ nhờ các số ρ Li , j ρΔy
Δy
'
=
v
∂u ∂u ,
rΔt 2 ~2 1/ 2
i, j
hạng phi tuyến dạng u , v ... và cho dao động mực nước cùng bậc 1 + ~ (vi , j + Si , j )
∂x ∂y
Li , j
với độ sâu biển nhờ sự thay thế độ sâu trung bình biển bằng D + ζ .
trong đó dùng các ký hiệu
Khi tích phân bằng số hệ phương trình này người ta hay sử dụng hệ
D + Di , j +1 + ζ i , j + ζ i , j +1
~
lưới sai phân so le, trong đó các điểm tính ζ , u , v dịch chuyển so với Di , j = i , j
2
nhau một nửa bước tính. Trị số của độ cao mực nước ζ được tính tại tâm
D + Di +1, j + ζ i , j + ζ i +1, j
của ô chữ nhật, các trị số của u và v được tính tại các điểm giữa của các ~
Li , j = i , j
cạnh ô chữ nhật (hình 2.4). Trong hệ lưới này các đạo hàm theo trục x 2
và y trong các phương trình vi phân cũng được xấp xỉ bằng sai phân hữu
v +v +v +v
~
K i , j = i , j i , j +1 i−1, j i−1, j −1
hạn trung tâm đối với những điểm tính bên trong vùng tính, sai phân hữu
4
hạn một chiều (tiến hoặc lùi) đối với các điểm trên biên cứng hoặc biên
u + ui , j −1 + ui −1, j + ui +1, j −1
~
lỏng. Còn đạo hàm thời gian được xấp xỉ bằng sai phân hữu hạn một Si , j = i , j
4
chiều tiến. Ở các điểm thuộc biên cứng kinh tuyến u = 0 và ở các điểm
thuộc biên cứng vĩ tuyến v = 0 theo điều kiện biên tương tự (2.25). các dấu phảy phía trên đại lượng chỉ trị số ở bước tính tiếp sau một thời
47
- gian Δt (bước thời gian) của đại lượng tương ứng. các tác giả Việt Nam cũng chủ yếu sử dụng mô hình này để nghiên cứu
những dạng dao động mực nước nguồn gốc khác nhau cho các vùng của
Trên đây mới chỉ giới thiệu một phương pháp giải số trị đơn giản
biển Đông. Thí dụ, bằng mô hình này Đỗ Ngọc Quỳnh (1982) [15] đã
nhất đối với hệ phương trình sóng dài trong nước nông dựa trên phương
nghiên cứu đặc điểm nước dâng trong bão ở biển Đông, Bùi Hồng Long
pháp sai phân hữu hạn và sử dụng sơ đồ hiện. Tính đơn giản của sơ đồ
(1987) [13] và Nguyễn Thọ Sáo (1988) [17] khảo sát những đặc điểm dao
giải này chủ yếu là ở chỗ những trị số của các hàm cần tìm tại mỗi điểm
động triều ở vịnh Bắc Bộ và toàn biển Đông nói chung, Phạm Văn Huấn
tính ở bước thời gian sau được tính chỉ dựa theo những trị số đã tính được
(1991) [12] khảo sát dao động tự do và dao động mùa do gió mùa của
của chúng ở bước tính trước và những trị số trên biên, chứ không phụ
mực nước ở biển Đông. Trong khuôn khổ đề tài cấp nhà nước "Thủy triều
thuộc vào chính những trị số cần tính tại bước tính đang xét của những
và sự dâng lên của mực nước biển Đông" (1991-1995) do Nguyễn Ngọc
điểm xung quanh. Do đó không đòi hỏi phải lập và giải hệ phương trình
Thụy làm chủ nhiệm, tập thể các tác giả như Đỗ Ngọc Quỳnh, Phạm Văn
đại số tuyến tính để tính đồng thời các trị số của các hàm chưa biết.
Ninh, Nguyễn Việt Liên, Đinh Văn Mạnh [16], Lê Trọng Đào, Nguyễn
Thọ Sáo cũng sử dụng mô hình vừa giới thiệu với những sơ đồ giải số trị
khác nhau để nghiên cứu thủy triều và dòng triều chi tiết cho vùng biển
này.
Hình 2.4. Vị trí các điểm tính ζ , u , và v trên lưới so le
Hiện nay mô hình dao động mực nước tổng cộng trên đây với những
sơ đồ giải số trị khác nhau là công cụ chủ yếu dùng để tính toán thủy
triều, nước dâng, dao động dâng rút do gió hoặc dao động tổng cộng của
mực nước trong các biển ven, những thủy vực nước nông ven biển và
vùng cửa sông (xem German, Levikov (1988), Koutitas (1988) [7]) trong
khuôn khổ bài toán truyền sóng dài hai chiều. Trong những năm gần đây
48
nguon tai.lieu . vn
