©2004 Przemysław Baran – www.ar.krakow.pl\~pbaran

1

Przykład rozwiązania równania Laplace’a metodą różnic skończonych

Nasze zadanie polegać będzie na znalezieniu wartości ciśnień i prędkości wewnątrz
nasypu (Rys. 1) piętrzącego wodę w warunkach filtracji ustalonej.

i

j

0 1 2 3 4 5 6 7 8 9 10

0

1

2

3

4

5

6

7

8

9

7,20

7,20

7,20

7,20

7,20

7,20

7,20

7,20

7,20

7,20

1,60

1,60

1,60

7,20 7,20

6,40

5,60

4,80

4,00

3,20

2,40

1,60

Rys. 1 Nasyp piętrzący wodę – przykład obliczeniowy dla MRS

Nasyp (k = 1x10

-4

[m/s]) jest prostopadłościanem o wymiarach 8m x 8m x 1m.

Spoczywa na podłożu nieprzepuszczalnym. Piętrzenie wody górnej 7,2m i dolnej
1,6m. Zakładamy, że znane jest położenie krzywej depresji w nasypie. Interesujący
nas obszar znajduje się pod krzywą filtracji i tam przeprowadzamy dyskretyzację
siatką kwadratową. Dyskretyzować należy w sposób, który pozwoli umieścić węzeł w
miejscu, gdzie znana jest wartość potencjału brzegowego. Jeżeli nie jest to możliwe,
stosowanie cytowanych wzorów relaksacyjnych jest niepoprawne. Ustalamy poziom
porównawczy u podstawy nasypu. Wartości potencjałów z lewej strony są równe
wysokości zwierciadła wody w piezometrze mierzonej od poziomu porównawczego.
Podobnie strona prawa. Na krzywej filtracji wartości potencjałów można wyznaczyć
jako wysokości zwierciadła, mierzone od poziomu porównawczego (jest to pewne
uproszczenie). Przy podstawie nasypu wprowadzamy warunek brzegowy dla
prędkości pionowej (v

z

= 0). Wartości potencjałów brzegowych są opisane kolorem

żółtym. Węzły koloru czerwonego posiadają nieznany potencjał i stanowią one
pierwszy przedmiot naszych obliczeń.

9

7,200

7,200

7,200

8

7,200

6,400

7

7,200

5,600

6

7,200

4,800

5

7,200

4,000

4

7,200

3,200

3

7,200

2,400

2

7,200

1,600

1,600

1

7,200

1,600

0

7,200

1,600

0

1

2

3

4

5

6

7

8

9

10

Tab. 1 Zestawienie numerów węzłów i ich wartości brzegowych

©2004 Przemysław Baran – www.ar.krakow.pl\~pbaran

2

W tym celu zestawmy potrzebne dane w formie tabelarycznej. Potrzebujemy numer
węzła w pionie i poziomie, oraz wartości potencjału na węzłach brzegowych (Tab. 1).
Następnym krokiem jest założenie wartości potencjałów wewnątrz obszaru analizy. I
tutaj, zgodnie z tym co było poruszane w części teoretycznej, istnieje w zasadzie
dowolność. Z tym, że im bardziej „realne” wartości zostaną przyjęte, tym szybciej
uzyskamy zbieżność iteracji. Proponuję zatem, aby ustalić wartości potencjałów na
zasadzie ich liniowego rozkładu pomiędzy wartościami brzegowymi. Zatem
wyjściowy zbiór danych przybierze poniższą postać (Tab. 2):

9

7,200

7,200

7,200

8

7,200

6,933

6,667

6,400

7

7,200

6,800

6,400

6,000

5,600

6

7,200

6,720

6,240

5,760

5,280

4,800

5

7,200

6,667

6,133

5,600

5,067

4,533

4,000

4

7,200

6,629

6,057

5,486

4,914

4,343

3,771

3,200

3

7,200

6,600

6,000

5,400

4,800

4,200

3,600

3,000

2,400

2

7,200

6,578

5,956

5,333

4,711

4,089

3,467

2,844

2,222

1,600

1,600

1

7,200

6,640

6,080

5,520

4,960

4,400

3,840

3,280

2,720

2,160

1,600

0

7,200

6,640

6,080

5,520

4,960

4,400

3,840

3,280

2,720

2,160

1,600

0

1

2

3

4

5

6

7

8

9

10

Tab. 2 Dane wyjściowe do obliczeń iteracyjnych

Teraz przystępujemy do wyznaczenia wartości potencjałów w pierwszym kroku
iteracyjnym, zgodnie z cytowaną w części teoretycznej formułą Jakobiego. Dla węzła
(3; 3) otrzymamy:

( )

( )

( )

( )

( )

( )

405

,

5

4

800

,

4

000

,

6

333

,

5

486

,

5

3

;

3

4

3

;

4

3

;

2

2

;

3

4

;

3

3

;

3

1

0

0

0

0

1

=

+

+

+

=

+

+

+

=

H

H

H

H

H

H

W Tab. 3 zestawiono wyniki pierwszej iteracji. Wartości zaznaczone kolorem żółtym
są wyliczane na podstawie warunku brzegowego dla prędkości (v

z

= 0). Dla węzła

(1; 0) otrzymamy:

( )

( )

( )

( )

( )

640

,

6

4

200

,

7

080

,

6

640

,

6

2

0

;

1

4

0

;

0

0

;

2

1

;

1

2

0

;

1

1

0

0

0

1

=

+

+

⋅

=

+

+

⋅

=

H

H

H

H

H

Przy założonej dokładności obliczeń równej 0,001, zbieżność uzyskano przy 34
iteracji (Tab. 4).

9

7,200

7,200

7,200

8

7,200

6,967

6,733

6,400

7

7,200

6,813

6,427

6,040

5,600

6

7,200

6,727

6,253

5,780

5,307

4,800

5

7,200

6,670

6,141

5,611

5,082

4,552

4,000

4

7,200

6,631

6,062

5,493

4,924

4,355

3,786

3,200

3

7,200

6,602

6,003

5,405

4,806

4,208

3,610

3,011

2,400

2

7,200

6,599

5,998

5,397

4,796

4,194

3,593

2,992

2,391

1,600

1,600

1

7,200

6,624

6,049

5,473

4,898

4,322

3,747

3,171

2,596

2,020

1,600

0

7,200

6,640

6,080

5,520

4,960

4,400

3,840

3,280

2,720

2,160

1,600

0

1

2

3

4

5

6

7

8

9

10

Tab. 3 Wyniki 1-szej iteracji

©2004 Przemysław Baran – www.ar.krakow.pl\~pbaran

3

9

7,200

7,200

7,200

8

7,200

7,021

6,791

6,400

7

7,200

6,892

6,544

6,113

5,600

6

7,200

6,804

6,381

5,910

5,381

4,800

5

7,200

6,743

6,269

5,763

5,217

4,630

4,000

4

7,200

6,702

6,191

5,659

5,097

4,503

3,871

3,200

3

7,200

6,673

6,138

5,586

5,013

4,414

3,784

3,117

2,400

2

7,200

6,655

6,103

5,540

4,960

4,360

3,736

3,084

2,390

1,600

1,600

1

7,200

6,644

6,084

5,514

4,931

4,332

3,719

3,096

2,478

1,916

1,600

0

7,200

6,641

6,078

5,506

4,921

4,324

3,716

3,104

2,511

1,985

1,600

0

1

2

3

4

5

6

7

8

9

10

Tab. 4 Wyniki 34-tej iteracji

Jak widać, procedura iteracyjna Jakobiego dała zgodność dopiero przy 34 iteracji.
Problem tutaj prezentowany opiera się na niewielkiej siatce (75 punktów). Przy
gęstszej dyskretyzacji mogłoby się okazać, że potrzeba o wiele większej ilości
przybliżeń, a co za tym idzie i czasu

1

. Poniżej zaprezentuję tą samą iterację, ale

metodą nad-relaksacji. Zgodnie z formułą, część danych węzłowych bierze się z
poprzedniej iteracji, a część z bieżącej. Tabela wyjściowa nie zmienia się (Tab. 2),
stąd np. dla węzła (1; 8) w pierwszej iteracji otrzymamy:

( )

( )

( )

( )

[

]

( ) (

)

(

)

993

,

6

4

200

,

7

667

,

6

825

,

6

200

,

7

5

,

1

933

,

6

5

,

1

1

8

;

1

5

,

1

4

9

;

1

8

;

2

7

;

1

8

;

0

)

8

;

1

(

)

1

(

)

8

;

1

(

1

0

0

1

1

0

1

=

+

+

+

⋅

+

⋅

−

=

=

+

+

+

⋅

+

⋅

−

=

H

H

H

H

H

H

H

ω

ω

ω

Patrząc na wzór, zobaczymy, że z bieżącej iteracji potrzebujemy wartości H(i,j-1) i
H(i-1,j). Zatem, aby móc przeprowadzić obliczenia musimy znać te wartości. Dobrze
jest zacząć od węzłów, w których cytowane wartości są narzucone przez warunki
brzegowe. Warunek brzegowy dla metody nad-relaksacji przy v

z

= 0 można opisać:

( )

(

) ( )

(

)

(

)

(

)

[

]

4

,

1

1

,

2

,

1

,

1

,

1

1

+

+

−

+

+

⋅

+

+

⋅

⋅

⋅

−

=

m

m

m

m

m

j

i

H

j

i

H

j

i

H

j

i

H

j

i

H

ω

ω

I tak dla węzła (1; 0) w pierwszej iteracji otrzymamy:

( ) (

) ( )

( )

( )

( )

[

]

( ) (

)

[

]

640

,

6

4

080

,

6

640

,

6

2

200

,

7

5

,

1

640

,

6

5

,

1

1

0

;

1

5

,

1

4

0

;

2

1

;

1

2

0

;

0

0

;

1

1

0

;

1

1

0

0

1

0

1

=

+

⋅

+

⋅

+

⋅

−

=

=

+

⋅

+

⋅

+

⋅

−

=

H

H

H

H

H

H

ω

ω

ω

W Tab. 5 pokazano wyniki pierwszej iteracji. Zbieżność uzyskano już przy 14 (!)
iteracji (Tab. 6).

1

Pamiętam z XIV Konferencji w Korbielowie pt. „Metody numeryczne w projektowaniu i analizie konstrukcji

hydrotechnicznych”, jak pewien naukowiec powiedział, iż amerykańscy numerycy w przypadku problemów z
długim czasem obliczeń, nie zmieniają procedury tylko komputery...

©2004 Przemysław Baran – www.ar.krakow.pl\~pbaran

4

9

7,200

7,200

7,200

8

7,200

6,993

6,812

6,400

7

7,200

6,825

6,462

6,105

5,600

6

7,200

6,733

6,273

5,817

5,363

4,800

5

7,200

6,675

6,156

5,639

5,124

4,610

4,000

4

7,200

6,636

6,078

5,522

4,968

4,414

3,861

3,200

3

7,200

6,611

6,028

5,448

4,868

4,289

3,710

3,131

2,400

2

7,200

6,601

6,007

5,413

4,820

4,226

3,633

3,039

2,445

1,600

1,600

1

7,200

6,617

6,025

5,429

4,833

4,236

3,638

3,041

2,444

1,846

1,600

0

7,200

6,640

6,080

5,520

4,960

4,400

3,840

3,280

2,720

2,160

1,600

0

1

2

3

4

5

6

7

8

9

10

Tab. 5 Wyniki 1-szej iteracji metodą nad-relaksacji

9

7,200

7,200

7,200

8

7,200

7,022

6,792

6,400

7

7,200

6,894

6,547

6,115

5,600

6

7,200

6,807

6,387

5,915

5,385

4,800

5

7,200

6,748

6,277

5,772

5,226

4,635

4,000

4

7,200

6,707

6,202

5,672

5,110

4,513

3,877

3,200

3

7,200

6,680

6,151

5,603

5,031

4,429

3,795

3,122

2,400

2

7,200

6,662

6,118

5,559

4,981

4,379

3,751

3,094

2,395

1,600

1,600

1

7,200

6,653

6,100

5,535

4,953

4,354

3,737

3,109

2,486

1,919

1,600

0

7,200

6,650

6,093

5,526

4,945

4,347

3,735

3,118

2,519

1,989

1,600

0

1

2

3

4

5

6

7

8

9

10

Tab. 6 Wyniki 14-tej iteracji metodą nad-relaksacji

Pojawiły się pewne rozbieżności między wartościami potencjałów wyliczonych
metodą Jakobiego, a metodą nad-relaksacji. Jest to objaw typowy dla metod
przybliżonych. W dalszej części obliczeń będziemy się koncentrować na wynikach z
metody nad-relaksacji.

Po wyliczeniu wartości potencjałów w punktach analizowanego obszaru. Można
przystąpić do wyznaczenia prędkości wody. Na prędkość w danym węźle składają
się dwa wektory prędkości: pionowy v

z

i poziomy v

x

. Aby je wyznaczyć skorzystamy z

cytowanego w części teoretycznej prawa Darcy, wyrażonego poprzez różnice
skończone:

( )

(

)

(

)

( )

(

)

(

)

x

j

i

H

j

i

H

k

j

i

v

x

j

i

H

j

i

H

k

j

i

v

z

x

∆

⋅

−

−

+

⋅

−

=

∆

⋅

−

−

+

⋅

−

=

2

1

,

1

,

,

2

,

1

,

1

,

I tak, dla węzła (5; 4) prędkości wyniosą:

( )

( )

( )

( )

( )

( )

( )

(

) (

)

s

m

v

v

v

s

m

H

H

v

s

m

H

H

v

m

x

z

x

xz

z

x

/

10

815

,

7

10

280

,

1

10

710

,

7

4

;

5

/

10

280

,

1

6

,

1

429

,

4

635

,

4

10

1

8

,

0

2

3

;

5

5

;

5

10

1

4

;

5

/

10

710

,

7

6

,

1

110

,

5

877

,

3

10

1

8

,

0

2

4

;

4

4

;

6

10

1

4

;

5

8

,

0

5

2

5

2

5

2

2

5

4

4

5

4

4

−

−

−

−

−

−

−

−

−

⋅

=

⋅

−

+

⋅

=

+

=

⋅

−

=

−

⋅

⋅

−

=

⋅

−

⋅

⋅

−

=

⋅

=

−

⋅

⋅

−

=

⋅

−

⋅

⋅

−

=

=

∆

©2004 Przemysław Baran – www.ar.krakow.pl\~pbaran

5

A kierunek wektora prędkości dla tego samego węzła:

( )

°

−

=

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

⋅

⋅

−

=

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

=

⇒

=

−

−

−

−

428

,

9

10

710

,

7

10

280

,

1

tg

4

;

5

tg

tg

5

5

1

1

α

α

α

x

z

x

z

v

v

v

v

Poniżej zestawiono prędkości w nasypie (Tab. 7-Tab. 9) oraz kąty nachylenia
wektorów prędkości (Tab. 10).

8

2,549E-05 3,885E-05

7

4,081E-05 4,865E-05 5,919E-05

6

5,083E-05 5,576E-05 6,259E-05 6,969E-05

5

5,766E-05 6,096E-05 6,572E-05 7,112E-05 7,661E-05

4

6,236E-05 6,471E-05 6,822E-05 7,245E-05 7,710E-05 8,207E-05

3

6,554E-05 6,730E-05 7,000E-05 7,337E-05 7,726E-05 8,171E-05 8,721E-05

2

6,756E-05 6,895E-05 7,109E-05 7,378E-05 7,685E-05 8,031E-05 8,480E-05 9,342E-05

1

6,868E-05 6,985E-05 7,164E-05 7,381E-05 7,606E-05 7,787E-05 7,827E-05 7,440E-05 5,536E-05

0

6,904E-05 7,014E-05 7,180E-05 7,378E-05 7,569E-05 7,685E-05 7,599E-05 7,057E-05 5,749E-05

vx

1

2

3

4

5

6

7

8

9

Tab. 7 Wartości [m/s] prędkości poziomych v

x

8

-1,913E-05 -4,081E-05

7

-1,340E-05 -2,534E-05 -3,031E-05

6

-9,116E-06 -1,685E-05 -2,143E-05 -2,339E-05

5

-6,215E-06 -1,153E-05 -1,517E-05 -1,716E-05 -1,793E-05

4

-4,215E-06 -7,873E-06 -1,055E-05 -1,215E-05 -1,280E-05 -1,279E-05

3

-2,774E-06 -5,202E-06 -7,014E-06 -8,073E-06 -8,343E-06 -7,825E-06 -6,582E-06

2

-1,677E-06 -3,147E-06 -4,231E-06 -4,789E-06 -4,672E-06 -3,586E-06 -8,309E-07 5,358E-06

1

-7,885E-07 -1,477E-06 -1,972E-06 -2,182E-06 -1,968E-06 -1,017E-06 1,487E-06 7,795E-06 2,433E-05

0

0,000E+00 0,000E+00 0,000E+00 0,000E+00 0,000E+00 0,000E+00 0,000E+00 0,000E+00 0,000E+00

vz

1

2

3

4

5

6

7

8

9

Tab. 8 Wartości [m/s] prędkości pionowych v

z

8

3,187E-05 5,634E-05

7

4,295E-05 5,486E-05 6,650E-05

6

5,164E-05 5,825E-05 6,616E-05 7,351E-05

5

5,800E-05 6,204E-05 6,745E-05 7,316E-05 7,868E-05

4

6,250E-05 6,519E-05 6,903E-05 7,346E-05 7,815E-05 8,306E-05

3

6,559E-05 6,750E-05 7,035E-05 7,382E-05 7,771E-05 8,209E-05 8,746E-05

2

6,758E-05 6,902E-05 7,122E-05 7,394E-05 7,699E-05 8,039E-05 8,480E-05 9,357E-05

1

6,869E-05 6,987E-05 7,167E-05 7,385E-05 7,608E-05 7,788E-05 7,828E-05 7,481E-05 6,047E-05

0

6,904E-05 7,014E-05 7,180E-05 7,378E-05 7,569E-05 7,685E-05 7,599E-05 7,057E-05 5,749E-05

vxz

1

2

3

4

5

6

7

8

9

Tab. 9 Wartości [m/s] prędkości wypadkowych v

xz

8

-36,885

-46,412

7

-18,175

-27,515

-27,117

6

-10,167

-16,818

-18,900

-18,551

5

-6,151

-10,707

-12,993

-13,563

-13,172

4

-3,867

-6,937

-8,790

-9,516

-9,428

-8,858

3

-2,423

-4,420

-5,722

-6,279

-6,163

-5,470

-4,316

2

-1,422

-2,613

-3,406

-3,714

-3,479

-2,557

-0,561

3,282

1

-0,658

-1,211

-1,577

-1,694

-1,482

-0,748

1,088

5,981

23,729

0

0,000

0,000

0,000

0,000

0,000

0,000

0,000

0,000

0,000

alfa

1

2

3

4

5

6

7

8

9

Tab. 10 Wartości [

°] kątów

α nachylenia wektorów prędkości

©2004 Przemysław Baran – www.ar.krakow.pl\~pbaran

6

Pozostaje jeszcze ustalenie wartości ciśnień porowych u. Aby tego dokonać
potrzebujemy, oprócz wartości potencjału w rozpatrywanym punkcie, wysokość
danego punktu mierzoną od założonego poziomu porównawczego. I tak dla węzła
(3; 3) otrzymamy następującą wartość ciśnienia porowego:

( )

( )

( )

[

]

( )

[

]

( )

kPa

u

u

H

z

u

w

430

,

31

3

;

3

81

,

9

604

,

5

40

,

2

3

;

3

3

;

3

3

;

3

3

;

3

−

=

⋅

−

=

⋅

−

=

γ

Znak „minus” traktuje ciśnienie jako parcie (w przeciwieństwie do ssania – znak
„plus”). Analogia do ustrojów prętowych (np. kratownic), gdzie ściskanie jest „ujemne”
a rozciąganie „dodatnie”. Poniżej (Tab. 11) zestawiono obliczone wartości ciśnień
porowych.

9

7,20

0,000

0,000

8

6,40

-6,097

-3,847

0,000

7

5,60

-12,694

-9,291

-5,057

0,000

6

4,80

-19,690

-15,565

-10,938

-5,741

0,000

5

4,00

-26,959

-22,341

-17,390

-12,025

-6,227

0,000

4

3,20

-34,411

-29,452

-24,254

-18,744

-12,882

-6,643

0,000

3

2,40

-41,993

-36,802

-31,430

-25,815

-19,913

-13,688

-7,088

0,000

2

1,60

-49,672

-44,331

-38,849

-33,173

-27,268

-21,110

-14,663

-7,800

0,000

1

0,80

-57,426

-52,004

-46,462

-40,759

-34,876

-28,822

-22,653

-16,537

-10,975

0

0,00

-65,244

-59,796

-54,235

-48,526

-42,655

-36,647

-30,592

-24,719

-19,515

u

z

1

2

3

4

5

6

7

8

9

Tab. 11 Wartości [kPa] ciśnień porowych u