Arystoteles i modelowanie numeryczne

ArystotelesW dniach 9-11 listopada 2015 uczestniczyłem w szkoleniu w Salonikach na Uniwersytecie Arystotelesa zorganizowanym przez European Cooperation of Science and Technology. Spędziłem kilka słonecznych dni w śródziemnomorskim klimacie korzystając ze słonecznej pogody, chłonąc wiedzę wprost z Uniw. Arystotelesa i racząc się lokalnymi przysmakami.

Szkolenie było związane z oprogramowaniem gprMax do modelowania georadarowego. Program służy do symulacji rozchodzenia się fal georadarowych w ośrodku gruntowym. W ostatnim czasie cały kod źródłowy został przepisany na nowy język programowania Python i obecnie tylko poruszamy się w nowym środowisku programowania. Nie miałem do czynienia z poprzednimi wersjami więc nie jestem w stanie napisać wad i zalet pomiędzy wersjami.

Ekipa ze szkolenia
Międzynarodowa grupa biorąca udział w szkoleniu

Prowadzącym szkolenie był sam autor Dr Antonis Giannopoulos absolwent Uniwersytetu Arystotelesa w Salonikach, a od połowy lat 90 pracownik Uniwersytetu w Edynburgu. Program jest zupełnie darmowy i został napisany w latach 90 i do tej pory jest sukcesywnie rozwijany. Drugą osobą prowadzącą szkolenie był Dr Craig Warren, który w rozprawie doktorskiej zajmował się modelowaniem numerycznym wysokoczęstotliwościowych anten georadarowych, uzyskał także grant i kontrakt na współpracę z Departamentem Obrony w celu rozwijania technologii wykrywania niemetalicznych min lądowych. Szkolenie trwało 3 dni i można je podzielić na 2 części teoretyczną podstawy GPR, metody numeryczne i podstawy teorii Maxwella itp. Druga część praktyczna, w której zaczęliśmy przygodę z programem. Nie będę rozpisywać się o podstawach metody georadarowej oraz teoretycznych aspektach modelowań numerycznych tylko przejdę od razu do zastosowania praktycznego i przykładów.

GprMax niestety nie posiada GUI wszystkie operacje przeprowadzamy w terminalu, w MS Win w oknie poleceń, analogicznie w systemach Linux. W systemie musi być zainstalowany (najlepiej) najnowszy Python 3.5 i biblioteki wskazane przez autora na stronie internetowej. Modele budujemy w formie tekstowej czyli każdy element musimy opisać komendą i parametrami w kodzie. Każdą operację musimy wywoływać z okna poleceń tak jak wykonanie modelowania, podgląd geometrii, plot wyników itp. Jest to stosunkowo uciążliwe ale można zapamiętać wszystkie komendy i przyzwyczaić się.

W programie możemy modelować w 3 przestrzeniach 1D, 2D, 3D. W nomenklaturze anglojęzycznej mówimy A scan - czyli 1D, jeden impuls EM w jednej pozycji  i zarejestrowana odpowiedź ośrodka. B scan (2D) kolekcja wielu A scan np wzdłuż jakiegoś profilu z równym interwałem odległości pomiędzy trasami 1D. C scan analogicznie złożenie wielu B scan (profili 2D) także ortogonalnych w Polskiej nomenklaturze mówimy o pseudo 3D lub 2.5 D ponieważ nie jest to faktyczne 3D tak jak technologie stosowane w badaniach sejsmicznych 3D.

Poniżej prezentuję plik wsadowy modelowania 1D. W krótkim opisie jest to ośrodek gruntowy ograniczony od góry powierzchnią ziemi i powietrzem z obiektem cylindrycznym w środku (np. rura).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
#title: A-scan metalowego cylindra położony w dielektrycznej półprzestrzeni
#domain: 0.240 0.190 0.002
#dx_dy_dz: 0.002 0.002 0.002
#time_window: 3e-9
#time_step_limit_type: 2D
#pml_cells: 10 10 0 10 10 0

#material: 6 0 1 0 half_space

#waveform: ricker 1 1.5e9 my_ricker
#hertzian_dipole: z 0.100 0.170 0 my_ricker
#rx: 0.140 0.170 0

#box: 0 0 0 0.240 0.170 0.002 half_space
#cylinder: 0.120 0.080 0 0.120 0.080 0.002 0.010 pec

#geometry_view: 0 0 0 0.240 0.190 0.002 0.002 0.002 0.002 cylinder_half_space n

W pliku opisujemy cały model, który chcemy obliczyć.

#title: wpisujemy nazwę i przedmiot modelowania

#domain: definiujemy całą przestrzeń w jakiej będziemy modelować, wartości x,y,z

#dx_dy_dz: wielkość najmniejszej niepodzielnej komórki, zależne od częstotliwości anteny

#time _window: okno czasowe w jakiej będziemy modelować w tym przypadku 3 ns

#time_step_limit_type: rodzaj przestrzeni

#pml_cells: ilość komórek na brzegach modelu odpowiedzialnych za absorpcję fal

#material: parametry materiałowe przenikalność elektryczna, magnetyczna, oporność

#waveform: tutaj definiujemy rodzaj sygnału w tym przypadku Ricker'a i częstotliwość centralną naszej teoretycznej anteny georadarowej 1.5 GHz

#hertzian_dipole: pozycja anteny nadawczej sygnału x,y,z i rodzaj sygnału (wcześniej zdefiniowany)

#rx: pozycja anteny odbiorczej

#box: rozmiary naszej półprzestrzeni (ośrodek gruntowy/powietrze) i wsp. x,y,z początek układy i max zasięg modelu. Musi zawierać w naszej przestrzeni wcześniej zdefiniowanej #domain

#cylinder: parametry naszej rury (cylindra) x,y,z i na końcu rodzaj materiału (pec)

#geometry_view: funkcja ta generuje podgląd geometrii opisanej w powyższym tekście, na końcu podajemy nazwę pliku pod jakim ma być zapisany

Pliki wynikowe

Wygenerowany został sygnał Ricker'a o częstotliwości środkowej 1.5 GHz:

ricker
Źródło: www.gprmax.com

Podgląd geometrii modelu Tx - antena nadawcza. Rx antena odbiorcza;

geometria_ascan1
Model, opracowanie własne

Po obliczeniu teoretycznej odpowiedzi impulsowej ośrodka otrzymujemy wymodelowane 3 składowe elektryczne i magnetyczne:

Ascan_out
Opracowanie własne w gprMax

Podczas wykonywania obliczeń można dodać funkcję python #snapshot, dzięki której możemy co pewien ustalony ułamek czasu zrobić zrzut ekranu obecnej fazy modelowania. Dzięki temu  po połączeniu w klatki otrzymujemy wizualizację rozchodzenia się fali w rozpatrywanym ośrodku. W animacji pokazana jest tylko składowa Ez.

Animacja wygenerowana za pomocą gprMax i ParaView.