G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

1 

 

 

SPIS TREŚCI 

 

 

1. Metody przybliżone w mechanice konstrukcji 

    2 

2. Metoda Różnic Skończonych  

  9 

3. Metoda Elementów Brzegowych  

17 

 

3.1 

MEB dla równania Poissona  

17 

 

3.2 

Zagadnienia teorii sprężystości (poza programem) 

21 

4. Koncepcja MES na przykładzie równania Poissona 

35 

5. MES w analizie konstrukcji prętowych 

38 

5.1.Belki 

 

38 

5.2. Pręty rozciągane i skręcane. Sprężyny 

52 

5.3. Kratownice i ramy płaskie 

59 

5.4. Przestrzenne kratownice i ramy 

63 

6. Dwuwymiarowe i trójwymiarowe zadania teorii sprężystości 

66 

 

 6.1. Element skończony trójkątny CST  (constant strain triangle) 

   77 

 

 6.2. Izoparametryczny 8-węzłowy element skończony 

   84 

 

6. 3. Całkowanie numeryczne 

 

89 

7. MES w praktyce inżynierskiej 

 

   91 

7.1.Uwagi o dokładności obliczeń MES ... 

91 

7.2. Wykorzystanie profesjonalnych systemów obliczeniowych ..... 

 96 

 

Literatura   

 

 

 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

2 

1.  METODY PRZYBLIŻONE W MECHANICE 

KONSTRUKCJI  

 

Konstrukcje  odkształcalne  mogą  być  badane  doświadczalnie  lub  metodami 

teoretycznymi; analitycznymi i numerycznymi. Poważną wadą analiz eksperymentalnych jest 

ich  koszt  i  czasochłonność.  Dotyczy  to  zarówno  większości  badań  modelowych  jak  i  badań 

obiektów rzeczywistych. Pracochłonność metod doświadczalnych jest szczególnie odczuwana 

w  trakcie  prac  projektowych,  gdzie  analizie  poddawane  są  różne  warianty  konstrukcji. 

Dlatego  rozwój  metod  teoretycznych,  przede  wszystkim  metod  numerycznych,  wpływał 

zawsze na jakość projektowanych konstrukcji inżynierskich. Przykładem może być postęp w 

konstrukcjach statków, dźwigów, wysokich budynków, samochodów, samolotów. 

Badania  teoretyczne  polegają  na  sformułowaniu  odpowiedniego  opisu  matematycznego  i 

następnie  rozwiązaniu  postawionego  problemu.  Niestety,  dla  bardzo  wielu  praktycznych 

problemów  mechaniki  konstrukcji  można  zbudować  wiarygodny  model  matematyczny,  ale 

nie  są  znane  odpowiednie  ścisłe  rozwiązania  analityczne.  Prostym  przykładem  takich 

zagadnień  jest  wyznaczanie  współczynników  koncentracji  naprężeń.  Tylko  w  bardzo 

szczególnych przypadkach  znane są rozwiązania dokładne. 

 

 

Rzeczywiste zjawisko 

Model matematyczny 

ci

ą

gły 

Model 

dyskretny 

Rzeczywisty wynik 

Rozwi

ą

zanie 

ś

cisłe  

modelu matematycznego 

Rozwi

ą

zanie dokładne 

modelu dyskretnego 

Prawa fizyki 
Własno

ś

ci materiałowe, geometria, 

Dyskretyzacja 

Aproksymacja  

Realizacja oblicze

ń

 

Wynik numeryczny 

Warunki brzegowe 

M

 

E

 

T

 

O

 

D

 

A

 

P

 

R

 

Z

 

Y

 

B

 

L

 

I

 

Ż

 

O

 

N

 

A

 

 

 

Rys. 1. Rozwiązywanie zagadnień analizy ośrodków ciągłych metodami przybliżonymi. 

Schemat postępowania 

 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

3 

Dlatego równolegle do rozwoju analitycznych metod znajdowania rozwiązań ścisłych 

były  opracowywane  i  doskonalone  metody  przybliżone.  Większość  zadań  analizy 

wytrzymałościowej konstrukcji rozwiązuje się współcześnie za pomocą metod przybliżonych, 

przy wykorzystaniu obliczeń komputerowych. 

Ogólny schemat postępowania przy zastosowaniu metod przybliżonych przedstawiony jest na 

rysunku 1. 

Pierwszym  krokiem  na  drodze  poszukiwania  rozwiązania  jest  budowa  modelu 

matematycznego.  Potrzebna  jest  do  tego  znajomość  odpowiednich  praw  fizyki, 

sformalizowany opis własności materiałowych, kształtu konstrukcji, jej warunków podparcia i 

obciążenia.  Dla  każdego  rzeczywistego  problemu  możemy  zbudować  wiele  różnych  modeli 

matematycznych.  Bardzo  prostą  ilustracją  może  być  tu  zadanie  określenia    ugięcia  pod 

wpływem sił ciężkości  swobodnie podpartej drewnianej deski (rys. 2). 

 

 

a) model belki

 

b) model płyty

 

c) model trójwymiarowy bryły

 

q

 

0

 

N

 

m

 

p

 

0

 

m

 

N

 

2

 

N

 

m

 

0

 

3

 

 

 

Rys. 2. Różne modele w analizie wytrzymałościowej zginanej deski 

 

 

Dla  takiego  zagadnienia  można  przedstawić  typowy  jednowymiarowy  model 

matematyczny 

belki, 

dwuwymiarowy 

model 

zginanej 

płyty 

lub 

pełny 

model 

trójwymiarowego  zadania  mechaniki  ciał  odkształcalnych.  Własności  materiału  i  równania 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

4 

opisujące  konstrukcję  mogą  w  najprostszym  przypadku  zakładać  izotropową  liniową 

sprężystość.  Mogą  również  uwzględniać  anizotropię,  lepkosprężystość  czy  duże  ugięcia  i  w 

konsekwencji  nieliniowość  zjawisk.  Również  warunki  brzegowe  można  przyjąć  w  postaci 

uproszczonej  lub  uwzględnić  na  przykład  zjawiska  kontaktowe.  Przykład  ten  pokazuje,  że 

zwykle  nie  ma  jednego,  najlepszego  modelu  obliczeniowego  dla  analizy  wytrzymałości 

elementów konstrukcyjnych. Właściwy model obliczeniowy zależy od celu analizy, wymagań 

stawianych  konstrukcji,  żądanej  dokładności  wyników,  dostępności  danych  materiałowych, 

ale także - w dużym stopniu - od dostępnych narzędzi obliczeniowych.  

Ponadto  opis  matematyczny  wybranego  modelu  można  przedstawić  w  różnych 

sformułowaniach,  na  przykład  równań  różniczkowych,  odpowiadających  im  równań 

całkowych  lub  w  ujęciu  wariacyjnym,  w  postaci  problemu  minimalizacji  odpowiedniego 

funkcjonału  [5].  Budowa  modelu  matematycznego  stanowi  najważniejszy  element  analizy 

obliczeniowej,  w  którym  bardzo  trudno  zastąpić  bezpośrednią  decyzję  człowieka  przez 

sformalizowane algorytmy postępowania. 

W  metodach  przybliżonych  problem  poszukiwania  nieznanych  funkcji  (opisujących 

pole  przemieszczeń,  odkształceń  i  naprężeń)  zastępujemy  przez  problem  poszukiwania 

skończonej liczby parametrów, za pomocą których można opisać – w pewnym przybliżeniu- 

poszukiwane funkcje.  

x

f(x)

x

1

x

2

b

x =a

h

0

x  =a+ih

i

h

h

h

f

1

f

2

f

0

3

f

f

i

 

Rys. 3. Dyskretyzacja i aproksymacja na przykładzie funkcji jednej zmiennej 

 

 

Dokonujemy  tego  poprzez  dyskretyzację,  czyli  wybór  skończonej  liczby  parametrów 

opisujących w przybliżeniu analizowany problem oraz aproksymację poszukiwanych funkcji 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

5 

za  pomocą  innych  z  góry  założonych  prostych  funkcji,  które  uzależnione  są  od  wybranych 

parametrów. 

Najprostszym  przykładem  dyskretyzacji  i  aproksymacji  jest  numeryczna  reprezentacja 

dowolnej  funkcji  f(x)  w  przedziale  <a,b>.  Zrealizować  możemy  to  przez  podział  przedziału 

<a,b> na n  równych podprzedziałów o długości h=(b-a)/n.  

Funkcję  ciągłą  f(x)  reprezentować  będzie  wtedy  zbiór  (n+1)  wartości  f(a+ih),  i=0,1,2,  n. 

Przybliżone wartości funkcji dla dowolnej wartości x z przedziału <a,b> wyznaczać możemy 

za pomocą aproksymacji funkcją w podprzedziałach stałą (schodkową), liniową (łamaną) lub 

na przykład wielomianowymi funkcjami sklejanymi (rys. 3). 

Metody przybliżone analizy ośrodków ciągłych stanowią intensywnie rozwijaną część 

współczesnej  matematyki.  Wiążą  się  z  problemami  analizy  funkcjonalnej,  metod 

numerycznych  i  rachunku  wariacyjnego.  Prezentuje  się  je  czasem  jako  szczególne  warianty 

tak  zwanych  technik  residuów  ważonych  dla  przybliżonego  rozwiązywania  zagadnień 

opisywanych przez równania różniczkowe cząstkowe [5,7]. 

W naukach technicznych metody przybliżone opisywane są zazwyczaj od strony zastosowań, 

jako  wyspecjalizowane  procedury  obliczeniowe  dla  określonych  zagadnień,  na  przykład  pól 

deformacji  sprężystej,  pól  temperatury,  pół  elektrycznych  i  magnetycznych.  Istnieje  wiele 

metod przybliżonych, stosowanych od dawna w mechanice ciała odkształcalnego i mechanice 

konstrukcji, np. metoda Ritza, Galerkina, Trefftza, kollokacji, Kantorowicza [2,6]. Metody te 

różnią się postacią modelu matematycznego, sposobem dyskretyzacji i aproksymacji, a także 

stosowanym  algorytmem  obliczeń.  Jednak  współcześnie  dominują  te  metody  przybliżone, 

które  charakteryzuje  łatwość  pełnej  automatyzacji  i  które  rozwinęły  się  wraz    z  rozwojem 

technik  komputerowych.  Najbardziej  znane  z  nich  to  metoda  różnic  skończonych  (MRS), 

metoda elementów skończonych (MES) i metoda elementów brzegowych (MEB). Nazywamy 

je  metodami  numerycznymi  lub  komputerowymi.  Każda  z  tych  metod  ma  bardzo  bogatą 

literaturę i prezentowana może być w wielu odmiennych sformułowaniach. 

Metoda  różnic  skończonych  bezpośrednio  wykorzystuje  model  matematyczny  w 

postaci  równań  różniczkowych  opisujących  rozwiązywane  zagadnienie  [4,5].  Dyskretyzacja 

polega  na  ustaleniu  w  analizowanym  obszarze  siatki  węzłów,  np.  w  przypadku 

dwuwymiarowym  prostokątnej  lub  trójkątnej.  Przyjmujemy,  że  w  każdym  węźle  pochodne 

można  przybliżyć  przez  tak  zwane  ilorazy  różnicowe,  czyli  wyrażenia  algebraiczne 

uzależnione  od  wartości  funkcji  w  węzłach  przyległych.  Równanie  różniczkowe  w  każdym 

węźle zastępowane jest przez odpowiednie równanie algebraiczne wiążące wartości funkcji w 

najbliższych  węzłach.  W  rezultacie  otrzymujemy  układ  równań,  których  liczba  odpowiada 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

6 

liczbie  węzłów.  Niewiadomymi  są  wartości  poszukiwanej  funkcji  w  tych  węzłach.  Warunki 

brzegowe  przedstawiane  są  jako  dodatkowe  związki  łączące  wartości  funkcji  w    węzłach 

przylegających do konturu analizowanego obszaru. 

W  metodzie  elementów  skończonych  [6,7]  zamiast  rozwiązywać  równania 

różniczkowe  poszukujemy  przybliżonego  rozwiązania  sformułowania  alternatywnego  w 

postaci  problemu  minimalizacji  odpowiedniego  funkcjonału

1

.  Zagadnieniami  minimalizacji 

funkcjonałów i ich związkami z odpowiednimi równaniami różniczkowymi zajmuje się dział 

matematyki  nazywany  rachunkiem  wariacyjnym.  W  MES  analizowany  obszar  dzielony  jest 

na podobszary, nazywane elementami skończonymi. W przypadku dwuwymiarowym mogą to 

być np. trójkąty lub czworokąty. Elementy skończone łączą się ze sobą w węzłach. Położenie 

węzłów elementu wyznacza jego kształt. Przebieg poszukiwanej funkcji u wewnątrz każdego 

elementu  wyznacza  się  za  pomocą  z  góry  ustalonych  funkcji  aproksymujących  nazywanych 

funkcjami kształtu: 

( )

( )

i

i

i

u x

N x u

=

∑

 

gdzie 

i

u   stanowią  wartości  funkcji  w  węzłach,  x   jest  wektorem  współrzędnych,  a 

i

N

  są 

funkcjami kształtu. 

Minimalizowany  funkcjonał  przedstawić  wtedy  można  jako  funkcję  wielu  zmiennych 

Niewiadomymi są zwykle wartości poszukiwanej funkcji w węzłach całego obszaru. Warunek 

minimum  prowadzi  do  układu  równań  liniowych,  którego  rozwiązanie  określa  nieznane 

wartości funkcji w węzłach. 

Metoda  elementów  brzegowych 

[5]  wymaga  znajomości  tak  zwanych  brzegowych 

równań  całkowych,  które  w  innej  formie  wyrażają  związki  znane  w  postaci  równań 

różniczkowych lub w postaci problemu minimalizacji odpowiedniego funkcjonału. Dla wielu 

zagadnień    tak  zwane  brzegowe  sformułowanie,  które  wyraża  całkowe  zależności  między 

wartościami  funkcji  i  jej  pochodnych  dla  punktów  brzegowych  jest  znane  z  analizy 

matematycznej. 

Dyskretyzacja polega na podziale brzegu na małe segmenty zwane elementami brzegowymi. 

W  przypadku  dwuwymiarowym  mogą  to  być  odcinki  prostych  lub  też  krzywych,  których 

kształt  opisany  jest  przez  wielomiany.  Zapisanie  dyskretnego  odpowiednika  równania 

całkowego dla każdego  węzła prowadzi do układu równań algebraicznych, z którego można 

wyznaczyć wartości poszukiwanej funkcji na brzegu. 

 

                                                           

1

 Metoda elementów skończonych może być prezentowana jako oparta na innych sformułowaniach] 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

7 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

 

Rys. 4. Ogólny schemat postępowania przy obliczeniach metodą różnic skończonych, metodą 

elementów brzegowych i metodą elementów skończonych 

 

 

Równania różniczkowe 

cząstkowe 

Budowa siatki węzłów i 

przyjęcie wybranych 

schematów różnicowych 

Całkowe równania brzegowe 

Problem minimalizacji 

funkcjonału 

Podział brzegu na segmenty 

(elementy brzegowe) i założenie 

odpowiednich funkcji aproksymują-

cych na elementach – funkcji 

kształtu 

Podział  obszaru na małe 

podobszary (elementy skończone) i 

przyjęcie odpowiednich funkcji 

aproksymujących na elementach –

funkcji kształtu 

Zastąpienie równań 

różniczkowych przez równania 

różnicowe dla kolejnych 

węzłów obszaru. 

Formowanie układu 

równańliniowych

 

Budowa dyskretnej reprezentacji 

równania całkowego dla kolejnych 

węzłów brzegu. 

Formowanie układu równań 

liniowych 

Budowa macierzy sztywności 

kolejnych elementów. 

Formowanie układu równań 

liniowych 

Rozwiązanie układu równań 

liniowych (macierz rzadka, 

pasmowa, zwykle 

symetryczna) 

Rozwiązanie układu równań 

liniowych (macierz pełna, 

niesymetryczna) 

Rozwiązanie układu równań 

liniowych (macierz rzadka, 

pasmowa, zwykle symetryczna) 

Modyfikacja układu równań – 

wprowadzenie warunków 

brzegowych 

Modyfikacja układu równań – 

wprowadzenie warunków 

brzegowych 

Modyfikacja układu równań – 

wprowadzenie warunków 

brzegowych 

Obliczenia uzupełniające, np. 

poszukiwanych funkcji i ich 

pochodnych w wybranych punktach 

obszaru 

Obliczenia uzupełniające, np. 

funkcji i jej pochodnych funkcji 

wewnątrz elementów skończonych 

Obliczenia uzupełniające, np. 

pochodnych poszukiwanych 

funkcji w węzłach 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

8 

Rysunek  4  przedstawia  ogólny  szkic  algorytmów  MRS  MEB  i  MES  dla  typowego 

problemu  liniowego,  na  przykład  dla  zagadnienia  Laplace’a,  zagadnienia  Poissona  lub 

liniowego zadania teorii sprężystości. 

Porównanie  powyższych  trzech  metod  jest  w  ogólnym  przypadku  bardzo  trudne. 

Istnieją  pewne  szczególne  zagadnienia,  dla  których  za  najbardziej  efektywną  można  uznać 

metodę  różnic  skończonych  lub  metodę  elementów  brzegowych.  Jednak  niewątpliwie 

najbardziej  uniwersalną  jest  metoda  elementów  skończonych.  Fakt  ten  potwierdza 

powszechność zastosowań inżynierskich programów komputerowych MES.  

Metoda elementów skończonych powstała w latach pięćdziesiątych XX  wieku jako technika 

obliczeniowa stosowana praktycznie, bez formalnych podstaw teoretycznych. Podwaliny dało 

opracowanie  metody  analizy  wytrzymałościowej  poprzez  podział  złożonych  konstrukcji 

nośnych  na  skończoną  liczbę  uproszczonych  elementów  składowych.  Związki,  które  musiał 

spełniać  tak  zbudowany  model  zapisywano  w  postaci  macierzowej.  Otrzymywano  układy 

równań  z  wieloma  niewiadomymi,  których  liczba  związana  była  z  dokładnością  przyjętego 

modelu. 

Czynnikiem wpływającym na gwałtowny rozwój MES był postęp w technice komputerowej i 

potrzeby  intensywnie  rozwijanego  wówczas  przemysłu  zbrojeniowego  i  lotniczego.  Analiza 

podstaw  teoretycznych  metody,  jej  systematyzowanie  i  zwiększanie  efektywności 

doprowadziło  do  korzystania  z  MES  jako  ogólnej  metody  obliczeń  problemów  ciągłych  w 

matematyce,  a  także  do  szerokich  zastosowań  w  praktyce  inżynierskiej    przy  analizie  pól 

naprężeń  odkształceń  i  przemieszczeń,  pól  elektrycznych,  magnetycznych,  termicznych, 

zagadnień przepływowych a także pól sprzężonych (np. zagadnień naprężeń temperaturowych 

czy problemów  elektryczno-cieplnych). 

 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

9 

2.  METODA RÓŻNIC SKOŃCZONYCH 

 

 

Metoda  różnic  skończonych  (MRS)  jest  jedną  z  najstarszych  metod  numerycznych 

rozwiązywania 

zagadnień 

opisywanych 

równaniami 

różniczkowymi 

cząstkowymi. 

Przyjmujemy  w  niej,  że  poszukiwana  funkcja  określona  jest  przez  zbiór  jej  wartości  w 

wybranych  punktach  (tzw.  węzłach),  dostatecznie  gęsto  rozłożonych  w  analizowanym 

obszarze.  Punkty  te  dobierane  są  zazwyczaj  tak.  że  tworzą  regularne  siatki,  z  których 

najprostszymi są siatka prostokątna w przypadku dwuwymiarowym lub prostopadłościenna w 

przypadku trójwymiarowym (rys. 5).  

x

y

y

a)

r

θ

1

2

3

4

5

6

h

g

0

x

b)

x

y

c)

x

y

d)

e)

l

 

 

Rys. 5. Przykłady regularnych siatek MRS: a) prostokątna, b) trójkątna, c) sześciokątna, d) w 

układzie biegunowym, e) prostopadłościenna 

 

W  metodzie  różnic  skończonych  operatory  różniczkowania  funkcji  zastępowane  są 

odpowiednimi  operatorami  różnicowymi.  Oznacza  to,  że  wartości  pochodnych  w  punkcie 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

10 

zastępowane  są  przez  odpowiednie  przyrosty  (różnice)  funkcji  w  sąsiadujących  węzłach. 

Zamieniamy  w  ten  sposób  równania  różniczkowe  na  równania  algebraiczne,  tak  zwane 

równania  różnicowe.  W  efekcie,  zamiast  równania  różniczkowego,  otrzymujemy  układ 

równań  z  niewiadomymi  będącymi  wartościami  funkcji  w  węzłach.  Każde  równanie 

odpowiada  odpowiedniemu  węzłowi  siatki.  Otrzymany  układ  równań  jest  modyfikowany  w 

wyniku  uwzględnienia  warunków  brzegowych.  Rozwiązaniem  układu  są  wartości 

poszukiwanej funkcji we wszystkich węzłach siatki. 

 

Rys. 6. Dwuwymiarowa siatka prostokątna metody różnic skończonych 

 

 

Załóżmy, że mamy do czynienia z zagadnieniem dwuwymiarowym we współrzędnych 

kartezjańskich xy i stosujemy typową siatkę prostokątną (rys. 6) opisaną przez formuły: 

 

,

,

i

o

k

o

x

x

ih

y

y

kg

= +

=

+

 

(1) 

gdzie 

,

o

o

x y

 są współrzędnymi wybranego punktu odniesienia. 

Przyjmijmy oznaczenie 

 

,

( ,

)

i k

i

k

u

u x y

=

, 

(2) 

gdzie  ( , )

u x y

 jest poszukiwaną funkcją. 

Możemy wówczas definiować różne schematy różnicowe, na przykład pierwszą pochodną

u

y

∂

∂

 

można zastąpić przez trzy różne operatory różnicowe: 

h

 

h

 

g

 

g

 

u

 

u

 

u

 

u

 

i,k+1

 

i,k

 

i,k-1

 

i-1,k

 

u

 

i+1,k

 

x

 

x

 

x

 

y

 

y

 

y

 

0

 

i

 

0

 

k

 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

11 

 

,

1

,

,

,

1

,

1

,

1

a)

,

b)

,

c)

.

2

i k

i k

i k

i k

i k

i k

u

u

u

u

y

y

g

u

u

u

u

y

y

g

u

u

u

u

y

y

g

+

−

+

−

−

∂

∆

≈

=

∂

∆

−

∂

∆

≈

=

∂

∆

−

∂

∆

≈

=

∂

∆

 

(3) 

Pierwszy  z  operatorów  (3)  nazywany  ilorazem  różnicowym  przednim,  drugi  ilorazem 

różnicowym  tylnym,  a  trzeci  ilorazem  różnicowym  centralnym.  Ilorazy  różnicowe 

odpowiadające  drugim  pochodnym  możemy  otrzymać  stosując  analogiczne  schematy 

powtórnie,  tym  razem  w  stosunku  do  pierwszych  pochodnych.  Typowe  formuły  to  na 

przykład: 

 

2

2

1,

,

1,

2

2

2

2

2

,

1

,

,

1

2

2

2

2

,

2

.

i

k

i k

i

k

i k

i k

i k

u

u

u

u

u

x

x

h

u

u

u

u

u

y

y

g

+

−

+

−

−

+

∂

∆

≈

=

∂

∆

−

+

∂

∆

≈

=

∂

∆

 

(4) 

Podobnie  uzyskamy  również  wyrażenia  na  ilorazy  różnicowe  odpowiadające  innym 

pochodnym np. 

 

4

4

2,

1,

,

1,

2,

4

4

4

4

6

4

i

j

i

j

i j

i

j

i

j

u

u

u

u

u

u

u

x

x

h

+

+

−

−

−

+

−

+

∂

∆

≈

=

∂

∆

 

(5) 

Szczegóły dalszego postępowania czytelnik znajdzie w przykładach. 

 

 

 

P

RZYKŁAD 

1 

Belka o długości l jest utwierdzona sztywno na obu końcach. Znane jest początkowe (dla chwili 

0

t

=

) ugięcie 

belki 

( )

o

w x

 oraz rozkład prędkości początkowej 

0

0

( , )

( )

( ,

0)

t

w x t

v x

v x t

t

=

∂

=

= =

∂

. Należy obliczyć jak 

zmienia się linia ugięcia belki dla 

0

t

>

. 

Rozwiązanie. 

Funkcja 

( , )

w x t

 opisująca ugięcie belki spełnia równanie różniczkowe. 

 

4

2

4

2

2

1

0

w

w

x

a

t

∂

∂

+

=

∂

∂

, 

(6) 

gdzie 

EJ

a

A

ρ

=

 oraz warunki brzegowe 

 

0

(0, )

( , )

0

dla

0,

dla

0,

x

x l

w

t

w l t

t

w

w

t

x

x

=

=

=

=

≥

∂

∂









=

≥









∂

∂









 

(7) 

i warunki początkowe 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

12 

 

0

0

0

( , 0)

( )

dla

(0, ),

( )

dla

(0, ).

t

w x

w x

x

l

w

v x

x

l

t

=

=

∈

∂





=

∈





∂





 

(8) 

Analizujemy obszar 

(0, ) (0, )

l

T

×

, gdzie 

T

 oznacza czas końcowy procesu 

(

0)

T

>

. Przyjmijmy prostokątną 

siatkę węzłów dla metody różnicowej 

 

,

.

i

j

x

ih

t

jg

=

=

 

Pochodnym w równaniu (6) odpowiadają ilorazy różnicowe: 

 

4

2,

1,

,

1,

2,

4

4

2

,

,

1

,

2

2

2

4

6

4

,

2

.

i

j

i

j

i j

i

j

i

j

i j

i j

i j

w

w

w

w

w

w

x

h

w

w

w

w

t

g

+

+

−

−

−

−

−

+

−

+

∆

=

∆

−

+

∆

=

∆

 

(9) 

W  ten  sposób  otrzymujemy  układ  równań  liniowych  stanowiący  przybliżoną,  różnicową  postać  równania 
różniczkowego (6): 

 

2

2

4

2,

1,

,

1,

2,

,

,

1

,

2

(

4

6

4

)

(

2

2

)

0,

0,1, 2,..., ,

0,1, 2,... .

i

j

i

j

i j

i

j

i

j

i j

i j

i j

a g w

w

w

w

w

h w

w

w

i

n

j

m

−

−

+

+

−

−

−

+

−

+

+

−

+

=

=

=

 

(10) 

Z warunków brzegowych mamy 

 

,

,

0

dla

0,...,

o j

n j

w

w

j

m

=

=

=

 

(11) 

i z warunków początkowych 

 

,

,1

( )

1, 2,...,

1,

( )

( )

1, 2,...,

1.

i o

o

i

i

o

i

o

i

w

w x

i

n

w

w x

v x

g

i

n

=

=

−

=

+

⋅

=

−

 

(12) 

Układ  równań  (10)  wraz  z  warunkami  (11)  i  (12)  może  być  rozwiązywany  dla  kolejnych  kroków  czasowych  
j=2,3,4....[10]. 

Rozwiązaniem jest wektor 

j

w

 opisujący w przybliżeniu linię ugięcia belki w chwili 

j

t

 

 

(

)

2

2

1

2

4

2

j

j

j

a

g

w

A I

w

w

n

−

−





⋅

=

+

−









, 

(13) 

gdzie 

I

 jest macierzą jednostkową, 

 

6

4

1

0

0

...

0

0

4

6

4

1

0

...

0

0

1

4

6

4

1

...

0

0

0

1

4

6

4 ...

0

0

0

0

0

0

0

...

6

4

0

0

0

0

0

...

4

6

A

−









−

−









−

−





=

−

−













−









−





⋮

⋮

⋮

⋮

⋮

⋮

⋮

⋮

, 

(14) 

jest macierzą nieosobliwą a wektor 

j

w

 ma postać 

 

2,

3,

2,

j

j

j

n

j

w

w

w

w

−













=

















⋮

. 

(15) 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

13 

Wektory 

j

w

  oblicza  się  korzystając  ze  znajomości  wektorów  linii  ugięcia  dla  poprzedzających  punktów 

czasowych 

1

j

w

−

 i 

2

j

w

−

, przy czym 

0

w

 i 

1

w

 wynikają z warunków początkowych 

 

0

2

0

2

0

2

0

3

0

3

0

3

0

1

0

2

0

2

0

(

)

(

)

(

)

( )

( )

( )

(

)

(

)

(

)

n

n

n z

w x

w x

v x

g

w x

w x

v x

g

w

w

w x

w x

v x

g

−

−

−

+

⋅

















+

⋅









=

=

























+

⋅









⋮

⋮

 

(16) 

Wymiar  wektorów 

j

w

  jest  zmniejszony,  ponieważ  warunki  brzegowe  dają  nam 

0,

0

j

w

=

,

1,

0

j

w

=

, 

1,

0

n

j

w

−

=

, 

,

0

n j

w

=

. 

 

P

RZYKŁAD 

2 

Przedstawić równania różnicowe dla siatki prostokątnej opisanej przez związki (1) dla 
 
a)

 

równania Poissona 

 

2

2

2

2

( , )

0

u

u

f x y

x

y

∂

∂

+

+

=

∂

∂

, 

(17) 

 
b)

 

równania Helmholtza 

 

2

2

2

2

2

0

u

u

k u

x

y

∂

∂

+

+

=

∂

∂

, 

(18) 

 

c)

 

równania przewodnictwa ciepła 

 

2

2

2

1

0

u

u

x

a

t

∂

∂

+

=

∂

∂

.   

(19) 

Rozwiązanie: 
a)

 

Zastępując pochodne w równaniach różniczkowych przez odpowiednie ilorazy różnicowe otrzymamy: 

 

(

)

(

) (

)

1,

,

1,

,

1

,

,

1

2

2

1

1

2

2

,

0

i

j

i j

i

j

i j

i j

i j

i

j

u

u

u

u

u

u

f x y

h

g

+

−

+

−

−

+

+

−

+

+

=

. 

(20) 

Jeśli przyjmiemy 

h

g

=

 (siatka kwadratowa) i 

0

f

≡

 (równanie Laplace’a) to uzyskamy wynik: 

 

1,

1,

,

1

,

1

,

4

i

j

i

j

i j

i j

i j

u

u

u

u

u

+

−

+

−

+

+

+

=

. 

(21) 

Wartość funkcji w węźle jest wtedy równa średniej arytmetycznej z wartości funkcji w 4 najbliższych węzłach, 
oddalonych o 

h

. 

 
b)

 

Sprowadzenie do postaci różnicowej daje  

 

(

)

2

2

,

1,

1,

,

1

,

1

4

i j

i

j

i

j

i j

i j

u

g h

u

u

u

u

+

−

−

+

⋅ −

=

+

+

+

. 

(22) 

 
c)

 

Przyjmując 

x

h

∆ =

, 

t

g

∆ =

 otrzymamy 

 

(

)

2

2

2

2

2

1,

,

1,

,

1

2

i

j

i j

i

j

i j

a gu

h

a g u

a gu

h u

+

−

−

−

+

+

= −

. 

(23) 

 

P

RZYKŁAD 

3 

Jednowymiarowy model zginania belki jest opisany równaniem różniczkowym (patrz cz. I, wzór (1.75)) 

 

4

4

d w

p

dx

EJ

=

, 

(24) 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

14 

gdzie 

( )

w x

  jest  funkcją  ugięcia  belki, 

p

  obciążeniem  ciągłym  (N/m), 

E

  modułem  Younga  a 

J

 

odpowiednim momentem bezwładności przekroju. 
Po  zastosowaniu  dyskretyzacji:  określeniu  położenia  węzłów  dla  metody  różnic  skończonych  wg  formuły 

i

x

a ih

= +

,  możemy  zastosować  centralny  iloraz  różnicowy  dla  czwartej  pochodnej  i  równanie  różniczkowe 

(24) zastąpimy przez równanie 

 

2

1

1

2

4

4

6

4

i

i

i

i

i

i

w

w

w

w

w

p

h

EJ

−

−

+

+

−

+

−

+

=

, 

(25) 

gdzie 

( )

i

i

p

p x

=

. 

Możliwe są różne typy warunków brzegowych, które zapisujemy, uwzględniając, że : 
 
a)

 

wielkość ugięcia w punkcie odpowiada wartości funkcji 

( )

w x

, 

b)

 

kąt ugięcia, w zakresie małych deformacji odpowiada wartości  pochodnej 

dw

dx

, 

c)

 

moment gnący w przekroju 

x

 jest związany z funkcją ugięcia wzorem 

 

2

2

( )

d w

M x

EJ

dx

=

, 

(26) 

d)

 

siła tnąca w przekroju 

x

 jest związana z linią ugięcia 

( )

w x

 wzorem 

 

3

3

( )

d w

T x

EJ

dx

=

. 

(27) 

Uwzględniając  powyższe  stwierdzenia  możemy  rozpatrzyć  typowe  warunki  brzegowe  (rys.  7)  zakładając  w 
każdym przypadku, że węzeł 

k

 odpowiada przekrojowi belki z narzuconym warunkiem brzegowym  

 
Przekrój utwierdzony 
Odpowiednie warunki brzegowe mają postać: 

 

(

)

0,

(

)

0,

k

k

dw

w x

x

dx

=

=

 

co po dyskretyzacji prowadzi do zależności 

 

1

1

0,

0.

2

k

k

k

w

w

w

h

+

−

−

=

=

 

(28) 

Jeżeli  punkt 

k

  jest  początkiem  (końcem)  rozpatrywanego  przedziału 

x

  możemy  przyjąć  odpowiednio 

1

0

k

w

+

=

  (

1

0

k

w

−

=

). 

 
Przekrój przegubowo podparty, nieobciążony 
Odpowiednie warunki brzegowe mają postać 

(

)

0

k

w x

=

     oraz     

(

)

0

k

M x

=

. 

Po dyskretyzacji otrzymamy 

 

1

1

2

2

0

i

0

k

k

k

k

w

w

w

w

h

+

−

−

+

=

=

 

(29) 

 

Przekrój swobodny i obciążony momentem 

g

M

 

Warunki brzegowe odpowiadające temu przypadkowi to 

(

)

0

k

T x

=

     i     

(

)

k

g

M x

M

=

. 

Odpowiednie równania różnicowe mają postać: 

 

2

1

1

2

3

2

2

0

k

k

k

k

w

w

w

w

h

+

+

−

−

−

+

−

=

     i     

1

1

2

2

k

k

k

g

w

w

w

EJ

M

h

+

−

−

+

=

. 

(30) 

 

 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

15 

Mg

k - 1

k

k + 1

k + 2

k

k + 1

k + 2

k - 1

k - 2

k

k + 1

k + 2

k - 1

k - 2

 

 

Rys. 7. Warianty warunków brzegowych w modelu MRS belki zginanej (przykład 3) 

 

 

W metodzie różnic skończonych siatka węzłów obejmuje często zarówno węzły leżące 

wewnątrz  analizowanego  obszaru,  jak  i  dodatkowe  węzły  zewnętrzne,  które  ułatwiają 

wprowadzanie warunków brzegowych. 

W  przypadkach  dwu-  i  trójwymiarowych  często  zdarza  się,  że  brzeg  obszaru  przebiega 

między  węzłami  regularnej  siatki  (rys.  8).  W  takim  przypadku  węzłom  leżącym  najbliżej 

konturu przypisujemy wartości wynikające z interpolacji lub ekstrapolacji. 

 

1

2

0

1

2

0

h

h

δ

δ

a)

b)

 

Rys. 8. Interpolacja i ekstrapolacja warunku brzegowego w MRS 

 

 

Jeśli ostatni węzeł siatki (1) leży wewnątrz konturu (rys. 8a), to zamiast narzuconego 

warunku na wartość funkcji 

0

u

u

=

 wprowadzamy związek 

 

0

2

1

hu

u

u

h

δ
δ

+

=

+

. 

(31) 

W przeciwnym przypadku (rys. 8b) mamy 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

16 

 

0

2

1

hu

u

u

h

δ
δ

−

=

−

. 

(32) 

 

Metoda różnic skończonych ma wiele zalet, dzięki którym w przeszłości uważana była 

za  najbardziej  efektywną  metodę  obliczeń  zautomatyzowanych.  Do  najważniejszych  z  nich 

należą prostota algorytmów i bezpośrednie wykorzystanie równań różniczkowych problemu. 

W pewnych obszarach jest do dziś praktycznie stosowana, np. w analizie tarcz i płyt, czy  w 

zagadnieniach przepływowych. 

 

MRS  napotyka  na  szereg  trudności  w  przypadku  obszarów  o  złożonych  kształtach  i 

skomplikowanych  warunków  brzegowych.  Dla  poprawienia  jej  efektywności  opracowano 

pewne  modyfikacje  i  udogodnienia  mające  na  celu  ułatwienie  budowy  modelu  dyskretnego, 

np. algorytmy MRS dla nieregularnych siatek węzłów.  

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

17 

 

3.  METODA ELEMENTÓW BRZEGOWYCH 

 

 

Metoda  elementów  brzegowych  (MEB)  jest  stosunkowo  nową  komputerową  metodą 

rozwiązywania  zagadnień  opisywanych  równaniami  różniczkowymi  cząstkowymi.  Jej 

podstawą  jest  sprowadzenie  zadania  brzegowego  lub  brzegowo-początkowego  opisanego  za 

pomocą  układu  równań  różniczkowych  cząstkowych  z  odpowiednimi  warunkami 

granicznymi do tak zwanych brzegowych równań całkowych. Równania te są rozwiązywane 

w  sposób  przybliżony  przez  podział  brzegu  na  segmenty  zwane  elementami  brzegowymi  i 

przyjęcie  odpowiedniej  aproksymacji  analizowanych  funkcji  na  elementach  brzegowych  za 

pomocą  z  góry  założonych  funkcji  modelujących.  Warto  podkreślić,  że  w  MEB  zwykle 

niezbędna  jest  dyskretyzacja  tylko  brzegu  obszaru  a  nie  całego  analizowanego  obszaru,  co 

prowadzi do znacznego  zmniejszenia liczby niewiadomych w porównaniu do metody różnic 

skończonych  i  metody  elementów  skończonych.  Mówimy,  że  MRS  i  MES  należą  do  tzw. 

metod  obszarowych  przybliżonego  rozwiązywania  zagadnień  brzegowo-początkowych  a 

MEB reprezentuje metody brzegowe. 

Metodę elementów brzegowych omówimy najpierw na przykładzie równania różniczkowego 

Poissona w przypadku dwuwymiarowym. W poprzednim rozdziale przedstawiony był sposób 

rozwiązania  tego  równania  metodą  różnic  skończonych.  Rozwiązanie  zagadnienia  Poissona 

metodą elementów skończonych omówione będzie w rozdziale następnym. 

 

 

3.1.  METODA ELEMENTÓW BRZEGOWYCH DLA RÓWNANIA 

POISSONA 

 

Rozpatrzmy równanie Poissona: 

 

2

2

1

2

2

2

1

2

( ,

)

0

u

u

f x x

x

x

∂

∂

+

+

=

∂

∂

, 

(33) 

gdzie 

1

2

( ,

)

x

x x

=

∈Ω

 

z warunkami brzegowymi 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

18 

 

0

0

( )

,

( )

( )

,

u

q

u x

u

x

u x

q x

q

x

n

=

∈Γ

∂

=

=

∈Γ

∂

 

(34) 

Równanie Poissona i Laplace’a (gdy 

1

2

( ,

)

0

f x x

=

) opisuje wiele zjawisk o dużym znaczeniu 

w  technice:  stacjonarny  przepływ  ciepła,  stacjonarny  bezwirowy  przepływ  nieściśliwej  i 

nielepkiej cieczy, proste  pola magnetyczne i elektryczne [4]. W teorii sprężystości  równanie 

Poissona  opisuje rozkłady naprężeń w przekroju pręta skręcanego .  

Dla  takiego  zagadnienia  przedstawić  można    brzegowe  równanie  całkowe

2

,  które  stanowi 

sformułowanie równoważne problemowi opisanemu przez związki (33) i (34).  

 

( )

( ) ( )

( )

( , )

( )

( , )

( )

( )

( , )

( )

u x

c

u

u x q

x d

x

u

x d

x

f x u

x dR x

n

ξ

ξ

ξ

ξ

ξ

∗

∗

∗

Γ

Γ

Ω

∂

=

Γ

−

Γ

+

∂

∫

∫

∫

 

(35) 

W równaniu (35)  ( )

c

ξ

 jest współczynnikiem liczbowym, którego wartość jest równa 

1

2

 jeśli 

ξ

  leży  na  gładkim  konturze,  a  1  dla 

ξ

  leżącego  wewnątrz  obszaru.  W  przypadku  naroża 

współczynnik ten należy specjalnie obliczać. 

Funkcje  u

∗

  i  q

∗

  zależą  od  położenia  dwóch  punktów:  punktu 

ξ

  zwanego  punktem 

ź

ródłowym i punktu  x  zwanego punktem obserwacyjnym (rys. 9).  

 

x

x

2

1

r

r

r

1

2

Ω

Γ

Γ

u

Γ    Γ =Γ

u

q

q

n

n

n

1

2

ξ

x

 

Rys. 9. Oznaczenia pomocnicze do brzegowego sformułowania zagadnienia Poissona 

 

Funkcja  u

∗

 znana w teorii równań całkowych, jako tak zwane rozwiązanie fundamentalne ma 

postać: 

 

1

1

( , )

ln

2

u

x

r

ξ

π

∗

 

=

=

 

 

, 

(36) 

                                                           

2

 wyprowadzenie równania (35) znaleźć można w monografiach teorii spręzystości 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

19 

gdzie 

2

2

1

1

2

2

(

)

(

)

r

x

x

ξ

ξ

=

−

+

−

. 

Funkcja  q

∗

 określona jest przez pochodną kierunkową   u

∗

 : 

 

( , )

( , )

u

x

q

x

n

ξ

ξ

∗

∗

∂

=

∂

. 

(37) 

Dokonując różniczkowania otrzymamy 

 

1

2

1

2

1

1

2

2

2

,

(

)

,

2

u

u

q

n

n

x

x

r n

r n

q

r

π

∗

∗

∗

∗

∂

∂

=

⋅ +

⋅

∂

∂

− ⋅ + ⋅

=

 

(38) 

gdzie 

,

i

i

i

r

x

ξ

= −

 

1, 2

i

=

,  a 

1

2

,

n

n n

=

  jest  jednostkowym  wektorem  zewnętrznie  normalnym 

do brzegu 

Γ

.  

Wzór (38) można łatwo otrzymać, zauważając, że 

 

i

i

i

i

x

r

r

x

r

r

ξ

−

∂ =

=

∂

. 

(39) 

Brzegowe  równanie  całkowe  (35)  wiąże  ze  sobą  nieznane  funkcje 

( )

u x   i  jej  pochodną 

normalną 

( )

( )

u x

q x

n

∂

=

∂

  na  konturze  obszaru.  W  trzeciej  całce,  obszarowej,  funkcje 

podcałkowe  są  znane.  Równanie  (35)  można  rozwiązać  numerycznie  w  sposób  przybliżony, 

realizując następujący schemat: 

1.

 

Brzeg 

Γ

  dzielimy  na  LE  elementów  brzegowych,  będących  odcinkami  lub 

krzywoliniowymi segmentami z węzłami wyznaczającymi ich kształt. 

2.

 

Na  każdym  elemencie  brzegowym  funkcje 

( )

u x   i 

( )

( )

u x

q x

n

∂

=

∂

  aproksymujemy  za 

pomocą wartości w węzłach i odpowiednich, z góry założonych funkcji interpolacyjnych. 

3.

 

Wykorzystując  przyjętą  dyskretyzację  przekształcamy  równania  całkowe  formułowane 

dla każdego węzła (punktu 

ξ

) w algebraiczne równanie liniowe. 

4.

 

Otrzymujemy  układ  równań  liniowych  wiążący  ze  sobą  wartości 

i

u   i 

i

q   funkcji  ( )

u x   i 

( )

q x   dla  wszystkich  węzłów  konturu.    Układ  ten  rozwiązujemy  po  wprowadzeniu 

warunków  brzegowych  (w  każdym  węźle  znamy  albo 

i

u   albo 

i

q ).  Rozwiązaniem  jest 

wektor nieznanych wartości funkcji  ( )

u x  i jej pochodnych w węzłach brzegu. 

5.

 

Jeśli  chcemy  wyznaczyć  wartości  funkcji  u  dla  punktów  wewnętrznych  wykorzystujemy 

ponownie  wzór  (35)  przenosząc  punkt 

ξ

  do  wewnątrz  obszaru.  Wtedy 

( ) 1

c

ξ

=

  i 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

20 

ponieważ znane są już wartości funkcji i jej pochodnych na konturze zadanie sprowadza 

się do zwykłego numerycznego całkowania. 

 

Rozpatrzmy  na  przykład  przebieg  obliczeń  dla  najprostszych  elementów  brzegowych,  tzw. 

elementów stałych (rys. 10). Element brzegowy jest w tym przypadku odcinkiem a węzeł tego 

elementu usytuowany jest w środku.  

x

x

2

1

P

P

j

i

r

 

 

Rys. 10. Podział brzegu na elementy stałe (o stałych wartościach funkcji aproksymowanych) 

 

Zakładamy,  że  na  każdym  elemencie  funkcje  ( )

u x   i  ( )

q x   przybliżamy  funkcjami  stałymi. 

Otrzymamy wtedy dla każdego węzła i związek 

 

1

1

1

( )

( , ) (

)

( , ) (

)

2

( )

( , )

1, 2,..

j

j

LE

LE

i

i

j

j

i

j

j

j

j

i

u P

u P x q P d

q P x u P d

f x u P x dR

i

LE

∗

∗

=

=

Γ

Γ

∗

Ω

=

Γ −

Γ

+

=

∑

∑

∫

∫

∫

 

(40) 

Dla  obliczenia  wartości  trzeciej  całki,  ze  znanymi  funkcjami  podcałkowymi,  niezbędna  jest 

pomocnicza dyskretyzacja obszaru. Przyjmijmy oznaczenie 

 

( )

( , )

( )

i

i

f

f x u P x d

x

∗

Ω

=

Ω

∫

. 

(41) 

Po numerycznym całkowaniu dla każdego punktu węzłowego otrzymujemy związek 

 

1

1

1

( )

(

)

(

)

,

1, 2...

2

LE

LE

i

ij

j

ij

j

i

j

j

u P

U

q P

Q u P

f

i

LE

∗

∗

=

=

=

⋅

−

⋅

+

=

∑

∑

. 

(42a) 

 

{ }

{ }

{ } { }

1

2

u

U

q

Q

n

f

∗

∗









=

−

+









. 

(42b) 

Otrzymamy  więc  LE  równań  liniowych,  w  których  niewiadome  (42b)  stanowią  nieznane 

wartości  (

)

j

u P   (jeśli  węzeł 

j

P   leży  na  części 

q

Γ

  konturu)  i  ( )

i

q P   (jeśli  węzeł 

i

P   leży  na 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

21 

części 

u

Γ

 konturu). Oznaczając te niewiadome jako 

{ }

y

 układ (42) możemy doprowadzić do 

postaci: 

 

[ ]

{ } { }

A

y

b

=

 

(43) 

gdzie macierz 

[ ]

A

 i wektor prawych 

{ }

b

stron  są znane. 

W rezultacie otrzymujemy kompletną informację o funkcji  ( )

u x  i jej pochodnych na brzegu. 

W  odróżnieniu  do  MRS  i  MES,  gdzie  analogiczna  macierz  jest  zazwyczaj  pasmowa  i 

symetryczna,  w  metodzie  elementów  brzegowych  otrzymujemy  macierz  A  pełną  i 

niesymetryczną.  

 

3.2.  DWUWYMIAROWE  ZAGADNIENIE  TEORII  SPRĘŻYSTOŚCI. 

BRZEGOWE RÓWNANIE CAŁKOWE 

 

 

Opis matematyczny MEB dla zagadnień teorii sprężystości może być formułowany na 

wiele sposobów [5] i jest bardziej złożony niż w metodzie różnic skończonych, czy metodzie 

elementów  skończonych.  W  tym  rozdziale  prezentowany  jest  skrótowy  obraz  tak  zwanego 

podejścia  bezpośredniego,  w  którym  niewiadomymi  w  równaniach  całkowych  są 

przemieszczenia i naprężenia. 

 

Rozpatrujemy izotropowe, jednorodne ciało liniowo sprężyste zajmujące w przestrzeni 

2

R   obszar 

Ω

  i  znajdujące  się  w  stanie  równowagi  pod  działaniem  sił  zewnętrznych: 

powierzchniowych 

1

2

1

2

( ,

)

(

,

)

p x x

p p

=

  i  objętościowych 

1

2

1

2

( ,

)

(

,

)

X x x

X X

=

.  Obciążenia  te 

powodują  powstanie  odpowiedniego  pola  przemieszczeń 

1

2

1

2

( ,

)

( ,

)

u x x

u u

=

  i  pola  naprężeń 

1

2

( ,

),

ij

x x

σ

  ,

1, 2

i j

=

 (rys. 11). 

 

X

p

22

11

12

11

22

x

x

2

1

Γ    Γ =Γ

u

q

Γ

u

Ω

Γ

q

σ

σ

σ

σ

σ

 

Rys. 11. Analizowane ciało sprężyste 

 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

22 

Polu przemieszczeń  u  odpowiada tensor odkształcenia: 

 

(

)

,

,

1

2

ij

i j

j i

u

u

ε

=

+

. 

(44) 

Różniczkowe równanie równowagi  w układzie kartezjańskim 

1 2

x x

 ma postać: 

 

,

0

ji j

i

X

σ

+

=

, 

(45) 

gdzie 

ij

σ

 jest symetrycznym tensorem naprężenia. 

Ciało  jest  podparte  na  części 

u

Γ

  brzegu  (znane  przemieszczenia 

1

2

ˆ ˆ

( ,

)

u u u

∆

)  i  obciążone  na 

Części 

p

Γ

 (znane obciążenia 

1

2

ˆ

ˆ

(

,

)

p

p p

∆

=

). Mamy więc warunki brzegowe:  

 

ˆ ,

(

)

i

ji

j

i

p

p

n

p

x

σ

=

=

∈Γ

, 

(46) 

gdzie 

1

2

( ,

)

n

n n

=

 jest wersorem normalnym zewnętrznie do konturu oraz 

 

ˆ ,

(

)

i

i

u

u

u

x

=

∈Γ

. 

(47) 

Związek  między  tensorami  odkształcenia  i  naprężenia  określają  poniższe  zależności, 

wyrażające uogólnione prawo Hooke’a. Dla płaskiego stanu naprężenia (PSN): 

 

2

1

1

ij

ij

ij

kk

E

Ev

v

v

σ

ε

δ ε

=

+

+

−

, 

(48) 

a dla płaskiego stanu odkształcenia (PSO): 

 

1

(1

)(1 2 )

ij

ij

ij

kk

E

vE

v

v

v

σ

ε

δ ε

=

+

+

+

−

, 

(49) 

gdzie E jest modułem Younga  a v – stałą Poissona. 

 

Zauważyć  można,  że  związki  konstytutywne  (48)  i  (49)  różnią  się  wyłącznie 

współczynnikami  określonymi  przez  stałe  materiałowe.  Przez  podstawienie  do  (48) 

zmodyfikowanych wartości stałych: 

 

2

,

1

1

E

v

E

v

v

v

′

′

=

=

−

−

, 

(50) 

otrzymamy związek (49), a przez podstawienie do (49) 

 

2

(1 2 )

,

(1

)

1

E

v

v

E

v

v

v

+

′

′

=

=

+

+

, 

(51) 

otrzymamy związek (48). 

 

W  przypadku  wykorzystania  w  prawie  Hooke’a  modułu  sprężystości  postaciowej 

2(1

)

E

G

v

=

+

 i stałej Poissona, podstawienia (50) i (51) przejdą odpowiednio w (52) i (53): 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

23 

 

,

1

v

v

G

G

v

′

′

=

=

−

, 

(52) 

 

,

1

v

v

G

G

v

′

′

=

=

+

. 

(53) 

Powyższe związki pozwalają na proste dostosowanie algorytmów numerycznych do realizacji 

obliczeń  zarówno  dla  PSN  jak  i  PSO  w  zależności  od  określającego  stan  parametru 

sterującego. Dalsze zależności przedstawiane będą dla płaskiego stanu odkształcenia. 

Podstawiając  (44)  do  (49)  a  następnie  do  (45)  otrzymujemy  przemieszczeniowe 

równania różniczkowe Naviera: 

 

,

,

0

1 2

j kk

k kj

j

G

Gu

u

X

v

+

+

=

−

. 

(54) 

 

Punktem  wyjścia  dla  bezpośredniej  metody  elementów  brzegowych  są  wzory 

Somigliany.  Wzory  Somigliany  wyprowadzić  można  ze  znanego  twierdzenia  Bettiego  o 

wzajemności prac: 

 

i

i

i

i

i

i

i

i

X u dV

p u dA

X u dV

p u dA

Ω

Γ

Ω

Γ

′

′

′

′

+

=

+

∫

∫

∫

∫

. 

(55) 

Twierdzenie to można wyrazić następująco: 

 

Jeśli na ciało działają kolejno dwa układy obciążeń to praca pierwszego układu obciążeń 

( , )

X p   na  przemieszczeniach  wywołanych  przez  drugi  układ  ( )

u

′

  jest  równa  pracy  drugiego 

układu obciążeń  (

,

)

X p

′ ′

 na przemieszczeniach wywołanych przez układ pierwszy  ( )

u . 

x

r

ξ

X =   (x

ξ

i

i k

(k=1)

P   (x

ξ

i

δ

x

2

x

1

δ

   -   )

(k)

,   )

 

Rys. 12. Pomocniczy stan równowagi: siła w przestrzeni nieograniczonej 

 

 

Aby  otrzymać  wzory  Somigliany  przyjmijmy,  że  dany  jest  układ  obciążeń  ciała  p , 

X .  Drugi  stan  równowagi  konstruujemy  wyodrębniając  w  nieograniczonej  przestrzeni 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

24 

sprężystej 

∞

Ω

  obszar 

Ω

  odpowiadający  obszarowi  zajmowanemu  przez  analizowane  ciało 

(rys.  12).  W  punkcie 

ξ

∈Ω

  przykładamy  jednostkową,  objętościową  siłę  skupioną  X

′

  o 

kierunku k: 

 

(

)

i

ik

X

x

δ

ξ δ

′ =

−

. 

(56) 

W powyższym wzorze  (

)

x

δ

ξ

−

 jest funkcją (dystrybucją) Diraca o własnościach: 

 

( )

0,

0,

( )

,

0,

( ) (

)

( )

( ).

x

x

x

x

f x

x

d

x

f

σ

δ

δ

ξ

ξ

Ω

=

≠

→ ∞

=

−

Ω

=

∫

 

(57) 

Wynika stąd, że 

i

X

′

 określone wzorem (56) jest wszędzie równe zeru poza punktem  x

ξ

=

, 

gdzie  ma  jedną  składową  (k)  o  wartości  nieskończonej.  Siłę  określoną  przez  funkcję 

i

X

′

 

możemy traktować jako siłę jednostkową, gdyż: 

 

( )

i

ik

X dV x

δ

Ω

′

=

∫

. 

(58) 

Rozwiązanie  układu  równań  przemieszczeniowych  (54)  dla  takiej  siły,  z  warunkiem 

( )

0

u x

→

  dla 

x

→ ∞

,  daje  w  rezultacie  w  przypadku  płaskiego  stanu  odkształceń  pole 

tensorowe przemieszczeń [2]: 

 

( )

( )

1

, ,

2

( , )

ln

k

i

i k

ik

U

x

C r r

C

r

ξ

δ





=

−



 , 

(59) 

gdzie: 

 

1

1

8

(1

)

C

G

v

π

=

−

, 

2

3 4

C

v

= −

, 

 

(

)

1/ 2

(

)(

)

i

i

i

i

r

x

x

ξ

ξ

=

−

−

, 

,

i

i

i

r

r

r

x

r

∂

=

=

∂

, 

i

i

i

r

x

ξ

= −

. 

Funkcje 

( )

k

i

U

  nazywane  są  rozwiązaniami  (funkcjami)  fundamentalnymi  Kelvina  równań 

elastostatyki lub funkcjami przemieszczeniowymi Greena dla przestrzeni nieograniczonej [4, 

16]. 

 

Znajomość  tensora  przemieszczeniowego  Greena  (59)  pozwala  na  określenie  pola 

naprężeń 

( )

( )

,

k

ji

x

σ

ξ

 i następnie obciążeń 

( )

k

i

P

 odpowiadających brzegowi 

Γ

. Korzystając ze 

związku (44) mamy: 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

25 

 

(

)

(

)

( )

( )

( )

,

,

( )

1

2

2

2

1

,

2

2

4

(

1)

.

2

k

k

k

ij

i j

j i

i j k

k

ij

ij k

jk i

ik j

U

U

r r r

C

r

C

r

r

r

r

ε

ε

δ

δ

δ

=

+





=

−

−

−

+









 

(60) 

Wyznaczony z prawa Hooke’a wzór określający składowe tensora naprężeń ma postać: 

 

(

)

( )

3

4

2

2

2

i j k

k

ij

ij k

ik i

ik j

r r r

C

C

r

r

r

r

r

σ

δ

δ

δ





=

−

−

−









, 

(61) 

gdzie: 

3

1

4 (1

)

C

v

π

=

−

, 

4

1 2

C

v

= −

. 

Mając  dany  tensor  naprężeń  odpowiadający  przyłożeniu  w  punkcie 

1

2

( ,

)

ξ

ξ ξ

=

  siły 

jednostkowej  X

′

 można wyznaczyć obciążenia linii konturu 

Γ

: 

 

( )

( )

( )

3

4

4

2

2

( , )

( , )

( ),

,

2

(

)

.

k

k

i

ij

j

k

i k

i

i k

k i

ik

l

l

P

x

x

n x

x

A

C

r r

P

C n r

n r

C

r n

r

r

ξ

σ

ξ

δ

=

∈









=

−

−

+

















 

(62) 

Wstawiając wyrażenia (59) i (62) do (55) otrzymujemy wzory Somigliany: 

 

(

)

( )

( )

( )

( )

( )

( , )

( , ) ( )

( )

( , )

k

k

k

k

i

i

i

i

i

i

u

p x U

x

P

x

u x d

X x U

x

d

ξ

ξ

ξ

ξ

Γ

Ω

=

−

Γ +

Ω

∫

∫

 

(63) 

Ze wzorów Somigliany wynika, że znając rozkład sił masowych 

( )

i

X x  oraz przemieszczenia 

( )

i

u x  i obciążenia 

( )

i

p x  brzegu ciała 

Γ

 możemy obliczyć przemieszczenia 

k

u  w dowolnym 

punkcie wewnętrznym tego ciała

3

. 

 

Praktyczne  zastosowanie  wzorów  Somigliany  stało  się  możliwe  przez  przeniesienie 

punktu 

ξ

  (leżącego  wewnątrz  obszaru  ciała  i  w  którym  określane  jest  przemieszczenie)  na 

brzeg ciała 

Γ

. Można wykazać, że dla punktu 

ξ

∈Γ

 otrzymujemy  

 

(

)

( )

( )

( )

( )

( )

( )

( , )

( )

( , )

( )

( )

( , )

( ),

k

k

kj

j

i

i

i

i

k

i

i

c

u

p x U

x

u x P

x

d

x

X x U

x

d

x

ξ

ξ

ξ

ξ

ξ

Γ

Ω

=

−

Γ

+

+

Ω

∫

∫

 

(64) 

gdzie 

kj

c  zależy od kształtu konturu w punkcie 

ξ

 i  od współczynnika Poissona. 

Gdy brzeg w punkcie 

ξ

 jest gładki, lewa strona równania (64) ma postać 

1

( )

2

k

u

ξ

. 

                                                           

3

 Konieczność znajomości zarówno przemieszczeń 

( )

i

u x  jak i obciążeń 

( )

i

p x  brzegu dla zastosowania zasady 

Somigliany  powodowała,  że  przed rozwojem  metody  elementów  brzegowych  uważano,  że  „wzory  Somigliany 
mają  jedynie  teoretyczne  znaczenie,  bowiem  na  brzegu  ciała  mogą  być  znane  albo  przemieszczenia  albo 
obciążenia, ale każdy z tych typów warunków brzegowych oddzielnie” (Nowacki W., Teoria sprężystości) 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

26 

 

W dalszych rozważaniach przyjmiemy 

0

i

X

≡

. Siły masowe mogą być uwzględnione 

w  MEB.  Jednak  w  takim  przypadku  wymagana  jest  pomocnicza  dyskretyzacja  obszaru  by 

obliczyć wartość trzeciej całki w równaniu (64). Przyjęcie odpowiedniej dyskretyzacji brzegu 

i  aproksymacji  warunków  brzegowych  sprowadza  zadanie  rozwiązania  równania  (64)  do 

układu  równań  liniowych  z  niewiadomymi  określającymi  nieznaną  część  warunków 

brzegowych.  Po  wykonaniu  tego  zadania  można  obliczać  przemieszczenia,  odkształcenia  i 

naprężenia  w  wybranych  punktach  wewnętrznych  obszaru.  Przemieszczenia  obliczane  są 

bezpośrednio ze zdyskretyzowanych wzorów Somigliany (63). 

 

Wyrażenia  na  składowe  tensorów  odkształcenia  i  naprężenia  otrzymamy  po 

zróżniczkowaniu  równania  względem 

i

ξ

  i  wykorzystaniu  wzorów  (44)  i  (49).  Otrzymujemy 

([2,4]): 

 

( )

( , )

( )

( )

( , )

( )

( )

ij

kij

k

kij

k

A

x

p x d

x

B

x

u x d

x

ε ξ

ξ

ξ

Γ

Γ

=

Γ

−

Γ

∫

∫

, 

(65) 

gdzie: 

 

(

)

( )

( )

1

, ,

,

2

,

,

1

( , )

,

2

2

4

(1

)(

) ,

2

i

j

k

k

kij

j

i

kij

ij jk

i

j

k

jk

i

ik

j

U

U

A

x

C

A

r

r r r

C

r

r

r

ξ

ξ

ξ

δ

δ

δ







∂

∂

=

+













∂

∂









= −

−

+ −

+



 

(66) 

 

(

)

(

) (

)

}

( )

( )

3

,

,

,

, ,

,

2

4

, ,

,

,

,

,

,

1

( , )

,

2

2

4

2

2

,

i

j

k

k

kij

j

i

kij

ij

k

ik

j

jk

i

i

j

k

ik

j

jk

i

ij

k

i

j

k

k

j

i

k

i

j

P

P

B

x

C

r

B

r

v

r

r

r r r

r

n

C

n

n

n

r r n

v r r n

r r n

ξ

ξ

ξ

δ

δ

δ

δ

δ

δ







∂

∂

=

+











∂

∂









∂









=

+

+

−

+









∂







+

+

−

+

+

+





 

(67) 

oraz 

 

( )

( , )

( )

( )

( , )

( )

( )

ij

kij

k

kij

k

A

A

D

x

p x dA x

S

x

u x dA x

σ ξ

ξ

ξ

=

−

∫

∫

, 

(68) 

gdzie 

 

(

)

3

4

,

,

,

, ,

,

2

kij

ik

j

jk

i

ij k

i

j

k

C

D

C

r

r

r

r r r

r

δ

δ

δ





=

+

−

+



 , 

(69) 

 

(

)

(

)

(

)

(

)

(

)

3

,

,

,

, ,

,

2

, ,

,

,

,

,

2

2

1 2 )

(

4

1 2

2

4

1

2

.

kijh

ij

k

ik

jk

jk

i

i

j

k

i

j

k

ik

j

jk

i

ij

k

k

j

i

k

i

j

C G

r

S

v

r v

r

r

r r r

r

n

v

r r n

n

n

v

n

v r r n

r r n

δ

δ

δ

δ

δ

δ

∂









=

−

+

−

+









∂





+ −

+

+

+

−

+

+

+

 

(70) 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

27 

We wzorach (3.21)÷(3.25)  

 

,

i

i

i

i

i

i

x

r

r

r

r

x

r

r

ξ

ξ

−

∂

∂

=

= −

=

=

∂

∂

. 

 

 

3.3.  ALGORYTMY NUMERYCZNE MEB W DWUWYMIAROWYM 

ZADANIU TEORII SPRĘŻYSTOŚCI 

 

 

Najważniejsze  elementy  postępowania,  jakie  możemy  wyróżnić  w  analizie  metodą 

elementów brzegowych to:  

a)

 

Dyskretyzacja  brzegu  na  elementy  brzegowe.  Przyjmijmy  oznaczenia:  LE  –  liczba 

elementów  brzegowych,  LW  –  liczba  węzłów  konturu,  LWE  –  liczba  węzłów  elementu 

brzegowego. 

b)

 

Przyjęcie  funkcji  modelujących  przebieg  przemieszczeń  ( )

i

u   i  obciążeń  (

)

i

p   na 

elemencie w zależności od ich wartości w punktach węzłowych, 

c)

 

Przedstawienie  związków  (64)  dla  wszystkich  punktów  węzłowych  jako  układu  2· LW 

równań liniowych typu 

[ ]{ } [ ]{ }

P u

U

p

=

 

z  nieznanymi  2· LW  składowymi  przemieszczeń  i  obciążeń  węzłowych  pozostałymi  po 

uwzględnieniu znanych warunków brzegowych, 

d)

 

Rozwiązanie przekształconego układu równań 

[ ]{ } { }

A x

b

=

 

a więc określenie nieznanych obciążeń i przemieszczeń w węzłach brzegowych, 

e)

 

Obliczenie  składowych  tensora  naprężeń  na  brzegu,  odpowiadających  znanym  już 

obciążeniom i przemieszczeniom, 

f)

 

Obliczenie  przemieszczeń  i  naprężeń  w  wybranych  punktach  wewnętrznych  ze 

zdyskretyzowanych wzorów (65) i (68).  

Aproksymacja kształtu elementu, przemieszczeń i obciążeń określone mogą być wzorami 

 

1

1

1

( )

( ) ( ),

( )

( ) ( ),

( )

( )

( ),

LWE

G

i

l

i

l

LWE

u

i

l

i

l

LWE

p

i

l

i

l

x

x l

u

u l

p

p l

η

ψ η

η

ψ η

η

ψ η

=

=

=

=

=

=

∑

∑

∑

 

(71) 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

28 

gdzie LWE – jest liczbą węzłów elementu, 

η

 – lokalną współrzędną określającą położenie na 

linii  elementu  (

1,1 )

η

∈< −

>

, 

( )

i

x l , 

( )

i

u l , 

( )

i

p l   oznaczają  współrzędne,  przemieszczenia  i 

obciążenia  węzła  l  a 

ψ

∗

∗

  są  założonymi  z  góry  funkcjami  modelującymi.  Zazwyczaj 

przyjmuje się ten sam sposób aproksymacji wszystkich przybliżanych funkcji:  

 

,

1, 2.

G

u

P

l

l

l

l

l

ψ

ψ

ψ

ψ

=

=

=

=

 

 

η = 1

1

2

η = −1

η 

n

Ψ  = 1

2

Ψ  = 1

1

Ψ  =    ( 1−η)
Ψ  =    ( 1+η)

η        −1,1

1
2

1

2

1
2

 

Rys. 13. Liniowy element brzegowy 

 

Rysunek  13  ilustruje  przypadek,  gdy  funkcje 

l

ψ

  są  liniowe.  W  takim  przypadku  równanie 

(64) dla wybranego węzła 

j

M  konturu przybiera postać 

 

( )

( )

1

(

) (

)

(

( )

( ,

)

( ,

)

( )

l

LE

k

k

ki

j

i

j

i

i

j

i

i

j

l

l

L

c

M u M

p x U

x M

u P

x M dL x

=





=

−













∑ ∫

, 

(72) 

gdzie LE jest liczbą elementów konturu a 

l

L  oznacza l-ty element konturu (rys. 13). 

Po uwzględnieniu związków (71) i odpowiednich przekształceniach mamy: 

 

(

)

(

)

(

)

}

1

1

2

1

1

1

1

1

2

1

1

(

)

( )

(

)

( )

(

)

( ),

( )

( ) (

)

( ) (

)

( ( ),

)

( )

LE

k

ki

i

j

i

l

i

l

i

j

l

k

i

l

i

l

i

j

c u M

p M

p M

U

x

M

J

d

u M

u M

P x

M

J

d

ψ η

ψ η

η

η η

ψ η

ψ η

η

η η

+

=

+

−



=

+

+





−

+

∑ ∫

∫

 

(73) 

gdzie 

( )

J

η

 jest jakobianem przejścia z układu 

1 2

x x

 na zmienną 

η

. Dla przypadku elementu 

liniowego jest 

( )

2

j

l

J

η

=

, gdzie 

j

l  oznacza długość elementu j. 

Po dalszych przekształceniach równanie (73) można doprowadzić do postaci: 

 

(

)

( )

(

)

( )

( )

( )

1

1

(

)

,

,

LW

LW

k

k

ki

i

j

i

n

j

i

n

i

n

j

i

n

n

n

c u M

U

M M

p M

P

M M

u M

=

=

=

∆

−

∆

∑

∑

 

(74) 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

29 

gdzie  LW  jest  liczbą  węzłów  konturu  a  współczynniki 

(

)

( )

,

k

i

n

j

U

M M

∆

  i 

(

)

( )

,

k

i

n

j

P

M M

∆

  są 

wynikiem  odpowiednich  całkowań  w  związku  (73).  Całkowania  te  przeprowadzane  są 

numerycznie.  Dla 

j

n

=

  całki  składające  się  na 

(

)

(

)

( )

( )

,

,

k

k

i

n

j

i

j

j

U

M M

U

M M

∆

= ∆

 

i zawierające  osobliwości  funkcji  podcałkowych  obliczane  są  ze  wzorów  analitycznych  a 

całki  składające  się  na 

(

)

(

)

( )

( )

,

,

k

k

i

n

j

i

j

j

P

M M

P

M M

∆

= ∆

  obliczane  są  łącznie  ze 

współczynnikami 

ki

c  w dalszym etapie z warunku ruchu ciała jako bryły sztywnej. 

Zapisując związek (74) w postaci macierzowej otrzymamy 

 

(1)

(1)

(1)

1

1

2

1

1

11

(2)

(2)

(2)

1

1

2

1

1

21

(

,

)

(

,

)

(

,

)

(

)

(

,

)

(

,

)

(

,

)

(

)

j

j

j

j

j

j

j

j

j

j

P

M M

P

M M

P

M M

c

M

P

M M

P

M M

P

M M

c

M



∆

∆

∆

+



∆

∆

∆

+



……

……

 

 

1

1

2

1

(1)

(1)

(1)

1

2

12

1

2

(2)

(2)

(2)

2

2

22

1

2

1

2

(

)

(

)

(

)

(

,

)

(

)

(

,

)

(

,

)

(

)

(

,

)

(

)

(

,

)

(

,

)

(

)

(

)

j

j

j

j

LW

j

LW

j

j

j

j

j

LW

j

LW

j

LW

LW

u M

u M

u M

P

M M

c

M

P

M

M

P

M

M

u M

P

M M

c

M

P

M

M

P

M

M

u M

u M























∆

+

∆

∆





=

 



∆

+

∆

∆

 

























⋮

……

……

⋮

 

 

(1)

(1)

(1)

(1)

1

1

2

1

1

1

2

(2)

(2)

(2)

(2)

1

1

2

1

1

1

2

(

,

)

(

,

)

(

,

)

(

,

)

(

,

)

(

,

)

(

,

)

(

,

)

j

j

j

j

j

j

j

j

j

j

U

M M

U

M M

U

M M

U

M M

U

M M

U

M M

U

M M

U

M M



∆

∆

∆

∆



∆

∆

∆

∆



……

…

……

…

 

 

1

1

2

1

(1)

(1)

1

1

2

(2)

(2)

2

1

2

1

2

(

)

(

)

(

)

(

,

)

(

,

)

(

)

(

,

)

(

,

)

(

)

(

)

j

LW

j

LW

j

j

LW

j

LW

j

LW

LW

p M

p M

p M

U

M

M

U

M

M

p M

U

M

M

U

M

M

p M

p M























∆

∆





 



∆

∆

 

























⋮

…

…

⋮

 

(75) 

Związek ten stanowi układ dwóch równań z  2 LW

⋅

 niewiadomymi sformułowany dla węzła 

j

M   (rys.  14).  Niewiadomymi  jest  część  elementów  wektora 

{ }

u

  i  część  elementów 

wektora

{ }

p

.  Jeśli  warunek  brzegowy  dotyczy  przemieszczenia  znane  jest 

i

u   a  niewiadomą 

stanowi 

i

p  a jeśli znane jest obciążenie 

i

p  do wyznaczenia zostaje 

i

u . 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

30 

x

x

2

1

M (  )

M

M

M

x

n

η

j

l-1

l

l+1

L

i

r (   ,  )

x

ξ

ξ

 

Rys. 14. Sposób przeprowadzania całkowania po konturze 

 

Formułując analogiczne związki dla pozostałych węzłów otrzymujemy układ  2 LW

⋅

 równań 

z  2 LW

⋅

 niewiadomymi: 

 

[ ]

{ }

[ ]

{ }

P u

U

p

=

. 

(76) 

Uwzględniając,  że  warunki  brzegowe  zadania  określają  nam  2 LW

⋅

  wartości  składowych 

wektorów 

{ }

u

 i 

{ }

p

 równanie (76) można przekształcić do postaci: 

 

[ ]

{ } { }

A

y

b

=

, 

(77) 

gdzie 

{ }

y

 zawiera nieznane wartości węzłowe przemieszczeń i obciążeń. 

 

Dla modelowania nieciągłości obciążeń na brzegu wprowadza się tzw. węzły podwójne 

(rys.  15).  Węzeł  podwójny  stanowią  dwa  węzły  o  identycznych  współrzędnych,  w  których 

niezależnie określane są obciążenia 

i

p  a przemieszczenia z założenia są identyczne. 

k+3

k+2

k+1,k

k-1

P

 

Rys. 15. Węzeł podwójny (k+1,k) na konturze analizowanego obszaru 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

31 

 

Określenie,  po  rozwiązaniu  układu  równań  (77),  wszystkich  składowych  obciążeń  i 

przemieszczeń  węzłowych  pozwala  na  łatwe  wyznaczenie  stanu  naprężeń  w  węzłach 

brzegowych.  Naprężenia  obliczane  są  w  układzie  lokalnym 

1

2

,

x x

′ ′

  związanym  z  elementem 

brzegowym (rys. 16). Tak na przykład dla węzła n mamy 

 

11

1

2

12

1

2

( ) cos

( ) sin ,

( ) sin

( ) cos .

p n

p n

p n

p n

σ

α

α

σ

α

α

′ =

+

′ = −

+

 

(78) 

Naprężenie  normalne  w  kierunku  stycznym  do  brzegu  obliczamy  wykorzystując  definicje 

odkształceń i prawo Hooke’a (44) i (49): 

 

2

22

11

2

2

1

1

E

u

v

v

x

v

σ

σ

∂

′

′

=

+

′

−

∂

−

. 

(79) 

Uwzględniając, że przemieszczenie jest zmienne liniowo wzdłuż elementu możemy zapisać 

 

2

2

22

11

2

(

1)

( )

( )

( )

1

1

1

E

u n

u n

v

n

n

v

v

σ

σ

′

′

+ −

′

′

=

+

−

+

. 

(80) 

Otrzymane  składowe  stanu  naprężenia  transformowane  są  następnie  do  układu  globalnego 

1

2

,

x x

. 

x

x

2

1

α

22

σ

22

σ

11

σ

12

σ

'

'

'

'

11

σ

'

x

2

x

1

'

'

n+2

n+1

n-1

 

 

Rys. 16. Wyznaczanie naprężeń w węzłach konturu 

 

 

Dla każdego węzła naprężenia liczone są dwukrotnie przy wykorzystaniu wyników z 

obu  elementów  zawierających  dany  węzeł.  Końcowe  naprężenia  w  węźle  stanowią  średnią 

arytmetyczną z rezultatów otrzymanych dla obu elementów. 

 

Określenie  naprężeń  brzegowych  kończy  zwykle  proces  obliczeń.  W  dowolnym, 

wybranym  punkcie  wewnętrznym  można  jednak  określić  przemieszczenie  i  naprężenia 

bezpośrednio  za  pomocą  wzorów  (63)  i  (68).  Wzory  te  są  dyskretyzowane  zgodnie  z 

opisanymi  powyżej  zasadami.  Odpowiedni  algorytm  automatyzujący  określenie  siatki 

punktów  wewnętrznych  może  umożliwiać  otrzymywanie  rozwiązania  dla  całego 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

32 

analizowanego  obszaru.  Schemat  działania  typowego  programu  metody  elementów 

brzegowych dla dwuwymiarowych zadań teorii sprężystości przedstawia rys. 17. 

 

1.  Wprowadzenie informacji o geometrii 

konturu, warunkach brzegowych i stałych 
materiałowych. 

 

2.  Interakcyjna prezentacja graficzna danych. 

 

3.  Budowa układu równań. 

 

4.  Rozwiązanie układu równań. 

 

5.  Obliczanie naprężeń brzegowych. 

 

6.  Interakcyjna prezentacja graficzna wyników. 

 

7.  Obliczenia dla wybranych punktów 

wewnętrznych lub dla automatycznie 
wygenerowanej siatki punktów 
wewnętrznych. 

 

8.  Prezentacja izolinii wybranych wielkości. 

 

Rys. 17. Schemat podstawowych działań w programie MEB dla dwuwymiarowych zagadnień 

teorii sprężystości 

 
P

RZYKŁAD 

4

    

Rura grubościenna

 

Analizie poddano zagadnienie płaskiego stanu odkształceń przedstawione na rys. 18.  

a = 1  10   m

b = 2.5  10   m

E = 2  10   MPa

    = 0.3

p  = 100 MPa

-2

5

-2

o

a

b

C

D

E

F

A

p

B

o

ν

x

2

x

1

 

 

Rys. 18. Rura grubościenna poddana ciśnieniu wewnętrznemu (przekrój) 

 
Obliczenia  przeprowadzono  dla  fragmentu  ABCD  obszaru  przyjmując  symetryczne  warunki  brzegowe  na  odcinkach  AB 

1

(

0)

u

=

 i CD 

2

(

0)

u

=

. Przyjęto trzy stopnie dyskretyzacji w modelu numerycznym MEB: 18 elementów brzegowych 

(A), 36 elementów brzegowych (B) i 72 elementy brzegowe (C). Wyniki obliczeń komputerowych (naprężenia zredukowane 

G

RZEGORZ 

K

RZESIŃSKI

.

  

MES_1.

  

C

ZĘŚĆ 

1.

  MATERIAŁY DO WYKŁADU

.

  

 

 

33 

Hubera-Misesa 

red

σ

 i przemieszczenia promieniowe 

r

u

 w punktach E i F) porównano z wynikami rozwiązania ścisłego, w 

którym mamy: 

 

(

)

(

)

2

2

0

2

2

( )

1 2

1

r

p

ab

r

b

u

r

v

v

v

b

a

E

b

r





=

=

−

−

+ +





−





, 

 

1

4

2

2

0

2

2

( )

3

1 4 (1

)

red

a

b

r

p

v

v

b

a

r

σ





 

=

+ −

−





 





−

 





. 

Otrzymane rezultaty prezentuje tablica 1. 

Tablica 1. 

Przykład 4 – wyniki obliczeń 

 

( )

red

E

σ

 

[MPa] 

( )

red

F

σ

 

[MPa] 

( )

r

u E  

4

[10

m]

−

 

( )

r

u F  

4

[10

m]

−

 

A 

211,6 

32,7 

0,845 

0,453 

B 

209,0 

33,6 

0,831 

0,441 

C 

207,2 

33,8 

0,826 

0,435 

Wartości 

dokładne 

206,3 

33,9 

0,823 

0,443 

Rysunek 19 prezentuje przyjęty podział brzegu i obraz rozkładu naprężeń zredukowanych na konturze dla przykładu C. 

 

 

Rys. 19. Rura grubościenna pod ciśnieniem wewnętrznym – dyskretyzacja i rozkład 

naprężenia zredukowanego Hubera-Misesa na konturze analizowanego obszaru 

 

 

Metoda  elementów  brzegowych  rozwijała  się  najbardziej  intensywnie  w  latach 

siedemdziesiątych i osiemdziesiątych ubiegłego wieku. Duże nadzieje wiązano z nią głównie 

z  powodu  możliwości  ograniczenia  dyskretyzacji  jedynie  do  brzegu  (konturu)  ciała.  Dzięki 

temu w porównaniu do MRS i MES dyskretny model obliczeniowy MEB opisany może być 

przez znacznie mniejszą liczbę niewiadomych. Tak na przykład podział sześcianu określony 

przez siatkę n 

x

 n 

x

 n=n

3

 węzłów w metodach obszarowych (MRS i MES) odpowiada liczbie 

2

6

12

8

n

n

−

+

  węzłów  w  metodzie  elementów  brzegowych.  Jeśli  przyjmiemy 

40

n

=

  mamy 

odpowiednio 64000 węzłów w porównaniu do 9128 . Korzyści z dyskretyzacji jedynie brzegu 

są jednak znacznie mniejsze dla obszarów wielospójnych lub obszarów o złożonym kształcie.