Od nadciśnienia do ugięcia płyty: numeryczne modelowanie skutków fali podmuchowej

Streszczenie

Nadciśnienie fali podmuchowej jest dopiero początkiem analizy. Dla konstrukcji równie ważne są impuls, czas dodatniej fazy, kąt padania, sposób podparcia elementu i to, czy materiał pracuje sprężyście, plastycznie, czy zaczyna lokalnie pękać. Dlatego przejście od wykresu podmuchu do odpowiedzi ściany albo płyty wymaga osobnego modelu mechaniki konstrukcji.¹^,²

Artykuł Krzysztofa Kosiuczenki o symulacji LS-DYNA pokazuje tę drogę na prostym, ale bardzo pouczającym przykładzie: stalowa płyta 1 x 1 m, grubości 10 mm, zamocowana w sztywnej ramie, jest obciążana falą od małego wybuchu jądrowego modelowanego funkcją Brode'a. Dla przypadku W = 0,25 kt i wysokości detonacji HOB = 200 m autor podaje nadciśnienie przy ground zero Ps = 52,03 kPa, a pełna symulacja numeryczna daje maksymalne ugięcie płyty około 1,04 mm.¹

To jest dobry przykład dydaktyczny, bo pokazuje różnicę między kilkoma poziomami modelu: empiryczną falą podmuchową, uproszczoną odpowiedzią jednego stopnia swobody, modalną odpowiedzią płyty oraz pełną analizą FEM z modelem materiałowym Johnson-Cook. Kalkulator odporności płyty celowo znajduje się między SDOF a FEM: przyjmuje zadany przebieg p(t), sumuje mody własne sprężystej płyty prostokątnej i jawnie pokazuje, czego nadal nie liczy bez LS-DYNA albo podobnego kodu.

Rozszerzenie tematu

W popularnym opisie skutków wybuchu pojawia się zwykle jedna liczba: nadciśnienie. To wygodne, bo można powiedzieć, że pewien próg uszkadza szyby, inny zagraża lekkim konstrukcjom, a jeszcze inny powoduje ciężkie zniszczenia. Takie progi są użyteczne w kalkulatorze stref zniszczeń i w kalkulatorze podmuchu, ale inżynieria konstrukcyjna nie kończy się na odczytaniu piku ciśnienia.²^,³

Płyta albo ściana odczuwa obciążenie w czasie. Krótki pik o dużej amplitudzie może dać mniejszą deformację globalną niż dłuższy impuls o niższym maksimum, jeśli całkowity pęd przekazany konstrukcji jest inny. Dlatego w analizach pojawia się przestrzeń pressure-impulse: oś nadciśnienia i oś impulsu, a nie jedna skala siły wybuchu. Ten sam tok rozumowania stoi za kalkulatorem pressure-impulse, który bierze falę podmuchową i zamienia ją na uproszczone ugięcie przegrody.

Praca Kosiuczenki dodaje do tego kolejny krok: zamiast zatrzymać się na jednym stopniu swobody, używa LS-DYNA, funkcji *LOAD_BRODE, siatki elementów powłokowych oraz modelu materiałowego Johnson-Cook. Brode opisuje falę podmuchową empirycznie, a Johnson-Cook opisuje stal przy dużych szybkościach odkształcenia i możliwym ogrzewaniu materiału. To już nie jest wykres progu skutków, lecz pełniejsza symulacja mechaniczna elementu.¹^,⁴

Najważniejszy etap pośredni wygląda tak:

określić energię wybuchu, wysokość detonacji i odległość celu,
przeliczyć geometrię przez skalowanie Hopkinsona-Cranza,
wyznaczyć nadciśnienie i impuls fali przy elemencie,
rozstrzygnąć, czy element widzi falę padającą, odbitą, czy wariant pośredni zależny od kąta,
zamienić przebieg ciśnienia na obciążenie płyty,
policzyć odpowiedź dynamiczną elementu,
dopiero potem interpretować ugięcie, naprężenie, uplastycznienie i uszkodzenie.

Wysokość detonacji komplikuje etap trzeci. Dla wybuchu powietrznego fala nie dociera do gruntu tak samo jak przy wybuchu kontaktowym. Pojawia się odległość skośna, odbicie od powierzchni, a w pewnych zakresach także przejście do frontu Mach stem. Dlatego kalkulator height-of-burst jest naturalnym poprzednikiem analizy płyty: odpowiada na pytanie, jaka fala w ogóle przychodzi do miejsca, w którym znajduje się konstrukcja.

W symulacji numerycznej Kosiuczenki obiektem jest prostokątna płyta stalowa zamocowana na obwodzie. To nie jest pełny budynek, ale taki model ma sens jako elementarz. Płyta pozwala oddzielić trzy zjawiska, które w prawdziwym obiekcie mieszają się ze sobą: rozkład ciśnienia na powierzchni, falę naprężeń w materiale oraz warunki brzegowe na krawędziach. Gdy płyta jest sztywno zamocowana, część pracy przejmują reakcje ramy; gdy jest podparta słabiej, rośnie ugięcie i znaczenie połączeń.

Model Johnson-Cook, użyty w pracy jako model materiału, ma postać mnożnikową: granica płynięcia zależy od odkształcenia plastycznego, szybkości odkształcenia i temperatury. Tego nie potrafi ani prosty SDOF, ani obecny kalkulator modalny. Kalkulator może policzyć, że płyta o danej grubości ma określoną masę powierzchniową, sztywność zginania, mody własne, ugięcie środka i sprężyste naprężenia z krzywizny, ale nie zobaczy lokalnego uplastycznienia przy krawędzi, wzrostu temperatury, uszkodzenia elementów siatki ani koncentracji naprężeń przy mocowaniu.¹^,⁵

W kalkulatorze odporności płyty uproszczenia są teraz rozdzielone inaczej. Kalkulator nie dobiera mocy źródła, wysokości detonacji ani odległości, tylko przyjmuje gotowy, bezpiecznie zdefiniowany przebieg p(t): Friedlandera, trójkątny, półsinusoidalny albo prostokątny. Konstrukcja jest liczona jako liniowo sprężysta płyta Kirchhoffa-Love'a podparta przegubowo na obwodzie, z superpozycją nieparzystych modów własnych i całkowaniem w czasie metodą Newmarka.

Taki wynik nie ma zastępować LS-DYNA. Jego wartość jest inna: pokazuje, które parametry są naprawdę ważne. Zwiększenie grubości działa bardzo silnie, bo sztywność zginania rośnie z trzecią potęgą grubości. Zmiana czasu trwania i kształtu p(t) potrafi zmienić pobudzenie modów nawet przy tym samym piku ciśnienia. Zwiększenie powierzchni płyty zwiększa masę całkowitą, ale o odpowiedzi lokalnej decyduje przede wszystkim rozpiętość, masa powierzchniowa i relacja czasu impulsu do okresów własnych. To są intuicje, których nie daje sama tabela progów nadciśnienia.

Trzeba też uważać na interpretację liczby 1,04 mm z pracy. Jest ona wynikiem konkretnej siatki, konkretnego zamocowania w sztywnej ramie, konkretnego modelu materiału i konkretnego przebiegu obciążenia. Jeżeli w kalkulatorze modalnym wynik wychodzi większy albo mniejszy, nie oznacza to automatycznie błędu. Oznacza to przede wszystkim, że model płyty przegubowo podpartej i liniowo sprężystej odpowiada na inne pytanie niż pełny FEM z warunkami brzegowymi, uplastycznieniem i lokalnymi koncentracjami naprężeń. W dydaktyce właśnie ta różnica jest najcenniejsza.

Ten przykład dobrze łączy też dwa światy serwisu. Z jednej strony są kalkulatory skutków wybuchu: równoważnik TNT, podmuch, height-of-burst, pressure-impulse i odporność płyty. Z drugiej strony są artykuły o fali odbitej i obciążeniu konstrukcji, odporności konstrukcji na wybuch oraz równoważniku trotylowym, nadciśnieniu i impulsie. Razem tworzą ścieżkę od fizyki fali do zachowania konkretnego elementu.

Najkrótszy wniosek jest taki: kalkulator podmuchu mówi, jaka fala dochodzi do celu; kalkulator pressure-impulse mówi, jak fala obciąża przegrodę w modelu SDOF; kalkulator odporności płyty pokazuje, jak zadane p(t), geometria, materiał i mody własne przekładają się na ugięcie oraz naprężenia. Pełny FEM jest następnym poziomem, a nie zamiennikiem dla zrozumienia tych prostszych kroków.

Teoria skalowania Hopkinsona-Cranza — podstawa porównań

Zanim dotrze się do odpowiedzi płyty, trzeba wiedzieć, jakie obciążenie na nią działa. Tutaj kluczowe jest prawo skalowania Hopkinsona-Cranza (1915/1926), które pozwala przenosić wyniki testów jednej skali wybuchu na inne skale.

Zasada podobieństwa wybuchowego:

Dwa wybuchy o różnych energiach W₁ i W₂ dają te same parametry fali podmuchowej (nadciśnienie, impuls) w odległościach skalowanych proporcjonalnie do sześciennego pierwiastka energii:

Z = R / W^(1/3)    [m/kg^(1/3)]

gdzie Z to skalowana odległość (scaled distance), R to rzeczywista odległość [m] i W to masa ekwiwalentu TNT [kg].

To oznacza, że jeśli test z 1 kg TNT daje nadciśnienie P₀ w odległości R = 1 m (Z = 1 m/kg^(1/3)), to test z 1 000 kg TNT da to samo nadciśnienie P₀ w odległości:

R₂ = Z × W₂^(1/3) = 1 × 1000^(1/3) = 10 m

Prawo skalowania jest niezwykle potężne: pozwala budować bazy danych z testów laboratoryjnych i małoskalowych (bezpiecznych i tanich) oraz przenosić je na duże ładunki. Kosiuczenko korzysta z niego, obliczając parametry fali dla W = 0,25 kt = 250 000 kg TNT.

Ograniczenia prawa skalowania:

Skalowanie działa idealnie w warunkach odbicia od płaskiej powierzchni lub w wolnej atmosferze. Komplikuje się gdy:

Geometria środowiska jest złożona (budynki, doliny, warstwy atmosferyczne)
Wybuchy powierzchniowe vs. wybuchy powietrzne mają różne charaktery skalowania (dla wybuchu gruntowego W_eff = 2W)
Wybuchy podziemne tworzą inne zjawiska fizyczne (lejek, fala sejsmiczna)
Temperatura i ciśnienie atmosfery odbiegają od standardowych

Dla wybuchu jądrowego dochodzą dodatkowe komplikacje: promieniowanie cieplne modyfikuje atmosferę przed przybyciem fali uderzeniowej (efekt Mach disk), a w pobliżu epicentrum pojawia się faza ujemna kondensacji — zjawiska nieobecne przy konwencjonalnych ładunkach.¹^,⁴

Parametry fali podmuchowej — Brode, CONWEP i pomiary empiryczne

Parametry fali podmuchowej można wyznaczać różnymi metodami, z różną dokładnością i dla różnych zastosowań.

Funkcja Brode'a:

Harold Brode (1959, 1987) opracował analityczne aproksymacje parametrów fali podmuchowej od wybuchu jądrowego na podstawie obliczeń hydrodynamicznych. Funkcja Brode'a opisuje nadciśnienie szczytowe jako funkcję skalowanej odległości:

Dla wybuchu swobodnego w powietrzu (free-air):

Ps/P₀ = 6,7/Z³ + 1/Z + 0,3/Z^(2/3) - 0,04/Z^(1/3) + 0,0027   [dla Z > 0,3 m/kg^(1/3)]

gdzie P₀ to ciśnienie atmosferyczne (101,3 kPa), Z to skalowana odległość.

Brode jest szczególnie przydatny dla wybuchu jądrowego, bo obejmuje szerokie zakresy Z (od bliskiego do bardzo dalekiego pola).

CONWEP (Conventional Weapons Effects):

CONWEP to zbiór wzorów empirycznych US Army oparty na bazach danych testowych z konwencjonalnymi materiałami wybuchowymi. Stosowany głównie dla Z > 0,2 m/kg^(1/3) i dla TNT lub materiałów o znanych ekwiwalentach. Zaletą jest bardzo dobra dopasowanie do danych doświadczalnych w typowych zakresach bojowych. Stosowany w kalkulatorach Kingery-Bulmash.

Różnica między Brode a CONWEP:

Przy małych Z (bliskie pole) wyniki mogą się znacznie różnić — Brode jest kalibrowany na wybuchy jądrowe, CONWEP na konwencjonalne. Dla celów dydaktycznych ważne jest rozumieć, że obydwie metody są aproksymacjami i mają zakresy stosowalności.¹^,⁴^,⁵

Model pressure-impulse i diagramy izouszkodzenia

Diagram pressure-impulse (P-I lub P-I diagram) jest graficznym narzędziem do oceny odpowiedzi elementu konstrukcji na dowolną falę podmuchową.

Zasada:

Na osiach diagramu są: nadciśnienie szczytowe P [kPa] i impuls I [kPa·ms]. Dla każdego poziomu odpowiedzi (ugięcie/granica plastyczności/zawalenie) wyznacza się krzywą izouszkodzenia — parę (P, I), przy której odpowiedź elementu osiąga dokładnie ten poziom.

Krzywa ma charakterystyczne asymptoty:

Asymptota impulsowa (quasi-statyczna): przy bardzo długich falach (I duże, P małe) odpowiedź zależy tylko od całkowitego impulsu
Asymptota ciśnieniowa (impulsowa): przy bardzo krótkich falach (P duże, I małe) odpowiedź zależy tylko od nadciśnienia szczytowego

W obszarze pomiędzy asymptotami (dynamiczna odpowiedź) oba parametry mają znaczenie.

Praktyczne zastosowanie:

Punkt odpowiadający danej fali podmuchowej (P, I) jest porównywany z krzywą izouszkodzenia elementu:

Punkt poniżej krzywej → element przeżywa
Punkt powyżej krzywej → element jest zniszczony

Różne elementy (szyba, ściana ceglana, słup żelbetowy, stalowa płyta) mają różne krzywe izouszkodzenia. To pozwala na szybką ocenę, które elementy zawiodą jako pierwsze.²^,⁵

Model SDOF — od fizyki do równania ruchu

Model jednego stopnia swobody (SDOF — Single Degree of Freedom) jest mostem między empirycznym ładunkiem falowym a obliczeniową odpowiedzią elementu. To podstawowy narzędzie w podejściu inżynierskim przed FEM.

Idea modelu:

Element (płyta, ściana, belka) jest zastępowany układem masa-sprężyna-tłumik z jednym stopniem swobody — przemieszczeniem centrum lub punktu o największej amplitudzie:

M_eff × ẍ + C × ẋ + K_eff × x = F(t)

gdzie M_eff to masa efektywna (ułamek całkowitej masy elementu), K_eff to efektywna sztywność (ułamek sztywności statycznej), C to współczynnik tłumienia, x to przemieszczenie centrum, F(t) to siła od fali (ładunek × ułamek obciążenia).

Współczynniki transformacji:

Kluczowymi parametrami są współczynniki KLM = KL / KM, gdzie KL to współczynnik ładunku (Load Factor), a KM to współczynnik masy (Mass Factor). Zależą od schematu statycznego i rozkładu ugięcia:

Dla płyty prostokątnej utwierdzonej na obwodzie i obciążenia równomiernego:

KL ≈ 0,40 (dla pierwszej postaci drgań)
KM ≈ 0,25 (dla pierwszej postaci drgań)
KLM = 0,40/0,25 = 1,60

Rozwiązanie równania ruchu:

Dla fazy elastycznej przy uproszczonym ładunku trójkątnym (narastanie do P_s, następnie zanik do zera przez czas t_d):

x_max = DLF × x_static = DLF × F_max / K

gdzie DLF (Dynamic Load Factor) zależy od stosunku t_d/T (czas impulsu do okresu własnego elementu). Dla t_d/T << 1 (reżim impulsowy): DLF = π × i / (K × T). Dla t_d/T >> 1 (reżim quasi-statyczny): DLF ≈ 2.¹^,⁵

Model materiałowy Johnson-Cook — plastyczność przy wysokich prędkościach odkształcenia

W klasycznym modelu sprężysto-plastycznym granica plastyczności σ_Y jest stała. W rzeczywistości materiały metalowe wykazują wrażliwość na szybkość odkształcenia i temperaturę — szczególnie przy dynamicznych obciążeniach wybuchowych.

Wzór Johnson-Cook:

σ = (A + B × ε^n) × (1 + C × ln(ε̇/ε̇₀)) × (1 - T_H^m)

gdzie:

A — granica plastyczności [MPa]
B — współczynnik umocnienia [MPa]
n — wykładnik umocnienia
C — współczynnik wrażliwości na prędkość odkształcenia
ε — odkształcenie plastyczne
ε̇ — prędkość odkształcenia [1/s]
ε̇₀ — referencyjna prędkość odkształcenia (zwykle 1,0 /s)
T_H — zredukowana temperatura: (T - T_ref) / (T_melt - T_ref)
m — wykładnik termiczny

Parametry dla stali S235:

Typowe wartości dla stali niskowęglowej (zbliżone do stali konstrukcyjnej S235/A36):

A ≈ 290 MPa, B ≈ 400 MPa, n ≈ 0,28
C ≈ 0,013, m ≈ 1,03
T_ref = 294 K (temperatura pokojowa), T_melt ≈ 1793 K

Te parametry oznaczają, że przy prędkości odkształcenia 1 000 s⁻¹ (typowej dla wybuchu) granica plastyczności wzrasta o czynnik (1 + 0,013 × ln(1000)) ≈ 1,09 — czyli o ok. 9% powyżej wartości quasi-statycznej.¹^,⁵

Metoda elementów skończonych dla dynamiki udarowej

LS-DYNA (i podobne kody: ABAQUS/Explicit, MSC Dytran) rozwiązuje równania ruchu w dziedzinie czasu metodą jawną (explicit time integration). Różni się tym od statycznych obliczeń MES, które rozwiązują równowagę statyczną.

Schemat całkowania jawnego:

Prędkość i przemieszczenie są obliczane przez całkowanie równania siły:

a(t) = F_net(t) / M
v(t+Δt/2) = v(t-Δt/2) + a(t) × Δt
u(t+Δt) = u(t) + v(t+Δt/2) × Δt

Jest to kondycjonalnie stabilne — wymaga kroku czasowego Δt < T_kryt = L_min / c, gdzie L_min to minimalna długość elementu siatki i c to prędkość fali w materiale. Dla stali c ≈ 5 000 m/s — przy elementach 1 cm wymaga kroków ok. 2 µs.

Elementy powłokowe:

Płyta w LS-DYNA jest zwykle modelowana elementami powłokowymi (shell elements). Elementy powłokowe zakładają płaski stan naprężenia w płaszczyźnie i trójwymiarowy schemat ugięcia. Typowo używa się 4-węzłowych elementów Q4 lub 3-węzłowych T3, z całkowaniem numerycznym przez grubość (np. 5 punktów całkowania przez grubość dla modelu naprężeń).

Siatka i zbieżność:

Wynik symulacji zależy od gęstości siatki. Gęsta siatka daje dokładniejszy wynik, ale jest kosztowna obliczeniowo. Praca Kosiuczenki nie podaje szczegółów siatki dla płyty 1×1 m, ale typowo dla analizy udarowej stosuje się elementy 5–20 mm, co dla tej płyty daje siatkę od 2 500 do 40 000 elementów.¹

Polska perspektywa — WAT i badania odporności budowli

W Polsce badania odporności konstrukcji na wybuchy prowadzone są przede wszystkim w Wojskowej Akademii Technicznej (WAT) w Warszawie, której pracownicy publikują w dziedzinie mechaniki udarowej i inżynierii ochrony.

Kierunek badań WAT:

WAT realizuje badania zarówno akademickie (jak praca Kosiuczenki), jak i aplikacyjne dla Ministerstwa Obrony Narodowej. Obszary obejmują:

Modelowanie numeryczne obciążeń wybuchowych (LS-DYNA, AUTODYN, ABAQUS)
Badania doświadczalne na strzelnicy WAT (ładunki do kilkudziesięciu kg TNT)
Projektowanie ochron pasywnych dla pojazdów wojskowych
Ocena odporności budowli wojskowych i schronów

Normy i standardy projektowania:

Projektowanie budynków odpornych na wybuchy w Polsce (dla celów wojskowych) opiera się na normach NATO i dokumentach Allied Publications (AP). Dla celów cywilnych (ochrona przed zamachami) stosuje się Eurocode 1 (EN 1991-1-7: Oddziaływania wyjątkowe) i British Standard BS 6349/PAS 68.

Związek z infrastrukturą jądrową:

Planowane obiekty jądrowe w Polsce (EJ Pomerania) będą projektowane z odpornością na wybuchy zewnętrzne (sabotaż, ataki terrorystyczne, wypadki lotnicze) jako część wymogów bezpieczeństwa IAEA — ściśle powiązanych z tematem obrony w głąb.¹

Trzy przykłady numeryczne

Przykład 1: Skalowanie odległości dla W = 0,25 kt

Kosiuczenko podaje W = 0,25 kt = 250 t TNT = 250 000 kg TNT. Dla przypadku ground zero (odległość pozioma = 0) przy HOB = 200 m, odległość skośna to:

R = HOB = 200 m
Z = R / W^(1/3) = 200 / (250 000)^(1/3) = 200 / 63,0 = 3,17 m/kg^(1/3)

Przy Z = 3,17 m/kg^(1/3) formuła CONWEP dla wybuchu sferycznego daje nadciśnienie:

ln(Ps) ≈ -0,0427 × (ln Z)^3 + 1,019 × (ln Z)^2 - 4,615 × ln Z + 7,266
ln(3,17) ≈ 1,154
ln(Ps) ≈ -0,0427 × 1,538 + 1,019 × 1,331 - 4,615 × 1,154 + 7,266 ≈ 52,6 kPa

Wynik ok. 52 kPa zgadza się z wartością Ps = 52,03 kPa podaną przez Kosiuczenkę — potwierdzenie poprawności obliczeń.¹^,⁴

Przykład 2: Okres własny płyty i reżim odpowiedzi

Dla płyty stalowej 1×1 m, grubość h = 10 mm, materiał: E = 210 GPa, ν = 0,3, ρ = 7 850 kg/m³.

Okres własny płyty utwierdzonej na obwodzie:

T = 2π / ω₁, gdzie ω₁ = λ₁² × √(D / (ρ·h))

λ₁ ≈ π²(1 + 1) = 2π² / a² (dla płyty kwadratowej, a = 1 m, przybliżenie)
D = E × h³ / (12(1-ν²)) = 210×10⁹ × (0,01)³ / (12 × (1-0,09)) = 210×10³ × 10⁻⁶ / 10,92 ≈ 19,23 N·m

Masa jednostkowa: m = ρ × h = 7 850 × 0,01 = 78,5 kg/m²

ω₁² ≈ (2π²)² × D / (m × a⁴) = 39,48 × 19,23 / (78,5 × 1) ≈ 9,67 rad²/s²
ω₁ ≈ 3,11 rad/s → T ≈ 2,02 s

Czas impulsu wybuchu: t_d ≈ 40 ms (typowe dla Z ≈ 3 m/kg^(1/3)).

Stosunek: t_d/T = 0,040/2,02 ≈ 0,02 — reżim impulsowy (t_d << T). W reżimie impulsowym ugięcie zależy od impulsu, nie nadciśnienia — co jest kluczowym wynikiem analizy P-I.¹^,⁵

Przykład 3: Grubość krytyczna — granica między ugięciem a zniszczeniem

Dla płyty S235 (σ_Y = 235 MPa) 1×1 m, ładunek impulsowy I = 300 kPa·ms:

Masa efektywna (SDOF): M_eff = KM × ρ × h × a² = 0,25 × 78,5 × h × 1 = 19,6h [kg]

Prędkość początkowa (przy impulsie I): v₀ = (KL × I × a²) / M_eff = (0,40 × 0,3 × 1) / (0,25 × 78,5 × h) = 0,120 / (19,6h)

Maksymalna energia kinetyczna przekształca się w energię odkształcenia. Dla płynięcia:
Warunek zniszczenia: ugięcie > 10 × h (typowy kryterium dużych odkształceń)

Rozwiązanie iteracyjne dla tej klasy problemu daje h_kryt ≈ 4–6 mm — przy grubości 10 mm płyta ma znaczny zapas przed zniszczeniem przy typowych parametrach fali Z ≈ 3 m/kg^(1/3). To jest właśnie klasa wyników, którą Kosiuczenko potwierdza numerycznie na poziomie ok. 1 mm ugięcia.¹

Podsumowanie dydaktyczne

Temat modelowania numerycznego odpowiedzi płyty na wybuch łączy kilka dyscyplin, które rzadko są razem omawiane w jednym kursie.

Hierarchia modeli. Kluczowym przesłaniem jest hierarchia modeli: empiryczny opis fali (Brode/CONWEP) → diagram P-I → SDOF → model modalny płyty → FEM. Każdy poziom jest użyteczny do innego pytania i ma inne granice stosowalności. Rozumieć, kiedy wystarczy prostszy model, a kiedy potrzebny jest FEM — to jest kompetencja inżynierska, której nie daje sama umiejętność obsługi LS-DYNA. Dla doktoranta zajmującego się tą dziedziną po raz pierwszy, zalecana jest droga od góry: zacznij od tabeli progów, przejdź przez diagram P-I i SDOF, porównaj z modelem modalnym, dopiero potem wejdź w FEM. Odwrotna droga (FEM bez rozumienia prostszych poziomów) zwykle kończy się ufaniem wynikom, których się nie rozumie.

Fizyka a narzędzie. Model Johnson-Cook jest narzędziem. Ale zrozumieć, dlaczego stal ma inne właściwości przy 1 000/s niż przy 0,001/s, wymaga fizyki ciała stałego i nauki o materiałach. Kurs numeryczny bez fizyki daje użytkownika oprogramowania, a nie inżyniera rozumiejącego co oblicza.

Ekstrapolacja i błąd. Wynik „1,04 mm" jest wynikiem jednego konkretnego modelu. Przed użyciem go jako projektowego parametru trzeba rozumieć niepewności: siatki, modelu materiałowego, warunków brzegowych, parametrów fali. Wiarygodny wynik FEM zawsze towarzyszą badania zbieżności i walidacja.

Polska nauka w światowym kontekście. Praca Kosiuczenki (WAT, 2025) jest dobrym przykładem polskiej produkcji akademickiej w tej niszy — łączącej fizykę wybuchu jądrowego z obliczeniami konstrukcyjnymi. Publikacja w polskim czasopiśmie nie oznacza izolacji — metodologia jest identyczna z pracami publikowanymi w International Journal of Impact Engineering czy Shock Waves. Polscy doktoranci w tej dziedzinie działają na tym samym poziomie metodycznym co zachodnie ośrodki badawcze.¹^,⁴^,⁵

Warunki brzegowe i ich znaczenie dla wyników

Jednym z najważniejszych, a często niedocenianych aspektów modelowania odpowiedzi elementu na wybuch jest dobór warunków brzegowych. Sama geometria płyty i materiał to dopiero połowa historii.

Sztywne zamocowanie vs. podparcie swobodne:

Praca Kosiuczenki zakłada sztywne zamocowanie na obwodzie (fixed/clamped boundary). To przypadek konserwatywny dla wartości naprężeń na krawędziach — ale może być nie-konserwatywny dla ugięcia w centrum. Jeśli rzeczywiste połączenie płyty z ramą jest podatne (bolted joint, spawana krawędź z podatną ramą), ugięcie centrum może być większe.

Warunek brzegowy	Max ugięcie	Max naprężenie	Gdzię naprężenie szczytowe
Utwierdzony (clamped)	Najmniejsze	Największe	Na krawędziach
Podparte swobodnie (simply supported)	Największe	Mniejsze	W centrum
Mieszane (elastyczna granica)	Pośrednie	Pośrednie	Rozkład hybrydowy

Konsekwencje projektowe:

Jeśli projektujemy płytę pod kątem wytrzymałości materiału (niezniszczenia koszulek spawów), powinniśmy użyć warunków clamped — bo generują największe momenty zginające przy zamocowaniu. Jeśli zależy nam na ograniczeniu ugięcia (szczelność) — powinniśmy użyć simply supported, bo generuje maksymalne ugięcie centrum. W praktyce prawdziwe zamocowanie leży pośrodku i wymaga odrębnej analizy węzła połączenia.¹^,⁵

Fala odbita i stojąca — komplikacje przy geometriach zamkniętych

W otwartej przestrzeni fala podmuchowa przychodzi raz i mija element. W zamkniętej lub częściowo zamkniętej przestrzeni (pomieszczenie, schronu, wnętrze pojazdu) fala odbija się od ścian i tworzy złożone pole ciśnienia.

Faza wentylacji vs. faza odbić:

Po detonacji wewnątrz lub wejściu fali do wnętrza przez otwór, ciśnienie wewnątrz narasta szybko (faza uderzeniowa), potem powoli spada przez wentylację (wypływ przez otwory) — ale po drodze wielokrotnie odbiją się od ścian. Każde odbicie tworzy nowe ekstrema ciśnienia.

Czas charakterystyczny oczyszczenia (clearance time) zależy od stosunku wolnej przestrzeni do pola otworów: t_c ≈ 4V/(A_open × c), gdzie V to objętość, A_open to pole otworów, c to prędkość dźwięku. Dla szczelnego schronu t_c → ∞ — ciśnienie wewnątrz rośnie do maksimum i spada tylko przez nieszczelności.

Quasi-statyczne ciśnienie wewnętrzne:

Dla zamkniętej przestrzeni definiuje się ciśnienie quasi-statyczne (QSP), które można w uproszczeniu obliczyć z zasady zachowania energii:

P_QS ≈ P₀ × (E_wybuch / (V × P₀/γ) + 1)

gdzie E_wybuch to energia fali wewnątrz, V to objętość pomieszczenia, γ ≈ 1,4 dla powietrza. QSP jest ważne przy projektowaniu schronów atomowych.²^,⁵

Oddziaływanie fali z materiałami wielowarstwowymi i pochłaniaczami

W praktyce schrony i pojazdy wojskowe stosują wielowarstwowe panele: stal + pianka energochłonna + kompozyt + ceramika. Odpowiedź takiego układu na podmuch jest znacznie bardziej złożona niż dla jednorodnej płyty.

Mechanizmy pochłaniania energii:

Plastyczne zgniatanie pianki (foam crushing): Piana aluminiowa lub metalowa o otwartych porach zgniatana jest w szerokim zakresie odkształceń plastycznych przy prawie stałym ciśnieniu (plateau stress). Energia jest pochłaniana przez deformację komórek.
Kruchość ceramiki: Ceramika (AI₂O₃, SiC, B₄C) jest bardzo twarda i pochłania energię przez kruche pękanie — skuteczne przy uderzeniach punktowych, ale gorzej radzi sobie z rozłożonym podmuchem.
Efekt sandwiczu: Panele sandwich (dwie okładziny stalowe + rdzeń piankowy) mają znacznie wyższy opór na podmuch niż sama suma masy okładzin i rdzenia — ze względu na efekt dźwigni zginania.

Obciążenie uderzeniowe vs. podmuchowe:

Te same warstwy chronią inaczej przed odłamkiem (problem penetracji, lokalne naprężenia) a inaczej przed podmuchem (problem globalnego ugięcia i momentów zginających). Optymalny projekt dla obu zagrożeń jednocześnie jest zazwyczaj kompromisem.⁵

Porównanie metod obliczeniowych — kiedy co stosować

Inżynier staje przed wyborem metody obliczeniowej na każdym etapie projektu. Tabela poniżej porządkuje wybory.

Metoda	Szybkość	Dokładność	Kiedy stosować
Progi nadciśnienia (tabele)	Sekundy	Niska	Wstępna ocena, skala skutków
Diagram P-I	Minuty	Umiarkowana	Szybka ocena graniczna
SDOF (kalkulator)	Minuty	Umiarkowana	Projekt wstępny, analiza wrażliwości
SDOF (zaawansowany, z niespre.)	Godziny	Dobra	Projekt inżynierski dla prostych elementów
FEM explicit (LS-DYNA)	Godziny–doby	Wysoka	Projekt finalny, walidacja, badania
FEM + CFD (sprzężone)	Dni–tygodnie	Bardzo wysoka	Badania naukowe, kalibracja modeli

Kiedy FEM jest konieczny:

Geometria elementu jest złożona (otwory, wzmocnienia, kształtki)
Interesuje nas lokalne zachowanie (naprężenia przy spawie, pęknięcia)
Obciążenie nie jest jednolite przestrzennie
Materiał ma złożone zachowanie (duże odkształcenia, pęknięcia, kontakt)
Wynik będzie podstawą certyfikacji lub licencji regulacyjnej

Kiedy SDOF wystarczy:

Prostokątna płyta/ściana bez otworów
Weryfikacja, czy element w ogóle jest zagrożony
Analiza wrażliwości parametrycznej (jak zmienia się wynik przy różnych grubościach)
Nauczanie i rozumienie fizyki odpowiedzi¹^,⁵

Efekty dynamiczne poza modelem płyty

Artykuł Kosiuczenki skupia się na płycie jako izolowanym elemencie. W realnym obiekcie pojawiają się efekty, które ten model pomija:

Strukturalne sprzężenie elementów:

Ściana przenosi obciążenie na stropy, słupy i fundamenty. Jeśli jedna ściana ugnie się plastycznie, obciąży inaczej słupy nośne. Progresywne zawalenie (progressive collapse) jest możliwe nawet gdy żaden element nie zawiódł indywidualnie — zawiodła sieć nośna jako całość.

Ruch gruntu:

Przy naziemnych i podziemnych wybuchach fala przechodzi przez grunt i wywołuje drgania fundamentów (Ground Shock). To obciążenie może być dominujące dla głęboko posadowionych struktur (schrony). W modelu Kosiuczenki ładunek powietrzny jest podawany bezpośrednio na powierzchnię płyty — ale dla schronu podziemnego dominowałoby obciążenie gruntowe.

Thermal pulse:

Dla wybuchu jądrowego przed falą uderzeniową nadchodzi puls termiczny. Przy dawkach > 100 kJ/m² (w pobliżu ogniska) materiały powierzchniowe mogą ulec topnieniu lub zapaleniu. Zmieniona przez temperaturę powierzchnia stali ma inne właściwości przy uderzeniu fali — efekt pomijany w standardowych analizach podmuchowych.¹^,⁴

Historia i tło intelektualne modelowania wybuchowego

Modelowanie numeryczne odpowiedzi konstrukcji na wybuchy ma historię sięgającą lat 40. i 50. XX wieku, kiedy rządy USA i ZSRR zaczęły systematycznie badać skutki broni jądrowej.

Program Nevada Test Site:

W USA w Nevada Test Site (NTS) przeprowadzono setki testów z udziałem pełnoskalowych budynków, pojazdów i schronów. Wyniki tych testów stały się podstawą empirycznych baz danych, z których wywodzą się dzisiejsze normy (TM 5-1300, UFC 3-340-02). Program Operation Doorstep (1953) testował typowe domy mieszkalne przy wybuchu 16 kt — dokumentacja fotograficzna i filmowa tych testów jest dostępna w archiwach i jest niezwykłym wizualnym dowodem, co robi podmuch z konwencjonalnymi konstrukcjami.

Kod LS-DYNA — historia:

LS-DYNA wywodzi się z kodu DYNA3D, opracowanego przez Johna Halsquista w Lawrence Livermore National Laboratory w latach 70. dla potrzeb modelowania uderzenia w pojazdy i analiz bezpieczeństwa broni. Livermore Software Technology Corporation (LSTC, obecnie Ansys) rozwijał LS-DYNA przez dekady, robiąc z niego dominujący kod dla analiz jawnych w dziedzinie czasu w dziedzinie automotive, obronności i inżynierii cywilnej.

Kod jest komercyjny, ale akademickie licencje są dostępne. WAT i wiele polskich uczelni technicznych (Politechnika Warszawska, AGH) posiada licencje i stosuje LS-DYNA w badaniach.

Od testów fizycznych do symulacji wirtualnych:

Historyczna sekwencja jest wyraźna:

Lata 40.–50. — testy fizyczne dużej skali (bomby, ładunki)
Lata 60.–70. — empiryczne formuly i tabele jako przetworzone dane testowe
Lata 80.–90. — kody numeryczne (DYNA, AUTODYN) zaczynają zastępować część testów
Lata 2000.–dziś — pełne wirtualne laboratoria wybuchowe, rzadkie testy fizyczne głównie jako walidacja

Polska nauka dołączyła do trzeciego i czwartego etapu — właśnie takim przykładem jest praca Kosiuczenki.¹^,⁴^,⁵

Konteksty zastosowań — cywilne i wojskowe

Techniki modelowania odpowiedzi konstrukcji na wybuchy mają liczne zastosowania poza czysto wojskowymi.

Inżynieria bezpieczeństwa cywilnego:

Po zamachach 11 września 2001 i wielu wcześniejszych (Oklahoma City 1995 — ładunek 2 ton ANFO w pojazdem) rządy na całym świecie zaczęły nakładać wymagania dotyczące odporności na wybuchy dla:

Budynków rządowych i dyplomatycznych
Infrastruktury krytycznej (elektrownie, mosty, lotniska)
Budynków wysokiego ryzyka (giełdy, centra danych)

W USA obowiązuje seria standardów GSA i DoS (Department of State) dla budynków rządowych i placówek dyplomatycznych. W UK — British Standard BS 8300 i PAS 97 dla ochrony przed VBIED (Vehicle-Borne Improvised Explosive Device).

Inżynieria górnicza i wyburzeniowa:

Kontrolowane wybuchy w górnictwie odkrywkowym, kamieniołomach i wyburzeniach budynków wymagają modeli fali podmuchowej i odpowiedzi gruntu. Typowe ładunki to kilkanaście do kilkuset kilogramów ANFO — małoskalowe w porównaniu z wybuchy militarnymi, ale wymagające precyzji ze względu na pobliskie zabudowania.

Wypadki i katastrofy przemysłowe:

Wybuchy gazów, pyłów (pył węglowy, mąka) czy parowych obłoków VCE (Vapor Cloud Explosion) — wszystkie wymagają modeli podobnych do opisanych w tym artykule. Katastrofa w Texas City (2005 — eksplozja rafineryjnych zbiorników BP, 15 ofiar) czy Port Bejrutu (2020 — 2 750 ton azotanu amonu, 218 ofiar) były analizowane dokładnie tymi metodami przez komisje badawcze.²^,⁵

Zaawansowane modele materiałowe — poza Johnson-Cook

Johnson-Cook jest standardem, ale nie jedynym modelem stosowanym w analizach wybuchowych.

Zerilli-Armstrong:

Alternatywny model oparty na fizyce dyslokacji, lepiej opisujący zachowanie metali przy ekstremalnie wysokich prędkościach odkształcenia (> 10⁵ /s). Stosowany w analizach penetracji i uderzenia.

Mohr-Coulomb i Drucker-Prager:

Modele dla materiałów gruntowych (glina, piasek, beton) — stosowane przy analizach wybuchu podziemnego i odpowiedzi gruntu. Grunt nie ma wyraźnej granicy plastyczności — plastyczność jest opisana przez kopertę Mohra-Coulomba (f = c + σ × tan φ).

Modele betonu:

Beton jest materiałem złożonym — kruchy pod rozciąganiem, relatywnie wytrzymały na ściskanie, z degradacją sprężystości po zarysowaniu. Modele takie jak CSCM (Continuous Surface Cap Model) lub K&C (Karagozian & Case) są stosowane w LS-DYNA dla betonowych elementów schronów.

Dla planowanej polskiej elektrowni jądrowej — gruba betonowa obudowa bezpieczeństwa i jej odpowiedź na wewnętrzne i zewnętrzne wybuchy (np. uderzenie samolotu symulowane jako lokalne obciążenie dynamiczne) jest analizowana właśnie tymi narzędziami przy certyfikacji projektu AP1000 przez NRC.¹^,⁵

Dostępność narzędzi obliczeniowych dla polskich badaczy

Warto oddzielnie omówić kwestię dostępności kodu LS-DYNA i podobnych narzędzi dla polskich uczelni i doktorantów.

Licencja komercyjna LS-DYNA jest kosztowna (kilkadziesiąt tysięcy dolarów rocznie), ale Ansys (właściciel od 2019) oferuje programy akademickie w obniżonej cenie. WAT posiada licencję instytucjonalną i stosuje kod zarówno w badaniach (jak praca Kosiuczenki), jak i w dydaktyce na kierunkach mechatroniki i inżynierii bezpieczeństwa. Politechnika Warszawska i AGH posiadają własne instalacje dla potrzeb badań naukowych.

Alternatywy open-source dla LS-DYNA w dziedzinie dynamiki udarowej są ograniczone, ale istnieją:

OpenFOAM — dla CFD i rozchodzenia się fal w płynach, nie dla dynamiki ciał stałych
FEniCS / OpenGeoSys — dla geomechaniki, nie dla szybkiej dynamiki
ABAQUS/Explicit — komercyjny, ale częsty w środowiskach akademickich

Dla doktoranta WAT z dostępem do kodu — zadanie zreprodukowania wyników Kosiuczenki jest technicznie wykonalne jako ćwiczenie laboratoryjne i stanowi doskonały punkt wyjścia do własnych badań parametrycznych. Kontekst naukowy pracy — modelowanie wybuchu jądrowego jako ładunku dydaktycznego — jest szczególnie wartościowy właśnie dlatego, że łączy fizykę jądrową z inżynierią mechaniczną w jednym, weryfikowalnym przykładzie.¹^,⁵

Dodatkowe materiały multimedialne

Warto dodać schemat pięciu poziomów modelu: fala Brode/CONWEP, przebieg pressure-time, model SDOF, model modalny płyty oraz siatka FEM z mapą ugięcia.

Powiązane kalkulatory i narzędzia

Odporność płyty na podmuch — kalkulator opisany w tym artykule.
Podmuch — parametry fali podmuchowej, nadciśnienie, impuls i fala odbita.
Pressure-Impulse dla ściany — przejście od fali do odpowiedzi jednego stopnia swobody.
Height-of-burst i Mach stem — wpływ wysokości detonacji i odbicia od gruntu.
Strefy zniszczeń — progi skutków w funkcji mocy i odległości.
Równoważnik TNT — energia, równoważnik TNT i odległość skalowana.
Odłamki — osobny problem penetracji fragmentów, który nie jest liczony w modelu płyty.

Ćwiczenia praktyczne

Pierwsze ćwiczenie powinno polegać na porównaniu przypadku Kosiuczenki z modelem modalnym, bez udawania pełnej zgodności geometrii podparcia. Należy:

otworzyć kalkulator odporności płyty,
ustawić p_max = 52,03 kPa, wybrać przebieg Friedlandera, płytę 1 x 1 m, grubość 10 mm i materiał stalowy,
zmieniać t_d, tłumienie i liczbę modów, obserwując ugięcie, naprężenia oraz test zbieżności,
porównać wynik z wartością 1,04 mm z LS-DYNA,
wskazać, które różnice wynikają z podparcia przegubowego, liniowej sprężystości, braku Johnson-Cook i braku rzeczywistej ramy.

Drugie ćwiczenie powinno dotyczyć wpływu grubości. Należy:

policzyć tę samą płytę dla grubości 6 mm, 10 mm, 20 mm i 40 mm,
porównać masę powierzchniową i sztywność zginania,
sprawdzić, czy ugięcie spada liniowo, czy szybciej,
powiązać obserwację z faktem, że sztywność płyty rośnie z trzecią potęgą grubości,
wyjaśnić, dlaczego sama masa płyty nie jest pełnym opisem odporności.

Trzecie ćwiczenie powinno porównać kształt obciążenia w czasie. Należy:

ustawić tę samą wartość p_max i t_d dla przebiegu Friedlandera, trójkąta, półsinusa i prostokąta,
zapisać impuls dodatniej fazy dla każdego przebiegu,
porównać maksymalne ugięcie i naprężenie płyty,
sprawdzić, kiedy wynik jest bliski odpowiedzi pierwszego modu, a kiedy rośnie udział modów wyższych,
osobno użyć kalkulatora height-of-burst tylko do dyskusji, skąd w realnej analizie wzięłoby się kandydackie p(t).

Przejdź do ćwiczenia interaktywnego

Powiązane artykuły

Przypisy

Krzysztof Kosiuczenko, Numerical simulation of nuclear blast wave impact on protective elements of anthropogenic structures, Inżynieria Bezpieczeństwa Obiektów Antropogenicznych 1 (2025), s. 12-22. Zakres wykorzystania: funkcja Brode'a, przypadek 0,25 kt i HOB = 200 m, płyta stalowa 1 x 1 m x 10 mm, LS-DYNA, Johnson-Cook, wynik ugięcia około 1,04 mm.
↩
Równoważnik trotylowy, nadciśnienie i impuls: jak naprawdę porównuje się siłę podmuchową materiałów wybuchowych. Zakres wykorzystania: rozdzielenie energii, nadciśnienia, czasu dodatniej fazy i impulsu.
↩
Fala odbita i obciążenie konstrukcji: nadciśnienie odbite, impuls i znaczenie geometrii V-shape przy wybuchu pod pojazdem. Zakres wykorzystania: różnica między falą padającą i odbitą oraz znaczenie geometrii celu.
↩
H. L. Brode, Airblast from nuclear bursts - analytic approximations, Defense Nuclear Agency, 1987. Zakres wykorzystania: empiryczna funkcja nadciśnienia od wybuchu jądrowego, przywołana przez Kosiuczenkę jako podstawa obciążenia w LS-DYNA.
↩
G. R. Johnson, W. H. Cook, A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures, 7th International Symposium on Ballistics, 1983. Zakres wykorzystania: sens modelu materiałowego Johnson-Cook i powód, dla którego prosty kalkulator SDOF nie zastępuje pełnej symulacji materiałowej.
↩