Metoda elementów skończonych w przeglądarce internetowej:
podstawy teoretyczne, budowa solvera TOMes® i weryfikacja
Praca opisuje budowę i weryfikację systemu PoSMES — samodzielnego programu do analizy wytrzymałościowej metodą elementów skończonych, działającego w przeglądarce internetowej bez instalacji, bez zależności zewnętrznych i bez przesyłania modelu na serwer. Częścią systemu jest silnik obliczeniowy TOMes®, napisany w całości od podstaw: od formatu macierzy rzadkiej i faktoryzacji LDLᵀ, przez elementy bryłowe TET4, TET10 i HEX8, po trzydzieści siedem rodzajów analiz obejmujących statykę, dynamikę, wyboczenie, plastyczność, kontakt, hipersprężystość, przewodnictwo ciepła, akustykę, zmęczenie, reologię oraz wypełnianie formy wtryskowej.
Omówiono sformułowanie metody, dobór elementów wraz z mechanizmami przeciwdziałającymi blokadzie objętościowej i ścinaniu pasożytniczemu, algorytmy rozwiązywania układów równań i zagadnienia własnego, oraz sformułowania nieliniowe. Osobny rozdział poświęcono weryfikacji: system zawiera zestaw 239 testów porównujących wynik z rozwiązaniami zamkniętymi, którego pełny przebieg jest przytoczony wraz z liczbami. Rozdział trzynasty przedstawia cztery przykłady rachunkowe policzone równolegle wzorem i metodą elementów skończonych.
Programy do analizy metodą elementów skończonych są dziś narzędziami dojrzałymi i drogimi. Licencja komercyjnego pakietu kosztuje tyle, co samochód, a jego instalacja wymaga stacji roboczej, uprawnień administratora i kilkudziesięciu gigabajtów miejsca. Dla studenta, konstruktora pracującego jednoosobowo i dla małego biura projektowego jest to bariera nieproporcjonalna do rzeczywistej potrzeby, którą zwykle da się opisać jednym zdaniem: czy ten wspornik wytrzyma i o ile się ugnie.
Celem pracy było zbudowanie systemu, który odpowiada na to pytanie rzetelnie, a jednocześnie jest dostępny natychmiast — bez instalacji, bez konta i bez wysyłania czyjejkolwiek dokumentacji na cudzy serwer. Zakres obejmuje kompletny łańcuch: od wczytania geometrii z systemu CAD, przez siatkowanie, sformułowanie i rozwiązanie układu równań, po prezentację wyniku i raport nadający się do dokumentacji.
Silnik obliczeniowy powstał w całości od podstaw. Nie jest nakładką na istniejący solver ani kompilacją cudzego kodu do postaci uruchamialnej w przeglądarce: format macierzy rzadkiej, faktoryzacja, algorytm przenumerowania, iteracja podprzestrzeni, całkowanie po czasie i modele materiałowe zostały napisane bezpośrednio w języku JavaScript i są przedmiotem rozdziałów od drugiego do dziesiątego.
Cztery decyzje podjęte na początku określiły kształt całego systemu.
Zadanie liczy się na komputerze użytkownika, a nie na serwerze. Konsekwencje są dalekosiężne i w większości korzystne: hosting nie musi mieć mocy obliczeniowej ani uruchamiać procesów w tle, więc systemem da się zarządzać na najtańszym koncie współdzielonym; koszt nie rośnie z liczbą użytkowników; a plik CAD — który w praktyce inżynierskiej bywa objęty umową o poufności — nigdzie nie jest wysyłany. Ceną jest wykorzystanie mocy jednej maszyny i jednego wątku, do czego wracamy w punkcie 11.2.
Silnik nie korzysta z żadnej biblioteki. Nie ma etapu budowania, transpilacji ani menedżera pakietów — pliki źródłowe są tym, co wykonuje przeglądarka. Ta decyzja kosztowała najwięcej pracy przy pisaniu solvera rzadkiego, ale daje własność, której nie da się kupić za żadną cenę: system otwarty za dziesięć lat nadal będzie działał, bo nie ma w nim niczego, co mogłoby przestać być utrzymywane.
Wewnątrz silnika obowiązuje konsekwentnie układ SI: metry, niutony, paskale, kilogramy. Przeliczenia na milimetry i megapaskale wykonuje warstwa interfejsu. W kodzie solvera nie ma ani jednego mnożnika 10⁻³ czy 10⁶, a to najczęstsze źródło błędów rzędu tysiąca w wynikach obliczeń wytrzymałościowych.
Każda analiza deklarowana jako działająca ma odniesienie w rozwiązaniu zamkniętym i test w zestawie weryfikacyjnym; analizy niegotowe są widoczne w katalogu jako zapowiedzi i nie da się ich wybrać. Wynik niesie ze sobą kontrolę równowagi globalnej, a raport — wykaz zastrzeżeń. Rozdział 12 pokazuje, jak ta zasada jest egzekwowana.
Punktem wyjścia jest układ równań opisujący równowagę ciała odkształcalnego. W obszarze Ω o brzegu ∂Ω poszukujemy pola przemieszczeń u(x) spełniającego równanie równowagi wewnętrznej wraz ze związkami geometrycznymi i fizycznymi:
do których dochodzą warunki brzegowe: kinematyczny u = ū na części brzegu Γu oraz statyczny σ·n = t̄ na pozostałej części Γt. Rozwiązanie ścisłe tego układu istnieje tylko dla nielicznych, prostych geometrii — i to jest cała przyczyna, dla której metoda elementów skończonych w ogóle powstała.
Zamiast żądać, by równanie (2.1) było spełnione w każdym punkcie obszaru, żądamy, by było spełnione w sensie całkowym — po przemnożeniu przez dowolne pole przemieszczeń wirtualnych δu znikające tam, gdzie przemieszczenia są narzucone. Po scałkowaniu przez części otrzymujemy zasadę prac wirtualnych:
Praca sił wewnętrznych na odkształceniach wirtualnych równa się pracy obciążeń zewnętrznych na przemieszczeniach wirtualnych. Postać ta ma dwie zalety rozstrzygające o jej użyteczności: wymaga od pola przemieszczeń różniczkowalności tylko pierwszego rzędu — zamiast drugiego, jak (2.1) — oraz włącza statyczne warunki brzegowe w sposób naturalny, jako składnik prawej strony, zamiast nakładać je osobno.
Obszar dzielimy na elementy, a pole przemieszczeń wewnątrz elementu przybliżamy wartościami w jego węzłach, interpolowanymi funkcjami kształtu Na:
Odkształcenie, jako pochodna przemieszczenia, wyraża się przez pochodne funkcji kształtu. Macierz łącząca przemieszczenia węzłowe z odkształceniem oznaczamy tradycyjnie B; dla elementu bryłowego o n węzłach ma ona wymiar 6 × 3n:
Podstawienie (2.5) i (2.6) do zasady prac wirtualnych (2.4) i wykorzystanie dowolności δu prowadzi wprost do układu równań algebraicznych, w którym macierz sztywności elementu i wektor sił węzłowych mają postać:
Agregacja polega na dodaniu wyrazów macierzy elementu do odpowiadających im pozycji macierzy
globalnej, wskazanych przez numery stopni swobody elementu. W kodzie odpowiada temu klasa
Assembler, a same całki (2.7) — funkcja elementSystem.
Całek (2.7) nie liczy się analitycznie. Zastępuje je sumowanie wartości podcałkowej w wybranych punktach z odpowiednimi wagami — kwadratura Gaussa:
Wyznacznik jakobianu przekształcenia z układu naturalnego (ξ, η, ζ) do globalnego (x, y, z) pełni tu podwójną rolę: jest miarą objętości i jednocześnie sprawdzianem poprawności siatki. Wartość ujemna oznacza element wywrócony — błąd w kolejności węzłów, a nie w solverze — i dlatego silnik przerywa wtedy obliczenia z jawnym komunikatem, zamiast liczyć dalej na ujemnej objętości.
Liczba punktów całkowania nie jest kwestią gustu. Musi wystarczyć do ścisłego scałkowania iloczynu BᵀDB, ale nadmiar kosztuje czas proporcjonalnie. Dla tetraedru o prostych krawędziach jakobian jest stały, więc dla TET4 wystarcza jeden punkt, a dla TET10 — reguła czteropunktowa dokładna dla wielomianów drugiego stopnia. HEX8 całkowany jest pełną kwadraturą 2×2×2.
System udostępnia trzy elementy bryłowe. Różnice między nimi nie są kosmetyczne — decydują o tym, czy wynik będzie miał błąd promila, czy stu procent, przy tej samej liczbie niewiadomych. Rozdział 13 pokazuje to na liczbach.
Funkcje kształtu są współrzędnymi barycentrycznymi, więc pole przemieszczeń jest liniowe, a odkształcenie — stałe w całej objętości elementu. Element o stałym odkształceniu nie potrafi opisać zginania inaczej niż schodkowo: zgięcie belki odwzorowuje ciągiem elementów o różnych, ale w środku jednorodnych stanach naprężenia.
Skutek jest drastyczny. W przykładzie z punktu 13.2 ten sam wspornik policzony elementami TET4 daje pierwszą częstość drgań własnych zawyżoną o 133%, podczas gdy TET10 na tej samej siatce myli się o 0,54%. TET4 pozostaje w systemie, bo jest tani i wystarcza tam, gdzie dominuje rozciąganie lub ściskanie — nie dlatego, że jest dobrym wyborem domyślnym.
Dziesięć węzłów — cztery narożne i sześć na środkach krawędzi — daje kwadratowe pole przemieszczeń i liniowe odkształcenie. To wystarcza, by element zginał się poprawnie zamiast ścinać pasożytniczo. Numeracja węzłów jest zgodna z kartami C3D10 pakietów Abaqus i CalculiX oraz CTETRA formatu NASTRAN, dzięki czemu porównanie wyniku z zewnętrznym solverem nie wymaga przestawiania indeksów.
TET10 jest w systemie elementem zalecanym i domyślnym. Kosztuje około 2,5 raza więcej stopni swobody na tej samej siatce, ale — jak pokazuje tabela 13.2 — daje wynik dokładniejszy o dwa rzędy wielkości, więc jest tańszy przy zadanej dokładności, nie droższy.
Ośmiowęzłowy sześcian o trójliniowych funkcjach kształtu jest tani i regularny, ale w czystej postaci cierpi na dwie dolegliwości, z których każda potrafi zawyżyć sztywność o rzędy wielkości.
Przy liczbie Poissona zbliżonej do 0,5 materiał staje się praktycznie nieściśliwy. Element o niskim rzędzie interpolacji nie potrafi jednocześnie spełnić warunku stałej objętości i odwzorować deformacji postaciowej, więc usztywnia się drastycznie. Środkiem zaradczym jest sformułowanie B-bar: część objętościowa macierzy B zostaje zastąpiona jej wartością uśrednioną po elemencie.
Element trójliniowy nie ma w polu przemieszczeń członu kwadratowego, więc poddany zginaniu nie zgina się, tylko ścina. Efekt narasta z wydłużeniem elementu — przy stosunku boków 4:1 potrafi zawyżyć sztywność kilkukrotnie. Lekarstwem są mody niezgodne Wilsona: trzy dodatkowe funkcje bąblowe dokładające brakującą krzywiznę.
Dziewięć wewnętrznych stopni swobody α (trzy funkcje × trzy kierunki) nie trafia do układu globalnego — usuwa je kondensacja statyczna wykonana wewnątrz elementu:
Element pozostaje więc ośmiowęzłowy z punktu widzenia solvera. Istotny szczegół: pochodne funkcji bąblowych przelicza się na współrzędne globalne z poprawką Taylora–Beresforda–Wilsona, biorąc jakobian ze środka elementu i skalując stosunkiem wyznaczników. Bez tej poprawki element przestaje przechodzić test łaty na siatce zniekształconej — czyli dokładnie wtedy, kiedy najbardziej go potrzebujemy.
Macierz konstytutywna materiału izotropowego wyraża się przez stałe Lamégo, wyliczane z modułu Younga i liczby Poissona:
Konwencja Voigta przyjęta w systemie: naprężenie [σxx, σyy, σzz, σxy, σyz, σzx] jako składowe tensorowe bez mnożnika, odkształcenie [εxx, εyy, εzz, γxy, γyz, γzx] ze ścinaniem inżynierskim γ = 2ε. Przy takim zestawieniu iloczyn σ·ε jest poprawną gęstością pracy, a macierz D nie wymaga dodatkowego skalowania.
Model plastyczności opiera się na warunku Hubera–Misesa z wzmocnieniem mieszanym — izotropowym i kinematycznym. Powierzchnia plastyczności ma postać:
gdzie s jest dewiatorem naprężenia, α — naprężeniem resztkowym wzmocnienia kinematycznego, a εᵖ — zastępczym odkształceniem plastycznym. Całkowanie prowadzi się algorytmem powrotu promienistego: najpierw liczymy naprężenie próbne przy zamrożonym odkształceniu plastycznym, a jeżeli wypada poza powierzchnią plastyczności — rzutujemy je z powrotem wzdłuż kierunku dewiatora. Dla wzmocnienia liniowego mnożnik plastyczny wychodzi wprost, bez iteracji:
Kluczem do zbieżności całej analizy nieliniowej jest moduł styczny spójny z algorytmem (Simo–Taylor), a nie ciągły moduł sprężysto-plastyczny. Tylko on zachowuje kwadratową zbieżność metody Newtona:
System zawiera bibliotekę 129 gatunków w czternastu grupach — od stali konstrukcyjnych, przez stopy lekkie i tworzywa, po ceramikę i materiały budowlane. Każdy wpis niesie komplet dziewięciu właściwości, których używa solver, wymienionych w tabeli 4.1.
| Właściwość | Jednostka | Gdzie decyduje o wyniku |
|---|---|---|
| Moduł Younga E | GPa | macierz konstytutywna — każda analiza |
| Liczba Poissona ν | — | macierz konstytutywna; przy ν → 0,5 blokada objętościowa |
| Gęstość ρ | kg/m³ | drgania własne, dynamika, ciężar własny |
| Granica plastyczności Re | MPa | współczynnik bezpieczeństwa, uplastycznienie |
| Wytrzymałość Rm | MPa | poprawki na naprężenie średnie w analizie zmęczeniowej |
| Moduł wzmocnienia H | MPa | nachylenie krzywej za granicą plastyczności |
| Rozszerzalność α | 10⁻⁶/K | naprężenia termiczne, sprzężenie termomechaniczne |
| Przewodność λ | W/(m·K) | przewodnictwo ustalone i nieustalone |
| Ciepło właściwe c | J/(kg·K) | pojemność cieplna w zadaniach nieustalonych |
Wartości są typowe dla gatunku, zebrane z norm przedmiotowych i kart katalogowych — nie są wynikami badań konkretnej partii wsadowej. Rzeczywista granica plastyczności wyrobu bywa wyższa od normowej o kilkanaście procent, a właściwości tworzyw zależą od wilgotności, prędkości odkształcania i kierunku wtrysku silniej niż od samego gatunku. Do obliczeń trafiających do dokumentacji dane bierze się z karty dostawcy — i system mówi to wprost w panelu materiałów, zamiast pozwalać traktować bibliotekę jako źródło danych odbiorowych.
Osobną wartość ma ostrzeżenie przypisane do materiału. Żeliwo szare, szkło, ceramika i beton są materiałami kruchymi o różnej wytrzymałości na rozciąganie i ściskanie — hipoteza Hubera-Misesa, którą system pokazuje jako pole domyślne, jest dla nich kryterium niewłaściwym i należy je oceniać naprężeniem głównym σ₁. Kompozyty i drewno są ortotropowe, więc podane wartości dotyczą kierunku włókien, a wynik izotropowy jest miarodajny wyłącznie dla obciążenia zgodnego z tym kierunkiem. Te uwagi są częścią karty materiału i przechodzą do raportu z analizy — bo błąd interpretacji poprawnie policzonego wyniku jest w praktyce częstszy od błędu obliczeń.
Reologia opisuje zachowanie materiału rozłożone w czasie: pełzanie, czyli narastanie odkształcenia pod stałym obciążeniem, oraz relaksację, czyli spadek naprężenia przy stałym odkształceniu. Dla stali w temperaturze pokojowej zjawiska te są pomijalne. Nie znaczy to jednak, że dotyczą wyłącznie przypadków skrajnych — przeciwnie:
System zbiera te dane w dwóch modelach:
W prawie Nortona wykładnik n mieści się zwykle w przedziale od 3 do 8 i mówi, jak silnie prędkość pełzania rośnie z naprężeniem, a energia aktywacji Q — jak szybko rośnie ona z temperaturą. W szeregu Prony’ego udział g∞ opisuje część sprężystą, która nie relaksuje; udziały g∞ i wszystkie gᵢ muszą sumować się do jedności, co system sprawdza przy wprowadzaniu danych.
Użytkownik może dodać własny materiał, a formularz wymaga deklaracji modelu reologicznego. Odpowiedź „sprężysty, bez efektów reologicznych" jest dopuszczalna i dla stali konstrukcyjnej poprawna — musi być jednak wybrana świadomie, a nie pominięta przez przeoczenie, bo dla tworzywa, cynku czy ołowiu byłaby nieprawdziwa.
Macierz sztywności zadania o stu tysiącach stopni swobody miałaby w postaci pełnej 10¹⁰ wyrazów, czyli 80 GB. Niezerowych jest wśród nich rzędu kilku milionów. System przechowuje więc wyłącznie je, w formacie CSR (compressed sparse row): tablica wartości, tablica indeksów kolumn i tablica wskaźników początków wierszy.
Istotne jest rozdzielenie wzorca od wartości. Wzorzec — który wyraz jest niezerowy — zależy wyłącznie od topologii siatki, więc liczy się raz. Wartości zmieniają się w każdej iteracji Newtona i w każdym kroku czasowym. Dzięki temu rozdziałowi zadanie nieliniowe nie płaci za analizę struktury przy każdym kroku obciążenia, a sesja licząca kilka analiz na jednej siatce buduje wzorzec tylko przy pierwszej.
Macierz sztywności jest symetryczna i — po nałożeniu wystarczających więzów — dodatnio określona. Rozkład Choleskiego byłby dla niej naturalny, ale wymaga pierwiastkowania, a przy zadaniach wyboczeniowych macierz bywa nieokreślona. Stosujemy więc rozkład LDLᵀ:
Zaimplementowano wariant „up-looking" Timothy'ego Davisa, w którym kolumna czynnika powstaje z przejścia w górę po drzewie eliminacji — dzięki czemu wzorzec niezerowych wyrazów wychodzi sam, bez wyszukiwania. Rozwiązanie układu sprowadza się do trzech podstawień:
Faktoryzacja daje przy okazji informację o wartości dla całego systemu: liczbę ujemnych wyrazów przekątnej D. Zgodnie z twierdzeniem Sylvestera o bezwładności jest to liczba wartości własnych leżących poniżej przesunięcia — używamy jej do sprawdzenia, czy solver zagadnienia własnego nie zgubił żadnej postaci. Odpowiadający temu test jest w zestawie weryfikacyjnym pod nazwą „Bezwładność macierzy".
Podczas faktoryzacji w czynniku L pojawiają się wyrazy niezerowe tam, gdzie macierz pierwotna miała zera. To zjawisko — wypełnienie — decyduje o zapotrzebowaniu na pamięć i czas. Jego rozmiar zależy dramatycznie od kolejności eliminacji niewiadomych, a dobór tej kolejności jest osobnym zagadnieniem kombinatorycznym. System stosuje algorytm przybliżonego minimalnego stopnia, a skuteczność przenumerowania jest przedmiotem testu „Wpływ przenumerowania na wypełnienie".
Powyżej progu 60 000 stopni swobody czynnik LDLᵀ przestaje mieścić się wygodnie w pamięci przeglądarki. System przechodzi wtedy automatycznie na metodę gradientów sprzężonych z prekondycjonerem niepełnego rozkładu Choleskiego IC(0), która nie tworzy wypełnienia wcale — kosztem tego, że rozwiązanie jest przybliżone z zadaną tolerancją i że każdy nowy wektor prawej strony wymaga pełnej iteracji od nowa.
Próg przełączenia dobrano z pomiarów, nie z teorii. Zgodność obu dróg — bezpośredniej i iteracyjnej — jest sprawdzana testem „Gradienty sprzężone wobec solvera bezpośredniego".
Wybór między rozkładem bezpośrednim a metodą iteracyjną opisuje się zwykle jednym progiem liczby niewiadomych. Jest to uproszczenie, które w tym systemie okazało się mylące, i warto powiedzieć dlaczego, bo przyczyna jest ogólna.
Rozkład LDLᵀ ma dwa koszty o zupełnie różnym charakterze: jednorazową faktoryzację i wielokrotne podstawienie wsteczne. Gradienty sprzężone mają jeden koszt, ponoszony za każdym razem od nowa. Próg dobrany dla statyki — gdzie układ rozwiązuje się raz — jest więc dobrany dla najgorszego możliwego przypadku faktoryzacji i najlepszego dla metody iteracyjnej.
| Etap | Rozkład LDLᵀ | Gradienty sprzężone |
|---|---|---|
| Przygotowanie | 280 s, 784 MB wypełnienia | 0,5 s (prekondycjoner IC(0)) |
| Jedno rozwiązanie układu | 0,53 s | 8,5 s (212 iteracji) |
| Statyka — jedno rozwiązanie | 281 s | 9 s |
| Drgania własne — około 180 rozwiązań | 376 s | 1530 s |
Konsekwencji praktycznej nie da się jednak wyciągnąć wprost, bo drugą granicą jest pamięć. Wypełnienie rośnie dla siatek przestrzennych mniej więcej jak n4/3: przy stu tysiącach stopni swobody to już półtora giga, czego karta przeglądarki nie udźwignie. Przy takich siatkach żadna droga nie jest dobra i jedyną uczciwą odpowiedzią jest powiedzenie tego użytkownikowi.
Dlatego zadanie własne powyżej progu robi trzy rzeczy, których nie robiło wcześniej: mierzy czas pierwszego rozwiązania i podaje z niego rzut czasu całości, zgłasza kolejne rozwiązania w logu, oraz ma budżet czasu, po którym zwraca częstości z ostatniej ukończonej iteracji. Wynik częściowy jest przy tym pełnowartościową informacją, a nie awaryjnym zastępnikiem: zbieżność iteracji podprzestrzeni jest jednostronna, więc częstości przerwanego zadania są oszacowaniem od góry, a nie liczbą o nieznanym błędzie. Bez tego zadanie na gęstej siatce wyglądało po prostu na zawieszony program — i było to jedyne, co użytkownik mógł o nim wywnioskować.
Drgania własne i wyboczenie liniowe sprowadzają się do uogólnionego zagadnienia własnego o tej samej postaci algebraicznej, choć o zupełnie różnym znaczeniu fizycznym:
Interesuje nas kilka najniższych par własnych z zadania o dziesiątkach tysięcy niewiadomych, więc metody wyznaczające całe widmo nie wchodzą w grę. Zastosowano iterację podprzestrzeni w ujęciu Bathego: zamiast szukać każdej postaci osobno, prowadzimy blok wektorów naraz. Mnożenie przez K⁻¹ tłumi składowe o wysokich wartościach własnych, a rzut Rayleigha–Ritza na bieżącą podprzestrzeń wydobywa z niej najlepsze przybliżenia:
Rozmiar bloku jest większy od liczby szukanych postaci — zgodnie z zaleceniem Bathego min(2q, q+8). Zapas jest tu konieczny, nie zalecany: bez niego iteracja gubi postacie o bliskich sobie częstościach, a takie występują wszędzie tam, gdzie konstrukcja ma symetrię. Przykład z punktu 13.3 pokazuje to wprost — pręt o przekroju kwadratowym ma dwie niemal identyczne postacie wyboczenia, po jednej na każdą płaszczyznę.
Przy wyboczeniu macierz po prawej stronie bywa nieokreślona, co uniemożliwia rozkład Choleskiego w małym zadaniu rzutowanym. Rozwiązanie jest proste: zamieniamy strony, bo macierzą dodatnio określoną jest zawsze macierz sztywności, i odwracamy otrzymane wartości własne.
Przy dużych przemieszczeniach liniowa miara odkształcenia przestaje być poprawna — obrót ciała sztywnego generowałby w niej naprężenia. Stosujemy sformułowanie Lagrange’a całkowite: odniesieniem pozostaje konfiguracja początkowa, miarą odkształcenia jest tensor Greena–Lagrange’a, a miarą naprężenia — drugi tensor Pioli–Kirchhoffa.
Sprawdzianem poprawności jest test niezmienniczości: pręt obrócony o 30° jako ciało sztywne musi wykazywać zerowe naprężenie. Zestaw weryfikacyjny podaje przy tym liczbę, która najlepiej uzasadnia potrzebę tego sformułowania — miara małych odkształceń daje przy samym obrocie 21,6 GPa naprężenia zastępczego, sformułowanie nieliniowe: 2·10⁻¹⁶ GPa.
Zadanie nieliniowe rozwiązujemy przyrostowo, szukając w każdym kroku obciążenia zera residuum:
Sterowanie obciążeniem uzupełniono automatycznym skracaniem kroku: gdy iteracja się nie zbiega, przyrost jest połowiony, a gdy zbiega w trzech iteracjach lub mniej — wydłużany. Jest to prostsze i pewniejsze niż zgadywanie z góry, ile kroków wystarczy. Norma residuum odnoszona jest do obciążenia, więc jest bezwymiarowa i niezależna od jednostek.
Sterowanie obciążeniem nie przejdzie przez punkt graniczny: w chwili, gdy konstrukcja traci nośność, macierz styczna staje się osobliwa i metoda się rozbiega. Zadania z przeskokiem — płaskie łuki, powłoki, zatrzaski — wymagają potraktowania mnożnika obciążenia jako dodatkowej niewiadomej i uzupełnienia układu równaniem więzu:
Kontakt realizowany jest metodą funkcji kary z warunkiem nieprzenikania i tarciem Coulomba. Stan tarcia — przyczepność albo poślizg — zapamiętywany jest dopiero po zbieżności kroku, ponieważ w trakcie iteracji węzeł może wielokrotnie przechodzić między tymi stanami. Tarcie działa wiarygodnie w zakresie przyczepności; poślizg właściwy bywa niezbieżny, co jest wymienione w wykazie ograniczeń.
Równanie ruchu układu dyskretnego ma postać, do której sprowadzają się wszystkie analizy dynamiczne — różnią się one wyłącznie sposobem jego rozwiązania:
System stosuje macierz mas skupionych, otrzymaną metodą HRZ: proporcje bierze się z przekątnej macierzy spójnej, a następnie skaluje tak, by ich suma dała dokładnie masę elementu. Macierz diagonalna pozwala odwracać masę bez rozwiązywania układu, co jest warunkiem opłacalności całkowania jawnego. Poprawność masy jest sprawdzana wprost — test porównuje masę modelu z iloczynem ρ·V i wymaga zgodności do 10⁻⁹.
Dla wymuszenia sinusoidalnego o ustalonej częstości rozwiązanie ustalone znajdujemy w dziedzinie zespolonej, bez całkowania po czasie:
Sprawdzianem jest wzmocnienie rezonansowe: przy tłumieniu ζ amplituda w rezonansie powinna być 1/(2ζ) razy większa od statycznej. Zestaw weryfikacyjny odtwarza tę zależność z błędem 0,001%.
Przebieg czasowy liczony jest metodą Newmarka z parametrami γ = ½ i β = ¼ (reguła średniego przyspieszenia). Schemat jest bezwarunkowo stabilny i nie wprowadza tłumienia sztucznego, więc krok czasowy dobiera się do dokładności, a nie do stabilności.
Test sprawdza klasyczny wynik teoretyczny: nagłe przyłożenie obciążenia stałego daje ugięcie dwukrotnie większe od statycznego. Solver odtwarza wartość 1,99 wobec 2,00 z teorii.
Przy zjawiskach szybkozmiennych opłaca się schemat różnic centralnych: nie wymaga rozwiązywania układu równań, bo macierz mas jest diagonalna. Cena to stabilność warunkowa — krok czasowy musi być mniejszy od czasu przebiegu fali przez najmniejszy element (warunek Couranta):
Silnik pilnuje tego warunku sam i raportuje zarówno krok użyty, jak i graniczny. Zgodność z całkowaniem niejawnym potwierdzono z dokładnością 0,6%.
Dla obciążeń sejsmicznych i losowych całkowanie po czasie jest niepotrzebnie kosztowne. Metoda spektrum odpowiedzi składa wkłady postaci własnych, a analiza drgań losowych operuje gęstością widmową mocy:
Wartość skuteczna z analizy drgań losowych zgadza się ze wzorem Milesa z błędem 0,6%.
Przewodnictwo ciepła, akustyka i filtracja prowadzą do formalnie tego samego zadania: jednego stopnia swobody w węźle i macierzy typu „sztywności" zbudowanej z gradientów funkcji kształtu. System rozwiązuje je wspólnym modułem pola skalarnego, różnicując wyłącznie znaczenie współczynników i warunków brzegowych.
Przewodnictwo nieustalone całkowane jest wstecznym Eulerem — schematem bezwarunkowo stabilnym, co przy zadaniach cieplnych jest ważniejsze od rzędu dokładności, bo stałe czasowe stygnięcia bywają o rzędy wielkości różne w jednym modelu.
Odkształcenie termiczne odejmowane jest od całkowitego przed wyznaczeniem naprężenia:
Dla pręta całkowicie skrępowanego daje to klasyczny wynik σ = −E·α·ΔT, który system odtwarza dokładnie — patrz przykład 13.4. Sprzężenie termomechaniczne jest sekwencyjne: najpierw rozwiązywane jest pole temperatury, potem zadanie mechaniczne. Odkształcenie nie wpływa zwrotnie na temperaturę, co dla zdecydowanej większości zadań inżynierskich jest założeniem poprawnym.
Przepływ tworzywa w cienkościennym gnieździe formy opisuje przybliżenie Hele-Shawa. Wynika ono z pominięcia bezwładności wobec sił lepkich — liczba Reynoldsa dla stopionego polimeru jest rzędu 10⁻³ — oraz z założenia, że grubość ścianki jest znacznie mniejsza od pozostałych wymiarów wypraski. Zadanie sprowadza się wtedy do jednego równania na ciśnienie, a więc do tego samego modułu pola skalarnego:
Wykładnik przy grubości wymaga komentarza, bo w piśmiennictwie występuje jako h³. Sześcian jest właściwy dla sformułowania powierzchniowego, w którym wypraskę odwzorowuje siatka dwuwymiarowa, a grubość wchodzi wyłącznie przez współczynnik: mieści on wtedy w sobie zarówno paraboliczny profil prędkości po grubości, jak i pole przekroju. W systemie równanie rozwiązywane jest na siatce przestrzennej, gdzie całkowanie objętościowe uwzględnia przekrój kanału samo z siebie; pozostawienie h³ liczyłoby grubość dwa razy. Pomyłka nie jest niewinna — dla ścianki 2,5 mm zawyżała ciśnienie wtrysku o dwa rzędy wielkości.
Czoło strugi wyznaczane jest z rozkładu ciśnienia: węzły porządkuje się według malejącego ciśnienia i napełnia w tej kolejności objętościami do nich przypisanymi, aż do wyczerpania strumienia. Daje to czas dopłynięcia w każdym węźle, a z niego dwie informacje, dla których tę analizę się w ogóle wykonuje. Linie łączenia powstają tam, gdzie sąsiednie elementy mają kierunki przepływu rozbieżne o więcej niż 120° — czyli tam, gdzie dwie strugi spotykają się czołami. Kryterium oparte na samym czasie dopłynięcia jest zwodnicze: wskazuje każdą izochronę, a więc setki miejsc także tam, gdzie wlew jest jeden i żadna linia łączenia powstać nie może. Miejsca zamknięcia powietrza to lokalne maksima czasu wypełnienia — obszary napełniane jako ostatnie, z których powietrze nie ma już którędy ujść.
Szybkość ścinania przyjmowana jest jako wartość charakterystyczna γ̇ ≈ 6·v/h, wyznaczona z prędkości średniej przy wlocie, a nie z rozwiązania. Sprzężenie lepkości z polem prędkości i temperatury należy do analizy pełnej; w tej wersji lepkość jest jedną liczbą na całą wypraskę. Rozkład ciśnienia wychodzi zatem ciśnieniem w samym gnieździe — bez oporu przewężki, kanałów doprowadzających i dyszy, które w rzeczywistej wtryskarce odpowiadają za większą część wskazania manometru. Do porównania wariantów rozmieszczenia wlewów jest to wystarczające, do doboru nastaw maszyny — nie.
Guma i elastomery wymagają opisu innego rodzaju niż metale. Naprężenie nie jest liniową funkcją odkształcenia w żadnym zakresie — wyprowadza się je z energii odkształcenia sprężystego. System używa modelu neo-Hooke’a w postaci ściśliwej:
Wybór tej postaci — zamiast częstszej w piśmiennictwie, z członem objętościowym (κ/2)(J−1)² — nie jest obojętny. Postać powyższa daje przy małych odkształceniach dokładnie stałe Lamégo λ i μ, więc ten sam materiał policzony liniowo i hipersprężyście zgadza się w zakresie, w którym oba modele obowiązują. Postać z (J−1)² różni się w module objętościowym o 2μ/3; rozjazd byłby widoczny już przy odkształceniu rzędu procenta i nie dałoby się rozstrzygnąć, czy bierze się z modelu, czy z błędu implementacji.
Ograniczeniem jest ściśliwość. Guma ma ν sięgające 0,4995, a elementy przemieszczeniowe blokują się objętościowo w tym zakresie: sztywność objętościowa przewyższa postaciową o cztery rzędy wielkości i zadanie staje się źle uwarunkowane. Sformułowanie mieszane, w którym ciśnienie jest osobną niewiadomą, nie zostało zaimplementowane — analiza sprawdza więc liczbę Poissona i mówi wprost, że wynik będzie za sztywny, zamiast zwrócić go bez komentarza. Drugim ograniczeniem jest sam model: neo-Hooke odtwarza krzywą rozciągania gumy do mniej więcej 100% wydłużenia, wyżej rzeczywisty materiał usztywnia się gwałtownie i potrzebny jest model Ogdena albo Arrudy-Boyce’a.
Napięcia wstępnego nie da się zadać jako obciążenia. Siła w śrubie jest wynikiem tego, że trzpień skrócono dokręceniem, a złącze się temu opiera; jej wartość zależy od podatności trzpienia i ściskanych części:
Przyłożenie pary sił ±F do łba i nakrętki jest błędem podwójnym: pomija sprężystość złącza, a przy późniejszym obciążeniu roboczym siła w śrubie nie zmienia się wcale — gdy tymczasem cały sens obliczania złącza polega na tym, żeby zobaczyć, jak ona rośnie i kiedy złącze się rozwiera.
Rozwiązanie idzie przez odkształcenie wstępne. Elementom trzpienia narzucany jest jednoosiowy skurcz ε₀ = −e·(n ⊗ n) — wyłącznie wzdłuż osi śruby, bez zmiany wymiarów poprzecznych, inaczej niż przy odkształceniu cieplnym, które jest objętościowe. Daje to zastępczy wektor obciążenia ∫Bᵀ·D·ε₀·dV, a naprężenie w trzpieniu wynosi σ = D·(ε(u) − ε₀). Zadanie jest liniowe względem e, więc jedno przejście próbne wystarcza do wyznaczenia współczynnika skali; drugie daje siłę żądaną dokładnie, a nie iteracyjnie. Współczynnik podatności złącza Φ wyznaczany jest z różnicy dwóch policzonych stanów, a nie ze wzorów tablicowych na sztywność ściskanego stożka — model bryłowy zna rzeczywisty rozkład sztywności, który te wzory tylko przybliżają.
Przebieg liniowy rozwiązuje równanie ruchu przy macierzy K stałej przez cały czas trwania zjawiska. Założenie to przestaje obowiązywać dokładnie tam, gdzie analiza dynamiczna jest najbardziej potrzebna: przy uderzeniu z uplastycznieniem, przy dużych ugięciach z usztywnieniem błonowym i przy zamykaniu się szczeliny kontaktowej. W każdym z tych przypadków wynik liniowy jest gładki, wiarygodnie wyglądający i nieprawdziwy — i nic tego nie sygnalizuje.
Całkowanie idzie schematem Newmarka w postaci przemieszczeniowej, z iteracją Newtona-Raphsona wewnątrz każdego kroku czasowego:
Macierz tłumienia liczona jest na sztywności początkowej, nie stycznej. Styczna przy uplastycznieniu maleje, a przy przejściu przez punkt graniczny zmienia znak; tłumienie policzone na niej rosłoby i malało razem z nią, a w skrajnym przypadku dostarczałoby energii zamiast ją rozpraszać. Tłumienie konstrukcyjne jest własnością konstrukcji, nie jej chwilowego stanu.
Naprężenie w karbie zależy od promienia zaokrąglenia, a promień bywa pięćdziesiąt razy mniejszy od gabarytu części. Siatka odwzorowująca karb wiernie ma w całej części tyle elementów, że zadania nie da się policzyć; siatka, którą da się policzyć, pokazuje w karbie wartość zależną od tego, gdzie akurat wypadł węzeł.
Podmodelowanie rozdziela te wymagania. Model globalny liczy rozkład sztywności zgrubną siatką, podmodel obejmuje sam karb z siatką dowolnie gęstą, a warunki na powierzchniach cięcia bierze z rozwiązania globalnego. Podstawą jest zasada św. Wenanta: rozkład naprężenia z dala od zaburzenia zależy tylko od wypadkowej. Przemieszczenia przenoszone są funkcjami kształtu elementu globalnego, po wyznaczeniu współrzędnych naturalnych iteracją Newtona na odwzorowaniu x(ξ) = Σ Nₐ(ξ)·xₐ — interpolacja odległościowa, choć prostsza, nie odtwarza nawet ruchu sztywnego.
Warunek „cięcie dostatecznie daleko od karbu" jest jedyną rzeczą, którą musi rozstrzygnąć człowiek, i decyduje o wiarygodności całości. Dlatego analiza sprawdza go po swojemu: porównuje naprężenie na powierzchniach cięcia policzone w obu modelach. Bez tej kontroli podmodelowanie daje wynik gładki i dowolnie zły. Osobną pułapką jest wykrywanie samej powierzchni cięcia: rozstrzyga o niej ściana, nie pojedynczy węzeł, bo węzeł na obwodzie cięcia leży jednocześnie na wolnej powierzchni bocznej i klasyfikowany osobno wypadłby z cięcia, zostawiając podmodel podparty tylko w środku ściany.
Model złożenia rośnie szybciej niż moc obliczeniowa, a pytanie bywa niezmienne: jakie są pierwsze częstości i jak konstrukcja odpowiada na wymuszenie. Odpowiedź dynamiczna leży w przestrzeni rozpiętej przez kilka najniższych postaci drgań — reszta stopni swobody potrzebna jest do tego, żeby te postacie policzyć, nie do tego, żeby je opisać. Redukcja wymienia więc miliony niewiadomych na kilkadziesiąt współrzędnych uogólnionych o tej samej treści fizycznej.
Punktem wyjścia jest podział stopni swobody na brzegowe (te, którymi podstruktura łączy się z resztą) i wewnętrzne. Kondensacja Guyana wyraża wnętrze przez brzeg:
Jest to zależność ścisła, jeżeli obciążenie przyłożone jest wyłącznie do brzegu: wnętrze nie ma wtedy własnej prawej strony, więc równanie równowagi wnętrza daje dokładnie wzór (9.7). Ta własność jest sprawdzalna bez tolerancji i w zestawie weryfikacyjnym sprawdzana właśnie tak — przemieszczenia modelu zredukowanego zgadzają się z pełnym z błędem 3·10⁻¹⁰ %, czyli na poziomie zaokrągleń podwójnej precyzji.
W dynamice kondensacja jest natomiast przybliżeniem, i to jednostronnym: pomija bezwładność wnętrza, więc konstrukcja wychodzi za sztywna, a częstości zawyżone. Metoda, która myli się zawsze w tę samą stronę, jest w obliczeniach wytrzymałościowych gorsza od mylącej się losowo — nie daje się skorygować zapasem, bo błąd rośnie z numerem postaci w sposób zależny od zadania.
Craig i Bampton dokładają do postaci więzowych postacie własne o zamocowanym brzegu — rozwiązania (K_ww − ω²M_ww)φ = 0 przy u_b = 0. Baza redukcji ma wtedy dwa rodzaje kolumn o wyraźnie rozdzielonych zadaniach:
Pierwszy blok przenosi ruch brzegu dokładnie tak, jak robi to statyka. Drugi dokłada zdolność wnętrza do drgania niezależnie od brzegu — czyli dokładnie to, czego kondensacji brakuje.
Zbieżność jest jednostronna i tym razem jest to zaletą: baza jest podprzestrzenią przestrzeni pełnej, a iloraz Rayleigha na podprzestrzeni nie schodzi poniżej minimum, więc każda częstość zredukowana jest górnym oszacowaniem dokładnej. Wynik ma więc znany kierunek błędu, a dokładanie postaci wnętrza może go tylko zmniejszyć. Test sprawdza to wprost: częstość, która wyszłaby niższa od dokładnej, oznaczałaby błąd w rzutowaniu, a nie lepsze przybliżenie.
| Postaci wnętrza | Współrzędnych | Redukcja | Błąd częstości 1 | Błąd częstości 3 |
|---|---|---|---|---|
| 0 (sam Guyan) | 48 | 90,9% | 0,49% | 44,85% |
| 2 | 50 | 90,5% | 0,04% | 0,39% |
| 5 | 53 | 90,0% | 0,008% | 0,12% |
| 20 | 68 | 87,1% | 0,0002% | 0,004% |
Ograniczeniem obecnej wersji jest to, że redukowana jest jedna podstruktura. Superelement powstaje i jego macierze są zwracane w całości, ale sklejenie dwóch takich elementów po wspólnym brzegu wymaga więzów wielopunktowych, których model jeszcze nie ma. Analiza jest więc redukcją struktury, a nie montażem złożenia z klocków — i tak jest opisana w interfejsie, zamiast obiecywać drugie.
Analizy zmęczeniowe nie rozwiązują nowego zadania brzegowego — interpretują policzone wcześniej pole naprężeń. Ich wynikiem jest liczba cykli, a nie rozkład przestrzenny, co ma swoje odzwierciedlenie w interfejsie: dla tych analiz moduł mapy wyników jest wyłączany, bo nie ma czego kolorować.
Wpływ naprężenia średniego uwzględniany jest jedną z trzech poprawek — Goodmana, Gerbera albo Soderberga — a sumowanie uszkodzeń regułą Palmgrena–Minera:
Równania (10.3) nie da się rozwiązać względem N w postaci zamkniętej. System stosuje metodę połowienia przedziału, doprowadzając resztę równania poniżej 10⁻⁶.
Całkowanie prowadzi się od wady początkowej do wymiaru krytycznego wyznaczonego odpornością na pękanie. Dla wykładnika m ≠ 2 istnieje rozwiązanie zamknięte, z którym wynik numeryczny zgadza się co do cyfry.
System dzieli się na trzy warstwy o wyraźnie rozdzielonych zadaniach. Podział ten nie jest ozdobnikiem — to on pozwala uruchomić ten sam solver w przeglądarce, w wątku roboczym i w środowisku Node bez jednej linii kodu warunkowego.
| Warstwa | Zawartość | Zależności |
|---|---|---|
| Silnik TOMes® | solver, elementy, analizy, siatkowanie | żadne |
| Aplikacja | interfejs, import CAD, warunki brzegowe, wyniki, raport | silnik |
| Serwer PHP | konta, projekty, katalog analiz, poziomy licencji | PDO (SQLite/MySQL) |
Silnik nie zna pojęcia „ściany" ani „krawędzi" — przyjmuje gotowe zbiory numerów węzłów. Wybór geometryczny należy do aplikacji, która ma własny model CAD i własny sposób zaznaczania. Dzięki temu biblioteka nie narzuca niczego co do formatu geometrii i pozostaje użyteczna poza tym konkretnym interfejsem.
Warstwa serwerowa jest świadomie wąska. Serwer przechowuje projekty i rozstrzyga o poziomach odpłatności analiz — nie liczy niczego. Praca bez konta jest w pełni możliwa: skoro analiza dzieje się w przeglądarce, serwer nie jest do niej potrzebny.
Pierwotne założenie — obliczenia w przeglądarce — ma jedną znaną wadę: zadanie zajmuje wątek karty, więc przy większych siatkach interfejs przestaje odpowiadać. System rozwiązuje to trójstopniowo, wybierając automatycznie najlepszą dostępną drogę.
| Droga | Moc | Interfejs | Warunek |
|---|---|---|---|
| Agent lokalny (Node) | wiele rdzeni, pełna pamięć | płynny | uruchomiony przez użytkownika |
| Wątek roboczy przeglądarki | jeden rdzeń | płynny | strona podana przez serwer |
| Wątek główny karty | jeden rdzeń | zablokowany | zawsze dostępna |
Agent liczący to ten sam silnik TOMes® uruchomiony poleceniem
node silnik/lokalnie/agent.js. Nasłuchuje wyłącznie na adresie pętli zwrotnej,
nie czyta plików z dysku i nie ma dostępu do bazy — jego jedynym wejściem jest definicja
zadania w treści zapytania. Zadania liczy w osobnych wątkach, dzięki czemu odpowiada na pytanie
o stan także w trakcie faktoryzacji, i obsługuje kilka zadań równolegle. Hosting nadal nie
liczy niczego: agent stoi na komputerze użytkownika.
To samo zadanie można policzyć wsadowo, bez przeglądarki:
Plik zadania jest tą samą strukturą, którą aplikacja zapisuje przyciskiem „Zapisz zadanie", więc wariant obliczony w nocy skryptem jest dokładnie tym zadaniem, które widać na ekranie.
Solver raportuje etapy pracy przez opcjonalny kanał postępu: budowę modelu, składanie macierzy, wybór solvera, kolejne iteracje i kontrolę równowagi. Przy liczeniu przez agenta zdarzenia wędrują do przeglądarki strumieniem NDJSON w trakcie obliczeń, a nie po nich. Najcenniejszy jest przy analizie nieliniowej — ciąg norm residuum Newtona jest najlepszym wskaźnikiem zdrowia zadania: spadek kwadratowy (10⁻², 10⁻⁴, 10⁻⁸) oznacza, że wszystko jest w porządku, spadek liniowy albo stanie w miejscu — że model wymaga uwagi.
System otwiera 17 formatów wymiany i siatkowych — od STEP i IGES, przez natywny FreeCAD, po formaty trójkątowe STL, 3MF, glTF i COLLADA. Pliku nie wysyła się nigdzie; czyta go przeglądarka. Formatów natywnych systemów CAD (SLDPRT, IPT, CATPart) system świadomie nie otwiera: ich specyfikacje nie są publiczne, a wszystkie dostępne czytniki są płatne, zamknięte i serwerowe. Zamiast udawać obsługę, program rozpoznaje plik po rozszerzeniu i pokazuje polecenie eksportu do formatu wymiany dla konkretnego programu.
Podział na warstwy ma cenę i warto ją nazwać, bo obie usterki opisane niżej wyszły z niego wprost. Żadna nie objawiła się wyjątkiem ani błędnym wynikiem — obie powodowały, że program wyglądał na działający i po cichu robił co innego, niż użytkownik zamierzał. To najkosztowniejsza klasa błędów w oprogramowaniu inżynierskim, bo nie zostawia po sobie śladu w wyniku.
Aplikacja jest zbiorem modułów ES ładowanych wprost przez przeglądarkę — bez
budowania i bez pakowania, co jest świadomym założeniem całego systemu. Znacznik
wersji stał jednak wyłącznie przy znaczniku <script src>, czyli
przy jednym pliku wejściowym. Wszystko, co ten plik importuje — a to praktycznie
cały system — przeglądarka pobierała spod adresu bez znacznika i przechowywała
według własnej heurystyki.
Po wgraniu poprawki dawało to stan gorszy od zwykłego nieodświeżenia: część modułów nowa, część sprzed zmiany. Program uruchamiał się bez błędu i zachowywał jak złożenie dwóch wersji. Objawem, który to wykrył, był nowy moduł interfejsu pytający stary katalog warunków brzegowych o warunki właściwe dla akustyki — i pokazujący w odpowiedzi zestaw mechaniczny, bo starsza wersja katalogu rodziny akustycznej jeszcze nie znała.
Rozwiązaniem jest mapa importów budowana po stronie serwera: każdemu modułowi przypisywany jest adres ze znacznikiem czasu jego pliku. Zmiana jednej linii unieważnia dokładnie ten jeden moduł i tylko jego. Jest to zarazem jedyne miejsce w systemie, w którym serwer w ogóle wie o istnieniu plików silnika.
Zestaw warunków brzegowych zależy od analizy, bo zależy od niej fizyka: przy czystym przewodnictwie utwierdzenie nie wpływa na wynik w żaden sposób, a przy statyce — temperatura powierzchni. Druga lista mówi, czego dana analiza wymaga, żeby dało się ją uruchomić. Obie powstały niezależnie i rozjechały się bez śladu.
| Analiza | Warunki dostępne | Warunki wymagane | Skutek |
|---|---|---|---|
| Filtracja w ośrodku porowatym | ciśnienie na wlocie i wylocie | dwa warunki temperatury | „Policz" wygaszone na zawsze |
| Przepływ potencjalny | potencjał prędkości | podparcie i obciążenie | „Policz" wygaszone na zawsze |
Obie analizy stały w katalogu jako działające, obie miały opisaną dokładność potwierdzoną testem solvera i obie były nieosiągalne z interfejsu. Test solvera nie mógł tego wykryć, bo solver liczył je bez zarzutu — zadanie po prostu nigdy do niego nie docierało.
Poprawka polega na przeniesieniu obu list do jednego pliku i sprawdzaniu ich zgodności testem: każdy wymagany typ warunku musi należeć do zestawu dopuszczonego przy tej samej analizie. Jest to warunek konieczny, nie dostateczny — ale wystarcza, żeby analiza nie stała się nieosiągalna wskutek niezależnej zmiany w drugiej liście. Wymagania nie dają się przy tym uśrednić: drgania własne nie potrzebują obciążenia, akustyka wnętrza zamkniętego nie potrzebuje żadnego warunku, bo ściana sztywna jest warunkiem naturalnym równania Helmholtza, a oba zadania polowe potrzebują dwóch warunków brzegowych, bo jeden daje pole stałe i zerowy przepływ.
Pytanie „skąd wiadomo, że to liczy poprawnie" pada przy każdym odbiorze i jest pytaniem zasadnym. Odpowiedzią jest zestaw weryfikacyjny uruchamiany jednym poleceniem:
Zestaw nie sprawdza, czy program się nie wywraca — sprawdza, czy wynik zgadza się z rozwiązaniem znanym niezależnie od tego programu. Każdy test podaje wartość otrzymaną, wartość odniesienia, różnicę i tolerancję, więc jego wynik da się ocenić bez zaglądania do kodu.
| Grupa | Testy | Przedmiot weryfikacji |
|---|---|---|
| Rdzeń algebraiczny | 17 | format CSR, LDLᵀ wobec eliminacji Gaussa, bezwładność Sylvestera |
| Solvery | 17 | PCG z IC(0), obroty Jacobiego, iteracja podprzestrzeni |
| Analizy liniowe | 28 | rozciąganie ścisłe, test łaty, belka Timoszenki, zbieżność siatki, wyboczenie Eulera |
| Analizy nieliniowe | 16 | niezmienniczość obrotu, plastyczność jednoosiowa, nośność graniczna, długość łuku |
| Kontakt | 8 | nieprzenikanie, jednostronność, nasycenie tarcia Coulomba |
| Analizy cieplne | 8 | Fourier, konwekcja, stygnięcie, naprężenia termiczne |
| Dynamika | 11 | rezonans, Newmark, dynamika jawna, spektrum, drgania losowe |
| Pola skalarne | 11 | naprężenie wstępne, rezonanse akustyczne, prawo Darcy’ego |
| Zmęczenie | 19 | Basquin, poprawki średnie, Coffin-Manson, Paris |
| Reologia | 9 | pełzanie Nortona, człon Arrheniusa, relaksacja lepkosprężysta |
| Wtrysk tworzywa | 10 | model Cross, Hele-Shaw, czas wypełnienia, linie łączenia |
| Hipersprężystość, dynamika nieliniowa, śruby, podmodelowanie | 24 | neo-Hooke wobec postaci zamkniętej, Newmark z Newtonem, siła w trzpieniu, ruch sztywny na cięciu |
| Wskazywanie warunków na modelu | 13 | tryby wyboru, tolerancje kierunkowe, trwałość wskazania po przegenerowaniu siatki |
| Podstruktury i CMS | 13 | kondensacja ścisła przy obciążeniu na brzegu, zbieżność jednostronna, symetria macierzy superelementu |
| Spójność interfejsu z silnikiem | 16 | wymagania analizy wobec dostępnych warunków, kompletność katalogu |
| Publiczne API | 19 | zgodność interfejsu, jednostki, sesja wieloanalizowa, narzędzia siatki |
| Sprawdzenie | Odniesienie | Wynik MES | Różnica |
|---|---|---|---|
| Rozciąganie jednoosiowe (TET4/TET10/HEX8) | 1,1905·10⁻⁴ m | 1,1905·10⁻⁴ m | 0,0000% |
| Test łaty — jednorodność σxx | 0 | 10⁻¹⁵ | dokładnie |
| Wspornik zginany, TET10 (Timoszenko) | 9,5981·10⁻⁴ m | 9,4314·10⁻⁴ m | 1,74% |
| Wspornik zginany, HEX8 | 9,5981·10⁻⁴ m | 9,7134·10⁻⁴ m | 1,20% |
| Wspornik zginany, TET4 | 9,5981·10⁻⁴ m | 6,4323·10⁻⁴ m | 32,98% |
| Pierwsza częstość, TET10 (Euler-Bernoulli) | 92,84 Hz | 93,34 Hz | 0,54% |
| Pierwsza częstość, TET4 | 92,84 Hz | 216,05 Hz | 132,72% |
| Siła krytyczna wyboczenia (Euler) | 1165,85 N | 1146,78 N | 1,64% |
| Obrót ciała sztywnego o 30° — naprężenie | 0 | 2,08·10⁻¹⁶ | dokładnie |
| Statyka nieliniowa wobec liniowej (małe obciążenie) | −1,8789·10⁻⁵ m | −1,8788·10⁻⁵ m | 0,0005% |
| Równowaga globalna (suma sił i reakcji) | 0 | 10⁻¹¹ | dokładnie |
Wynik metody elementów skończonych zależy od gęstości siatki i nie ma sposobu, by ocenić go na podstawie jednego przebiegu. Zestaw sprawdza to wprost, a system udostępnia badanie zbieżności użytkownikowi — liczy zadanie na trzech zagęszczeniach i pokazuje zmianę wyniku.
| Element | Stopnie swobody | Błąd wobec Timoszenki |
|---|---|---|
| TET4 | 297 | 64,9% |
| TET4 | 1575 | 33,0% |
| TET10 | 297 | 4,5% |
Element kwadratowy na siatce pięciokrotnie mniejszej daje wynik siedmiokrotnie dokładniejszy. To jest istota różnicy między rzędem interpolacji a zagęszczeniem — i powód, dla którego zagęszczanie siatki elementami liniowymi jest kosztowną drogą donikąd.
Cztery zadania policzone dwiema drogami: wzorem zamkniętym i metodą elementów skończonych. Wszystkie liczby pochodzą z rzeczywistych przebiegów opisywanego systemu — nie są przepisane z literatury. Materiałem jest stal konstrukcyjna S355: E = 210 GPa, ν = 0,3, ρ = 7850 kg/m³, α = 12·10⁻⁶ 1/K.
Wspornik o przekroju prostokątnym 40 × 30 mm i długości 500 mm, utwierdzony jednym końcem i obciążony na swobodnym końcu siłą P = 2000 N.
| Siatka TET10 | DOF | Ugięcie | Różnica | σzast maks. | Równowaga | Czas |
|---|---|---|---|---|---|---|
| 12 × 3 × 3 | 3 675 | 4,3594 mm | −1,41% | 128,47 MPa | 1,0·10⁻¹⁰ | 0,46 s |
| 20 × 4 × 4 | 9 963 | 4,3806 mm | −0,93% | 138,15 MPa | 1,9·10⁻¹⁰ | 1,71 s |
| Postać | MES | Wzór | Uwaga |
|---|---|---|---|
| 1 | 100,67 Hz | 100,26 Hz | zginanie w płaszczyźnie słabszej, +0,41% |
| 2 | 133,82 Hz | — | zginanie w płaszczyźnie sztywniejszej |
| 3 | 620,10 Hz | 628,33 Hz | druga postać zginania, −1,31% |
| 4 | 814,06 Hz | — | druga postać w płaszczyźnie sztywniejszej |
| 5 | 1390,84 Hz | — | skrętna |
| 6 | 1693,29 Hz | 1759,35 Hz | trzecia postać zginania, −3,75% |
Pręt o przekroju kwadratowym 20 × 20 mm i długości 500 mm, utwierdzony jednym końcem (długość wyboczeniowa 2L), obciążony siłą próbną 1000 N.
| Wielkość | Wzór Eulera | MES | Różnica |
|---|---|---|---|
| Mnożnik λ₁ | 27,635 | 27,753 | +0,43% |
| Siła krytyczna | 27 634,9 N | 27 752,9 N | +0,43% |
| Mnożnik λ₂ | 27,635 | 27,764 | postać w drugiej płaszczyźnie |
Przekrój kwadratowy ma jednakową sztywność w obu płaszczyznach, więc dwie pierwsze postacie wyboczenia są niemal identyczne (27,753 i 27,764). Jest to dokładnie ten przypadek, dla którego iteracja podprzestrzeni wymaga bloku wektorów większego od liczby szukanych postaci — bez zapasu jedna z tych dwóch bliskich postaci zostałaby zgubiona.
Pręt 300 × 20 × 20 mm, obustronnie skrępowany osiowo, ogrzany o ΔT = 60 K.
| Wielkość | Wzór | MES | Różnica |
|---|---|---|---|
| σxx w środku pręta | −151,20 MPa | −151,20 MPa | 0,00% |
| Naprężenie zastępcze | 151,20 MPa | 151,20 MPa | 0,00% |
Zgodność jest dokładna, i tak być powinno: stan naprężenia jest tu jednorodny, a element o dowolnym rzędzie interpolacji odtwarza pole stałe bez błędu. Jest to ten sam mechanizm, który sprawdza test łaty — i dlatego wynik niezgodny w tym zadaniu oznaczałby błąd w macierzy konstytutywnej, a nie niedokładność dyskretyzacji.
Wykaz poniższy jest częścią systemu, nie dodatkiem do dokumentacji: te same pozycje zwraca
programowo funkcja capabilities().knownLimitations i pojawiają się w raporcie
z każdej analizy.
Katalog systemu obejmuje 54 pozycje, z których 38 jest policzalnych i potwierdzonych testem. Pozostałe 16 jest widocznych jako zapowiedzi — nie da się ich wybrać ani wycenić w panelu administratora. Powody nie są jednorodne i warto je rozdzielić:
| Grupa | Pozycji | Przyczyna |
|---|---|---|
| Rodzina CFD (RANS, LES, ściśliwy, wielofazowy, maszyny wirnikowe, FSI) | 9 | wymaga siatki przyściennej, stabilizacji i całkowania o zupełnie innych wymaganiach niż mechanika ciała stałego |
| Przetwórstwo tworzyw: docisk i skurcz, paczenie, układ chłodzenia | 3 | nadbudowa nad policzalnym już wypełnianiem formy; wymaga sprzężenia przepływu z polem temperatury |
| Złącza i przeguby, połączenia śrubowe | 2 | więzy wielopunktowe (MPC), których model jeszcze nie zna |
| Dynamika wirników | 1 | macierze żyroskopowe i zagadnienie własne w dziedzinie zespolonej |
| Wibroakustyka | 1 | sprzężenie dwustronne konstrukcja–powietrze |
| Symbol | Znaczenie | Jednostka SI |
|---|---|---|
| σ, ε | tensor naprężenia i odkształcenia | Pa, — |
| D | macierz konstytutywna materiału | Pa |
| B | macierz odkształcenie–przemieszczenie | 1/m |
| N | funkcje kształtu elementu | — |
| K, K_g, K_t | macierz sztywności: sprężysta, geometryczna, styczna | N/m |
| M, C | macierz mas i tłumienia | kg, N·s/m |
| u, F | wektor przemieszczeń węzłowych i sił | m, N |
| J | jakobian przekształcenia geometrycznego | m |
| E, ν, ρ | moduł Younga, liczba Poissona, gęstość | Pa, —, kg/m³ |
| λ, μ | stałe Lamégo | Pa |
| κ | moduł ściśliwości objętościowej | Pa |
| σ_y, H | granica plastyczności, moduł wzmocnienia | Pa |
| Δγ | mnożnik plastyczny | — |
| ω, f | częstość kołowa i częstotliwość | rad/s, Hz |
| φ, λ | postać własna, mnożnik obciążenia krytycznego | — |
| α | współczynnik rozszerzalności cieplnej | 1/K |
| ΔK, C, m | zakres współczynnika intensywności naprężeń i stałe Parisa | Pa·√m |