MKE analiza

Teorija I reda

K - linearna matrica krutosti, u - vektor generalisanih pomeranja, F - vektor opterećenja. Pretpostavke: male deformacije i mala pomeranja - ravnoteža se postavlja na nedeformisanom sistemu, pa važi superpozicija.

Izvori: [1], [6].

Linearizovana teorija II reda (P-Δ i P-δ)

K_G - geometrijska matrica krutosti (matrica početnih napona), N_I - normalne sile po teoriji I reda istog slučaja. Šta je i kada: uticaj normalne sile na savijanje (P-Δ - obrtanje tetive, P-δ - krivljenje štapa) kod vitkih pritisnutih elemenata. Pretpostavke: ravnoteža na deformisanom sistemu, linearizacija (N se ne menja sa pomeranjima); rezultat je merodavan samo kada je λ_cr > 1 (K + K_G pozitivno definitna) - inače program upozorava.

Izvori: [1] §8.1, [5].

Bifurkaciona stabilnost i bočno-torziono izvijanje

λ - kritični faktor opterećenja (u EN 1993-1-1 α_cr; kritično opterećenje = λ_cr·zadato opterećenje), φ - oblik izvijanja (proizvoljno normiran), N, M₂, M₃ - presečne sile po teoriji I reda. Kod bočno-torzionog izvijanja K_G sadrži i članove od momenata savijanja (sprega savijanja i uvijanja), a kritični moment je M_cr = λ_cr·M; računa se sa vitoperenjem ili bez njega (prekidač u tabu Model), uz napadnu tačku opterećenja zadatu ekscentricitetom e₂, e₃ u odnosu na težište. Pretpostavke: linearna teorija stabilnosti savršenog sistema (bez imperfekcija), N i M proporcionalni opterećenju.

U programu: sopstveni problem rešava se gustim rešavačem ili iteracijom u potprostoru (subspace iteration) uz proveru Sturmovim nizom; ako iteracija ne uspe, koristi se gusti rešavač. Program sam bira rešavač, a korisnika pita kad procenjeno vreme premašuje 8 s.

Izvori: [2] §8.2.5, §11.6, [5], [6], [7], [8].

GNA, korak 1 - teorija

f_int - unutrašnje sile, F₀ - opterećenje prethodnih faza, q - referentno opterećenje faze, λ - faktor opterećenja, K_E - elastična matrica krutosti, K_T = ∂f_int/∂u - tangentna matrica krutosti, K_N - član od promene normalne sile sa pomeranjima. Pretpostavke pojedinih teorija date su ispod.

Nazivi analiza: opcija ima međunarodni naziv (first-order / second-order analysis - EN 1993-1-1 §5.2; co-rotational analysis - Crisfield, Vol. 1, §7), a u zagradi ustaljeni domaći naziv (teorija I, II, III reda). „Teorija III reda” u ovom programu znači analizu sa velikim pomeranjima i rotacijama i malim deformacijama (korotaciona formulacija), a ne teoriju velikih (konačnih) deformacija.

Pojedinačne teorije (materijal svuda linearno elastičan): Linearised second-order analysis (linearizovana teorija II reda) - male deformacije i rotacije (sin φ ≈ φ, cos φ ≈ 1), umerena pomeranja; ravnoteža na deformisanoj konfiguraciji, ali sa normalnim silama teorije I reda (S^II·D^II ≈ S^I·D^II; S - normalna sila, D - pomeranje, eksponent - red teorije po kojoj se veličina računa); promena N usled deformacije i skraćenje tetive se zanemaruju; [K_E + K_G(N^I)]·u = F je jedno linearno rešenje, a pošto N^I zavisi od opterećenja, superpozicija ne važi; det(K_E + λ·K_G) = 0 daje kritični faktor (bifurkacija). Second-order analysis - P-Δ only (teorija II reda - samo P-Δ): geometrijska matrica sadrži samo član N/L (N - normalna sila, L - dužina KE) uz relativno poprečno pomeranje krajeva KE; izostavljaju se lokalni efekat P-δ (doprinos normalne sile preko zakrivljenosti unutar KE) i skraćenje tetive, a elastično savijanje KE ostaje uključeno; pri podeli štapa na više KE deo efekta P-δ obuhvataju pomeranja međučvorova, pa razlika prema punoj teoriji II reda zavisi od mreže. Second-order analysis - P-Δ and P-δ (teorija II reda - P-Δ i P-δ): konzistentna geometrijska matrica k_G; zanemaruju se skraćenje tetive usled savijanja (član ½·ψ² u izduženju, ψ - rotacija tetive) i viši članovi rotacija. Co-rotational analysis (teorija III reda): f_int(u) = F₀ + λ·q na tekućoj deformisanoj konfiguraciji; korotaciona formulacija odvaja kruto kretanje svakog KE od malih lokalnih deformacija, za koje važi linearna teorija štapa ([3] §7.3) - nije teorija velikih (konačnih) deformacija, a tačnost ostaje uslovljena pretpostavkama elementa i mrežom. Savršen sistem (prav stub, nosač bez početne imperfekcije) u tački bifurkacije sam ne prelazi na izvijeni oblik: ostaje na primarnoj grani, koja je posle nje nestabilna - program javlja da je tačka bifurkacije prođena; za praćenje izvijenog oblika treba uvesti početnu imperfekciju.

Izvori: [1] §8.1, [2] §6.1, [3] pogl. 7, [9].

GNA, korak 2 - pristup rešavanju

n - redni broj inkrementa, Δu_n, Δλ_n - priraštaji, r - neuravnotežene (rezidualne) sile. Jednokoračni pristup postoji samo uz teoriju I reda i linearizovanu teoriju II reda. Bez iteracija rešenje se udaljava od ravnotežnog puta (greška je približno srazmerna 1/N, N - broj koraka); program na kraju ispisuje ‖r‖/‖F‖. Inkrementalno-iterativni pristup u svakom inkrementu iterira (prediktor + korektor) dok ravnoteža ne bude zadovoljena.

Kod čisto inkrementalnog pristupa u teoriji III reda greška je kod vitkih štapova velika i pri hiljadama koraka: tangentni korak izdužuje tetivu KE, pa nastaje fiktivna normalna sila ½·ψ²·E·A (ψ - rotacija tetive u koraku, E·A - aksijalna krutost). Pristup sa korekcijom neuravnoteženih sila odgovara jednoj Newton-Raphson iteraciji po koraku: odstupanje od ravnotežnog puta je mnogo manje, ali ravnoteža i dalje nije zadovoljena. Kod inkrementalno-iterativnog pristupa svaka prihvaćena tačka zadovoljava ravnotežu u zadatoj toleranciji; to samo po sebi ne dokazuje da je praćena željena grana niti da je ravnoteža stabilna.

Izvori: [3] §1.2, [4] §2.4, [1] (12.3)-(12.6).

GNA, korak 3 - kontrola puta

u_k - kontrolisana translacija U, V ili W (ona sa najvećim pomeranjem po teoriji I reda), Δl - dužina luka, ψ - težinski faktor opterećenja (0 - cilindrični, 1 - sferni uslov), u_I - pomeranje od opterećenja q na početku faze (skaliranje opterećenja): K_E^-1·q u teoriji II reda, K_T(u₀)^-1·q u teoriji III reda (u₀ - stanje na početku faze; u prvoj fazi K_T(0) = K_E), Δu₁ - prvi prediktor inkrementa (Riks: fiksna normalna ravan), Δu - tekući priraštaj (Ramm: ažurirana normalna ravan), δu - iterativna popravka. Norma luka se računa po translacijama. Dužina luka nije ponuđena uz pristupe bez iteracija.

Kontrola pomeranja ([10]; [3] §9.3.2.3): Δλ = Δu_c/δu_t,c, u iteracijama δλ = -δū_c/δu_t,c (u_c - kontrolisana translacija, δu_t = K_T^-1·q - tangentno pomeranje od referentnog opterećenja, δū = K_T^-1·r - popravka od neuravnoteženih sila, indeks c - komponenta u_c); kontrolisani stepen slobode navodi se u poruci šeme, a uz pristup bez iteracija izvršava se samo inkrementalni prediktor. Prolazi graničnu tačku po opterećenju ako je izabrano pomeranje pogodan parametar puta. Riks-Wempner ([3] §9.3.2.2, sl. 9.11a): δλ = -Δu₁^T·δū/(Δu₁^T·δu_t) - bez izbora korena, ali pri velikom koraku može da promaši put. Ramm ([3] §9.3.2.2, sl. 9.11b): δλ = -Δu^T·δū/(Δu^T·δu_t). Crisfield ([3] §9.3.2.1 i 9.4): δλ iz kvadratne jednačine, bira se koren sa najmanjim uglom prema dosadašnjem priraštaju; kod sfernog uslova (ψ = 1), ako jednačina nema realan koren, korak se prepolovljava. Dužina luka ne rešava izbor grane u tački bifurkacije. Kod kontrole opterećenja tangentna matrica K_T je u graničnoj tački singularna, a iza nje nema susedne ravnoteže za veće λ.

Izvori: [3] pogl. 9, [4] pogl. 4, [10], [11], [12], [13].

GNA, korak 4 - iterativni postupak

Newton-Raphson - K_T se obnavlja u svakoj iteraciji (lokalno kvadratna konvergencija - blizu rešenja, uz konzistentnu i nesingularnu tangentu; Bathe, 1996, §8.4.1); modified Newton-Raphson - K_T sa početka inkrementa (linearna); initial stiffness (metoda početne krutosti; u literaturi i „initial stress method”) - stalna K₀: u teoriji II reda K_E, u teoriji III reda K_T(u₀) sa početka faze (linearna; blizu granične tačke ne konvergira); quasi-Newton BFGS - sekantna dopuna inverzne matrice (lokalno superlinearna). Kriterijum konvergencije u teoriji II reda je po sili, pomeranju, energiji ili kombinovani, a u teoriji III reda po sili.

Newton-Raphson ima najskuplju iteraciju (nova tangenta i faktorizacija u svakoj iteraciji), a modified Newton-Raphson jeftinu (bez nove faktorizacije). Initial stiffness po pravilu najsporije konvergira i može da divergira; naziv „initial stress method” ne odnosi se na geometrijsku matricu K_G („matrica početnih napona”). Quasi-Newton BFGS ([14]) polazi od tangente sa početka inkrementa i ažurira njenu inverznu matricu iz promene unutrašnjih sila (sekantni uslov; [2] §8.4.2).

Izvori: [2] §8.4.1-8.4.4, [3] §1.2, §9.3-9.4, [14].

Faze i imperfekcije

Dejstva deluju u fazama, istovremeno ili jedno za drugim, sa postupkom po fazi. Oblik imperfekcije jednak je izabranom obliku izvijanja faze, a amplituda je e₀ (najveće translaciono pomeranje oblika, a kod pretežno torzionog oblika najveća komponenta); pomoćni proračun e₀ = L/α po EN 1993-1-1, Tab. 5.1 (L - dužina štapa, α = 350, 300, 250, 200, 150 za krive izvijanja a₀, a, b, c, d).

Imperfekcija ne menja geometriju, već ulazi kao ekvivalentno opterećenje (N - normalne sile na kraju faze, sa opterećenjem prethodnih faza, u₀ - imperfekcija): tačno u linearizovanoj teoriji II reda, približno u teoriji III reda. Imperfekcija iz ranije faze ostaje i u kasnijim fazama, a njeno ekvivalentno opterećenje raste sa normalnom silom svake naredne faze.

Izvori: [7], [15].

First-order theory

K - linear stiffness matrix, u - vector of generalised displacements, F - load vector. Assumptions: small strains and small displacements - equilibrium is written on the undeformed structure, so superposition holds.

Sources: [1], [6].

Linearised second-order theory (P-Δ and P-δ)

K_G - geometric (initial-stress) stiffness matrix, N_I - first-order axial forces of the same case. What and when: the effect of the axial force on bending (P-Δ - chord rotation, P-δ - member curvature) in slender compressed members. Assumptions: equilibrium on the deformed structure, linearised (N does not change with the displacements); the result is valid only when λ_cr > 1 (K + K_G positive definite) - otherwise the program warns.

Sources: [1] §8.1, [5].

Bifurcation stability and lateral-torsional buckling

λ - critical load factor (α_cr in EN 1993-1-1; critical load = λ_cr·applied load), φ - buckling mode (arbitrarily normalised), N, M₂, M₃ - first-order internal forces. For lateral-torsional buckling K_G also contains terms from the bending moments (coupling of bending and twisting), and the critical moment is M_cr = λ_cr·M; it is computed with or without warping (switch in the Model tab), with the load point given by its eccentricity e₂, e₃ from the centroid. Assumptions: linear stability theory of a perfect system (no imperfections), N and M proportional to the load.

In the program: the eigenproblem is solved by the dense solver or by subspace iteration with a Sturm sequence check; if the check fails, the dense solver. The program chooses by itself and asks the user when the estimated time exceeds 8 s.

Sources: [2] §8.2.5, §11.6, [5], [6], [7], [8].

GNA, step 1 - theory

f_int - internal forces, F₀ - load of the earlier phases, q - reference load of the phase, λ - load factor, K_E - elastic stiffness matrix, K_T = ∂f_int/∂u - tangent stiffness matrix, K_N - term from the change of the axial force with the displacements. Second-order analysis (P-Δ only, or P-Δ and P-δ): moderate displacements and small rotations; N depends on the displacements, so the problem is nonlinear; it does not describe the behaviour beyond a limit point, and outside its range of validity (rotation above 0.1 rad) the program warns. Co-rotational analysis (third-order theory in the local terminology; co-rotational formulation, planar or spatial, also with warping): the rigid-body motion of each element is separated by a co-rotating frame in which the strains are small, so the displacements and rotations of the structure may be large - provided that the strains and relative rotations within each element are small (a sufficiently fine mesh).

Analysis names: each option carries the international name (first-order / second-order analysis - EN 1993-1-1 §5.2; co-rotational analysis - Crisfield, Vol. 1, §7); the Serbian interface adds the established local name in parentheses (first-, second-, third-order theory). “Third-order theory” in this program means an analysis with large displacements and rotations and small strains (co-rotational formulation), not a large (finite) strain theory.

Individual theories (material linear elastic throughout): Linearised second-order analysis - small strains and rotations (sin φ ≈ φ, cos φ ≈ 1), moderate displacements; equilibrium on the deformed configuration, but with first-order axial forces (S^II·D^II ≈ S^I·D^II; S - axial force, D - displacement, superscript - order of the theory by which the quantity is computed); the change of N due to deformation and the chord shortening are neglected; [K_E + K_G(N^I)]·u = F is one linear solution, and since N^I depends on the load, superposition does not hold; det(K_E + λ·K_G) = 0 gives the critical factor (bifurcation). Second-order analysis - P-Δ only: the geometric matrix contains only the term N/L (N - axial force, L - element length) with the relative transverse displacement of the element ends; the local P-δ effect (axial-force contribution through the curvature within the element) and the chord shortening are omitted, while elastic bending of the element remains included; when a member is divided into several elements, part of the P-δ effect is captured by the intermediate node displacements, so the difference from full second-order theory depends on the mesh. Second-order analysis - P-Δ and P-δ: consistent geometric matrix k_G; the chord shortening due to bending (term ½·ψ² in the elongation, ψ - chord rotation) and higher-order rotation terms are neglected. Co-rotational analysis (third-order theory): f_int(u) = F₀ + λ·q on the current deformed configuration; the co-rotational formulation separates the rigid-body motion of each element from small local strains, for which linear beam theory holds ([3] §7.3) - it is not a large (finite) strain theory, and accuracy remains subject to the element assumptions and the mesh. A perfect system (straight column, member without initial imperfection) does not move to the buckled shape by itself at a bifurcation point: it stays on the primary branch, which is unstable beyond it - the program reports that the bifurcation point was passed; to trace the buckled shape an initial imperfection must be introduced.

Sources: [1] §8.1, [2] §6.1, [3] ch. 7, [9].

GNA, step 2 - solution approach

n - increment number, Δu_n, Δλ_n - increments, r - unbalanced (residual) forces. The single-step approach exists only with first-order and linearised second-order theory. Without iterations the solution drifts from the equilibrium path (the error is approximately proportional to 1/N, N - number of steps); the program prints ‖r‖/‖F‖ at the end. The incremental-iterative approach iterates within each increment (predictor + corrector) until equilibrium is satisfied.

With the purely incremental approach in third-order theory the error for slender members is large even with thousands of steps: the tangent step stretches the element chord, giving a spurious axial force ½·ψ²·E·A (ψ - chord rotation in the step, E·A - axial stiffness). The approach with correction of the unbalanced forces corresponds to one Newton-Raphson iteration per step: the the departure from the equilibrium path is much smaller, but equilibrium is still not satisfied. With the incremental-iterative approach every accepted point satisfies equilibrium within the given tolerance; this alone does not prove that the intended branch is followed or that the equilibrium is stable.

Sources: [3] §1.2, [4] §2.4, [1] (12.3)-(12.6).

GNA, step 3 - path control

u_k - controlled translation U, V or W (the one with the largest first-order displacement), Δl - arc length, ψ - load weighting (0 cylindrical, 1 spherical constraint), u_I - displacement due to the load q at the start of the phase (load scaling): K_E^-1·q in second-order theory, K_T(u₀)^-1·q in third-order theory (u₀ - state at the start of the phase; in the first phase K_T(0) = K_E), Δu₁ - first predictor of the increment (Riks: fixed normal plane), Δu - current increment (Ramm: updated normal plane), δu - iterative correction. The arc norm is computed over translations. The arc-length method is not offered with the approaches without iterations.

Displacement control ([10]; [3] §9.3.2.3): Δλ = Δu_c/δu_t,c, in the iterations δλ = -δū_c/δu_t,c (u_c - controlled translation, δu_t = K_T^-1·q - tangent displacement due to the reference load, δū = K_T^-1·r - correction due to the unbalanced forces, subscript c - component u_c); the controlled degree of freedom is given in the scheme message, and with an approach without iterations only the incremental predictor is performed. It passes a load limit point if the chosen displacement is a suitable path parameter. Riks-Wempner ([3] §9.3.2.2, Fig. 9.11a): δλ = -Δu₁^T·δū/(Δu₁^T·δu_t) - no root selection, but with a large step it can miss the path. Ramm ([3] §9.3.2.2, Fig. 9.11b): δλ = -Δu^T·δū/(Δu^T·δu_t). Crisfield ([3] §9.3.2.1 and 9.4): δλ from a quadratic equation, the root with the smallest angle to the previous increment is chosen; for the spherical constraint (ψ = 1) the step is halved when there is no real root. The arc-length method does not resolve branch selection at a bifurcation point. Under load control the tangent matrix K_T is singular at a limit point, and beyond it there is no adjacent equilibrium for a larger λ.

Sources: [3] ch. 9, [4] ch. 4, [10], [11], [12], [13].

GNA, step 4 - iterative procedure

Newton-Raphson - K_T is updated in every iteration (locally quadratic convergence - near the solution, with a consistent and nonsingular tangent; Bathe, 1996, §8.4.1); modified Newton-Raphson - K_T from the start of the increment (linear); initial stiffness (also called the initial stress method) - a constant K₀: K_E in second-order theory, K_T(u₀) from the start of the phase in third-order theory (linear; does not converge near a limit point); quasi-Newton BFGS - secant update of the inverse matrix (locally superlinear). The convergence criterion in second-order theory is force, displacement, energy or combined, and in third-order theory force.

Newton-Raphson has the most expensive iteration (a new tangent and factorisation in every iteration); quadratic convergence holds locally - when the iterate is sufficiently close to the solution and the tangent is consistent, nonsingular and does not change abruptly. Modified Newton-Raphson has a cheap iteration (no new factorisation). Initial stiffness has, as a rule, the slowest convergence and may diverge; K₀ in the first phase of third-order theory equals K_E. The name “initial stress method” does not refer to the geometric matrix K_G (the “initial-stress matrix”). Quasi-Newton BFGS ([14]) starts from the tangent at the beginning of the increment and updates its inverse from the change of internal forces (secant condition) ([2] §8.4.2).

Sources: [2] §8.4.1-8.4.4, [3] §1.2, §9.3-9.4, [14].

Phases and imperfections

Loads act in phases, simultaneously or one after another, with the procedure chosen per phase. An imperfection has the shape of a chosen buckling mode of the phase and the amplitude e₀ (largest translational displacement of the shape; for a predominantly torsional shape, its largest component); helper e₀ = L/α to EN 1993-1-1, Table 5.1 (L - member length, α = 350, 300, 250, 200, 150 for buckling curves a₀, a, b, c, d).

An imperfection does not change the geometry; it enters as an equivalent load (N - axial forces at the end of the phase, including earlier phases, u₀ - imperfection): exact in linearised second-order theory, approximate in third-order theory. An imperfection from an earlier phase remains in the later phases, and its equivalent load grows with the axial force of every following phase.

Sources: [7], [15].

Jednačina kretanja i pobuda

M, C, K - matrice mase, prigušenja i krutosti; u, u̇, ü - vektori pomeranja, brzina i ubrzanja (uz akcelerogram relativni u odnosu na tlo); F(t) - vektor opterećenja; ι - vektor uticajnih koeficijenata (uz model iz taba Model ι_j = 1 za čekirane stepene slobode, a 0 za ostale; uz ručno zadate matrice brojne vrednosti ι_j, podrazumevano 1); a_g - ubrzanje tla. Uz dejstva iz Modela statički i dinamički deo se sabiraju (superpozicija). Pretpostavke: linearan sistem (M, C, K ne zavise od vremena ni od odgovora), ista pobuda svih oslonaca, dinamičko dejstvo samo na stepenima slobode sa masom.

Izvori: [1] pogl. 9, §13.1, [2].

Direktna numerička integracija (DNI): Newmark i metoda centralnih razlika

Šta je i kada: jednačine M·ü + C·u̇ + K·u = F rešavaju se zajedno u svakom koraku Δt; važi i za neklasično prigušenje.

i - redni broj koraka, γ, β - parametri postupka. Newmark sa konstantnim srednjim ubrzanjem je bezuslovno stabilan i nema numeričkog prigušenja (amplituda ne opada), ali produžava period oscilovanja: T̄/T = (ω·Δt/2)/arctg(ω·Δt/2) ≈ 1 + (ω·Δt)²/12, gde je T̄ numerički dobijen period, T tačan sopstveni period, ω = 2π/T kružna frekvencija tona i Δt korak integracije (za Δt = T/10 produženje je oko 3,2 %, za Δt = T/20 oko 0,8 %; [1] §5.5, sl. 5.5.2, [3] §9.4.3, sl. 9.6). Metoda centralnih razlika (MCR) skraćuje period, T̄/T ≈ 1 - (ω·Δt)²/24, i uslovno je stabilna: Δt < T_n,min/π (T_n,min - najkraći sopstveni period); program to proverava i, ako uslov nije ispunjen, proračun se ne pokreće. Korak Δt zadaje korisnik, a program ga samo usklađuje sa korakom zapisa dejstva (v. ispod).

Izvori: [1] §5.3, §5.4, [3] §9.2.1, §9.2.3, [4].

Korak numeričke integracije i korak zapisa dejstva

Korak numeričke integracije Δt bira korisnik. Zapis dejstva (sila F(t) ili ubrzanje tla a_g(t)) ima svoj korak Δt_zapis: iz polja „Korak zapisa dejstava” za jednokolonske fajlove (prazno polje: Δt_zapis = Δt), a za dvokolonski fajl (t, vrednost) iz kolone vremena, koja mora rasti sa stalnim korakom (inače se fajl ne učitava). Integracija se uvek sprovodi sa stalnim korakom h, a vrednosti zapisa u trenucima proračuna t_j = j·h određuju se linearnom interpolacijom:

m - ceo broj (Δt_zapis = m·h; ⌈x⌉ - najmanji ceo broj ≥ x), F_i - i-ta vrednost zapisa, t_i - njen trenutak. Za Δt ≤ Δt_zapis trenutak svake vrednosti zapisa je i trenutak proračuna, pa se karakteristične vrednosti dejstva (npr. ekstremne) ne preskaču, a između njih se dejstvo menja linearno - to je ista pretpostavka koju koriste i postupci integracije (pobuda se unutar koraka menja linearno, [1] §5.2). Kada se h razlikuje od zadatog Δt, uz rezultate se ispisuje napomena sa usvojenim korakom. Za Δt > Δt_zapis vrednosti zapisa između trenutaka proračuna se preskaču; uz rezultate se ispisuje upozorenje sa najvećom vrednošću zapisa i najvećom vrednošću u trenucima proračuna. Ako zapisi imaju različite korake, h se usklađuje sa najmanjim; tačke zapisa čiji korak nije celobrojni umnožak tog koraka interpoliraju se (napomena uz rezultate). Posle kraja zapisa dejstvo je 0. Kada je Δt_zapis = Δt, zapis se prepisuje doslovno, pa su rezultati isti kao bez usklađivanja.

Preporuka: Δt ≤ T/10 (T - najkraći period tonova koji bitno doprinose odgovoru) i Δt ≤ Δt_zapis; za metodu centralnih razlika mora biti i Δt < T_min/π (T_min - najkraći sopstveni period sistema).

Izvori: [1] §5.2, §5.5.

Modalna superpozicija

q_n - modalna koordinata, φ̃_n - oblik normiran masom, ω_n, ζ_n - kružna frekvencija i relativno prigušenje tona n. U programu: dostupna za n ≥ 2, u ručnom režimu i uz dejstva iz Modela; koriste se svi tonovi, a svaka jednačina tona integriše se izabranim postupkom (Newmark ili MCR). Od Φ̃^T·C·Φ̃ zadržavaju se dijagonalni članovi, pa je postupak tačan samo za klasično C (program upozorava kada prigušenje nije klasično).

Izvori: [1] pogl. 12, [3] §9.3.

Rejlijevo prigušenje

a₀, a₁ - koeficijenti koji se dobijaju rešavanjem ove jednačine za dva izabrana tona i i j sa zadatim ζ_i, ζ_j (za n = 1 zadaje se jedno ζ). Umesto toga matrica C se može zadati ručno.

Izvori: [1] §11.4, [3] §9.3.3.

Kretanje tla i korekcija bazne linije

u_g, u̇_g - pomeranje i brzina tla (trapezno pravilo), u_g0, u̇_g0 - početni uslovi kretanja tla. Korekcija drifta brzine oduzima od u̇_g linearni član c₁·t (isto što i oduzimanje srednje vrednosti a_g), a korekcija drifta brzine i pomeranja još i parabolični član c₂·(t²/2 - T_z·t/2) (T_z - trajanje zapisa). Korekcija utiče samo na apsolutna pomeranja i brzine; relativni odgovor ne menja.

Izvori: [5], [6].

Equation of motion and excitation

M, C, K - mass, damping and stiffness matrices; u, u̇, ü - displacement, velocity and acceleration vectors (relative to the ground with an accelerogram); F(t) - load vector; ι - influence vector (with the model from the Model tab ι_j = 1 for the checked DOFs and 0 for the others; with manually entered matrices numerical values ι_j, default 1); a_g - ground acceleration. With Model load cases the static and dynamic parts are superposed. Assumptions: linear system (M, C, K independent of time and response), uniform support excitation, dynamic loads only on DOFs that carry mass.

Sources: [1] ch. 9, §13.1, [2].

Direct numerical integration (DNI): Newmark and the central difference method

What and when: the equations M·ü + C·u̇ + K·u = F are solved together at every step Δt; valid also for non-classical damping.

i - step number, γ, β - method parameters. Newmark with constant average acceleration is unconditionally stable and has no numerical damping (the amplitude does not decay), but it elongates the period of vibration: T̄/T = (ω·Δt/2)/arctan(ω·Δt/2) ≈ 1 + (ω·Δt)²/12, where T̄ is the numerically obtained period, T the exact natural period, ω = 2π/T the circular frequency of the mode and Δt the integration time step (for Δt = T/10 the elongation is about 3.2 %, for Δt = T/20 about 0.8 %; [1] §5.5, Fig. 5.5.2, [3] §9.4.3, Fig. 9.6). The central difference method (CDM) shortens the period, T̄/T ≈ 1 - (ω·Δt)²/24, and is conditionally stable: Δt < T_n,min/π (T_n,min - shortest natural period); the program checks this and, if the condition is not met, the analysis is not run. The step Δt is entered by the user and is not changed by the program; for accuracy Δt ≤ T/10 is usually taken for the highest mode that matters for the response.

Sources: [1] §5.3, §5.4, [3] §9.2.1, §9.2.3, [4].

Numerical integration step and load record step

The numerical integration step Δt is chosen by the user. The load record (force F(t) or ground acceleration a_g(t)) has its own step Δt_rec: from the field “Load record step” for one-column files (empty field: Δt_rec = Δt), and for a two-column file (t, value) from the time column, which must increase with a constant step (otherwise the file is not loaded). The integration always uses a constant step h, and the record is converted to the analysis instants t_j = j·h by linear interpolation:

m - integer (Δt_rec = m·h; ⌈x⌉ - the smallest integer ≥ x), F_i - the i-th record value, t_i - its instant. For Δt ≤ Δt_rec every record value is an analysis instant, so characteristic (e.g. peak) load values are not skipped, and between them the load varies linearly - the same assumption the integration methods use (excitation linear within a step, [1] §5.2). When h differs from the entered Δt, a note with the adopted step is shown with the results. For Δt > Δt_rec the record values between the analysis instants are skipped; a warning with the largest record value and the largest value at the analysis instants is shown with the results. If the records have different steps, h is aligned with the smallest; points of a record whose step is not an integer multiple of that step are interpolated (note with the results). After the end of the record the load is 0. When Δt_rec = Δt, the record is copied literally, so the results are the same as without alignment.

Recommendation: Δt ≤ T/10 (T - the shortest period of the modes contributing significantly to the response; Newmark elongates the period by ≈ (ω·Δt)²/12, ω = 2π/T, about 3 % for Δt = T/10) and Δt ≤ Δt_rec; for the central difference method also Δt < T_min/π (T_min - the shortest natural period of the system).

Sources: [1] §5.2, §5.5.

Modal superposition

q_n - modal coordinate, φ̃_n - mass-normalised mode shape, ω_n, ζ_n - circular frequency and damping ratio of mode n. In the program: for n ≥ 2, in manual mode and with Model load cases; all modes are used, and each modal equation is integrated by the chosen method (Newmark or CDM). Only the diagonal terms of Φ̃^T·C·Φ̃ are kept, so the method is exact only for a classical C (a warning is shown when damping is non-classical).

Sources: [1] ch. 12, [3] §9.3.

Rayleigh damping

a₀, a₁ - coefficients obtained by solving this equation for two chosen modes i and j with the prescribed ζ_i, ζ_j (one ζ for n = 1). Alternatively, the matrix C can be entered manually.

Sources: [1] §11.4, [3] §9.3.3.

Ground motion and baseline correction

u_g, u̇_g - ground displacement and velocity (trapezoidal rule), u_g0, u̇_g0 - initial conditions of the ground motion. The velocity drift correction subtracts a linear term c₁·t from u̇_g (the same as removing the mean of a_g), and the velocity and displacement drift correction also a parabolic term c₂·(t²/2 - T_r·t/2) (T_r - record duration). The correction affects only absolute displacements and velocities; the relative response is unchanged.

Sources: [5], [6].

Elastični i projektni spektar (EN 1998-1)

a_g = γ_I·a_gR - projektno ubrzanje tla tipa A (a_gR - referentno maksimalno ubrzanje, zadaje se a_gR/g; γ_I - faktor značaja, Tab. 4.3); S - faktor tla; T_B, T_C, T_D - granični periodi (Tab. 3.2 i 3.3 za tip tla A-E; Tip 3 - spektar nacionalnog priloga SRPS EN 1998-1/NA); q - faktor ponašanja; β = 0,2 - faktor donje granice; η = √(10/(5 + ξ)) ≥ 0,55 (η - faktor korekcije prigušenja, ξ = 100·ζ - relativno viskozno prigušenje u %). Vertikalni spektar (§3.2.2.3) ima faktor 3,0 umesto 2,5, sopstvene a_vg/a_g, T_B, T_C, T_D i poseban q. Parametri u tabeli mogu da se menjaju.

Izvori: [1] §3.2.2.2, §3.2.2.3, §3.2.2.5, [2] §6.6.

Modalne veličine i odgovor po tonovima

M*_k = φ̃_k^T·M·φ̃_k - generalisana masa, L_k - faktor pobude, Γ_k - faktor participacije, M_eff,k - efektivna modalna masa (i njen udeo u ukupnoj masi), S_k - seizmičke sile tona k, d_e,k - pomeranja tona k, ι - vektor uticajnih koeficijenata pravca pobude (uz Model automatski po pravcima X, Y, Z).

Odgovor: seizmičke sile i pomeranja po tonovima, a uz Model i presečne sile i reakcije. Pomeranja d_s = q_d·d_e (§4.3.4), gde je d_e pomeranje iz analize sa projektnim spektrom S_d (koji već sadrži redukciju faktorom ponašanja q) i q_d = q; za T* < T_C prikazuje se i q_d = μ_δ = 1 + (q - 1)·T_C/T* (Aneks B, (B.10); T* - period tona sa najvećom efektivnom modalnom masom u pravcu pobude, ne nužno T₁). Posle kombinovanja po tonovima d_s,i = min(q_d·d_e,i, 3·d_el,i), jer po Aneksu B ciljno pomeranje ne mora biti veće od trostrukog elastičnog; d_el - pomeranje iz elastičnog spektra S_e, prikazuje se i radi poređenja.

Izvori: [1] §4.3.3.3, §4.3.4, Aneks B, [2] §13.7, [3].

Kombinovanje tonova i pravaca (SRSS, CQC)

• Kombinovanje po tonovima: SRSS (koren zbira kvadrata) ili CQC (kompletna kvadratna kombinacija, koeficijenti korelacije ρ_jk po Der Kjuregijanu); po pravcima X, Y, Z uvek SRSS (§4.3.3.5); upozorenje ako relevantni tonovi u pravcu pobude (M_eff,k/M_tot ≥ 0,05, odnosno do zbira 0,90, §4.3.3.3.1(3)) nisu nezavisni (T_j > 0,9·T_i, §4.3.3.3.2), a izabran je SRSS.
• Koeficijent korelacije CQC: ρ_jk = 8ζ²·(1 + β_jk)·β_jk^3/2 / [(1 - β_jk²)² + 4ζ²·β_jk·(1 + β_jk)²], gde je β_jk = ω_k/ω_j odnos kružnih frekvencija tonova k i j, a ζ udeo kritičnog prigušenja (isti za sve tonove); ρ_jj = 1, ρ_jk = 1 i za β_jk = 1 (tonovi iste frekvencije), matrica je simetrična (ρ(β) = ρ(1/β)), a za ω ≈ 0 uzima se ρ_jk = 0 van dijagonale. Kombinacija: E_i = √(Σ_jΣ_k ρ_jk·E_j,i·E_k,i), E_k,i - doprinos tona k veličini i.

Pretpostavke: SRSS važi za dobro razdvojene tonove; CQC uzima u obzir korelaciju bliskih tonova. Isto ζ za sve tonove.

Izvori: [1] §4.3.3.3.2, §4.3.3.5, [2] §13.7, [4], [5].

Elastic and design spectrum (EN 1998-1)

a_g = γ_I·a_gR - design ground acceleration on type A ground (a_gR - reference peak ground acceleration, entered as a_gR/g; γ_I - importance factor, Table 4.3); S - soil factor; T_B, T_C, T_D - corner periods (Tables 3.2 and 3.3 for ground types A-E; Type 3 - National Annex spectrum SRPS EN 1998-1/NA); q - behaviour factor; β = 0.2 - lower bound factor; η = √(10/(5 + ξ)) ≥ 0.55 (η - damping correction factor, ξ = 100·ζ - viscous damping ratio in %). The vertical spectrum (§3.2.2.3) has the factor 3.0 instead of 2.5, its own a_vg/a_g, T_B, T_C, T_D and a separate q. The parameters in the table can be edited.

Sources: [1] §3.2.2.2, §3.2.2.3, §3.2.2.5, [2] §6.6.

Modal quantities and modal response

M*_k = φ̃_k^T·M·φ̃_k - generalised mass, L_k - excitation factor, Γ_k - participation factor, M_eff,k - effective modal mass (and its ratio to the total mass), S_k - seismic forces of mode k, d_e,k - displacements of mode k, ι - influence vector of the excitation direction (set automatically per X, Y, Z direction with the Model).

Response: modal seismic forces and displacements; displacements d_s = q_d·d_e (d_e - displacement from a linear analysis with the design spectrum S_d, which is already reduced by the behaviour factor q, before multiplication by q_d; q_d = q, and for T* < T_C q_d = μ_δ = 1 + (q - 1)·T_C/T* is also shown, T* - period of the dominant mode, i.e. the mode with the largest effective modal mass in the excitation direction (shown in the results as T* = T_k with the number k of that mode, not necessarily T₁); §4.3.4, Annex B, (B.10); after the modal combination d_s,i = min(q_d·d_e,i, 3·d_el,i), since by Annex B the target displacement need not exceed 3 times the elastic one - the bound is permitted, a result without it would be conservative) and, for comparison, d_el - displacement from the elastic spectrum S_e; with the Model, internal forces and reactions.

Sources: [1] §4.3.3.3, §4.3.4, Annex B, [2] §13.7, [3].

Combination of modes and directions (SRSS, CQC)

• Modal combination: SRSS (square root of the sum of squares) or CQC (complete quadratic combination, correlation coefficients ρ_jk after Der Kiureghian); over the X, Y, Z directions always SRSS (§4.3.3.5); a warning if the relevant modes in the excitation direction (M_eff,k/M_tot ≥ 0.05, or up to a sum of 0.90, §4.3.3.3.1(3)) are not independent (T_j > 0.9·T_i, §4.3.3.3.2) and SRSS is selected.
• CQC correlation coefficient: ρ_jk = 8ζ²·(1 + β_jk)·β_jk^3/2 / [(1 - β_jk²)² + 4ζ²·β_jk·(1 + β_jk)²], where β_jk = ω_k/ω_j is the ratio of the circular frequencies of modes k and j, and ζ the fraction of critical damping (the same for all modes); ρ_jj = 1, ρ_jk = 1 also for β_jk = 1 (modes of equal frequency), the matrix is symmetric (ρ(β) = ρ(1/β)), and for ω ≈ 0 ρ_jk = 0 off the diagonal. Combination: E_i = √(Σ_jΣ_k ρ_jk·E_j,i·E_k,i), E_k,i - contribution of mode k to quantity i.

Assumptions: SRSS holds for well-separated modes; CQC accounts for the correlation of closely spaced modes. One ζ for all modes.

Sources: [1] §4.3.3.3.2, §4.3.3.5, [2] §13.7, [4], [5].

Jednačina kretanja i spektri

u, u̇, ü - relativno pomeranje, brzina i ubrzanje; ω = 2π/T - kružna frekvencija, T - period; ζ - relativno prigušenje; S_d, S_v, S_a - spektri relativnog pomeranja, brzine i ubrzanja, S_aa - spektar apsolutnog ubrzanja, S_pv, S_pa - pseudo-spektri brzine i ubrzanja. Kontrola: S_aa = max|ω²·u + 2ζω·u̇| (uslov ravnoteže).

Izvori: [1] §6.6, §6.12, [2].

Postupak Nigam-Jennings (tačno rešenje uz linearnu interpolaciju pobude), Newmark i metoda centralnih razlika

Postupak Nigam-Jennings (podrazumevano; [1] §5.2, tab. 5.2.1, [2]): ubrzanje tla a_g menja se unutar koraka Δt linearno, a jednačina kretanja se za takvu pobudu rešava tačno, rekurentnim formulama u_i+1 = A·u_i + B·u̇_i + C·p_i + D·p_i+1, u̇_i+1 = A′·u_i + B′·u̇_i + C′·p_i + D′·p_i+1 (i - redni broj trenutka proračuna; p = -a_g - pobuda po jedinici mase; A, B, C, D, A′, B′, C′, D′ - koeficijenti koji zavise od ω, ζ i Δt i važe i za ζ = 0). Rešenje je tačno u trenucima proračuna za svaku veličinu koraka; podrazumevani korak Δt ≤ T/100 služi samo da maksimum odgovora bude dobro obuhvaćen. Maksimum može nastupiti i između dva trenutka proračuna (relativna brzina sistema dugog perioda, ≈ -u̇_g, i apsolutno ubrzanje sadrže visoke frekvencije kretanja tla), pa program u svakom koraku u kome izvod posmatrane veličine menja znak (u̇ za S_d, ü za S_v, dü/dt za S_a, dü_aps/dt za S_aa) nalazi trenutak u kome je izvod jednak nuli i u njemu računa stanje istim tačnim rešenjem; preostala relativna greška maksimuma je ≤ 2·10^-6.

Newmark (konstantno srednje ubrzanje, γ = 1/2, β = 1/4) i metoda centralnih razlika (MCR) ponuđeni su radi poređenja; podrazumevani korak integracije za njih je Δt ≤ T/20, pa je uslov stabilnosti MCR (Δt < T/π) ispunjen (periodi za koje uz zadati korak uslov nije ispunjen računaju se postupkom Newmark, uz napomenu). Oba postupka menjaju period oscilovanja: Newmark ga produžava, T̄/T = (ω·Δt/2)/arctg(ω·Δt/2) ≈ 1 + (ω·Δt)²/12, a MCR ga skraćuje, T̄/T = (ω·Δt/2)/arcsin(ω·Δt/2) ≈ 1 - (ω·Δt)²/24, gde je T̄ numerički dobijen period, T tačan sopstveni period, ω = 2π/T kružna frekvencija oscilatora i Δt korak integracije; za Δt = T/20 to je oko +0,8 %, odnosno -0,4 % ([1] §5.5, sl. 5.5.2; [4] §9.4.3, sl. 9.6). Greška perioda izaziva faznu grešku koja se bez prigušenja nagomilava iz ciklusa u ciklus, pa program uz rezultate ispisuje najveće izmereno odstupanje S_d od rešenja postupkom Nigam-Jennings sa istim korakom integracije i preporuku.

Korak numeričke integracije zadaje korisnik na jedan od dva načina; za svaki period T integriše se ceo akcelerogram. (1) Uslov Δt ≤ T/n: za svaki period T usvaja se Δt = Δt_zapis/m, m = ⌈n·Δt_zapis/T⌉, pa je Δt ≤ T/n i Δt ≤ Δt_zapis (prazno polje: preporučena vrednost n = 100 za Nigam-Jennings, n = 20 za Newmark i MCR; dozvoljeno je n ≥ 1 - za n = 1 je Δt = Δt_zapis za sve periode T ≥ Δt_zapis - ali se male vrednosti n ne preporučuju). (2) Zadat korak Δt, isti za sve periode: za Δt ≤ Δt_zapis usklađuje se sa korakom zapisa, m = ⌈Δt_zapis/Δt⌉, Δt = Δt_zapis/m; za Δt > Δt_zapis vrednosti zapisa između trenutaka proračuna se preskaču, pa karakteristične snimljene vrednosti mogu biti izostavljene (upozorenje uz rezultate). Uz usklađen korak trenutak svake snimljene vrednosti je i trenutak proračuna, a između snimljenih vrednosti ubrzanje a_g se linearno interpolira (m - ceo broj, ⌈x⌉ - najmanji ceo broj ≥ x). Manji n povećava grešku perioda kod postupaka Newmark i MCR, a za periode za koje je Δt ≥ T/π MCR nije stabilna, pa se ti periodi računaju postupkom Newmark (napomena uz rezultate).

Posle kraja zapisa a_g linearno opada do nule tokom jednog koraka Δt_zapis (zapis je dopunjen jednom nulom, pa pobuda ne zavisi od koraka integracije), a zatim se računa još jedan period T slobodnih oscilacija (a_g = 0), jer maksimum može nastupiti i posle kraja zapisa. Za vrlo kratke periode (T < Δt_zapis/100) sistem se smatra krutim: relativno kretanje je nula, a S_aa = S_pa = max|a_g|. Uslov Δt ≤ T/n se poštuje do granice obima proračuna od 10⁷ koraka integracije za jedan period (ceo akcelerogram i jedan period slobodnih oscilacija); kada ta granica smanji broj podkoraka m, uz rezultate se ispisuju traženi i primenjeni korak.

Izvori: [1] §5.2-5.4, [2], [3].

Equation of motion and spectra

u, u̇, ü - relative displacement, velocity and acceleration; ω = 2π/T - circular frequency, T - period; ζ - damping ratio; S_d, S_v, S_a - spectra of relative displacement, velocity and acceleration, S_aa - absolute acceleration spectrum, S_pv, S_pa - pseudo-velocity and pseudo-acceleration spectra. Check: S_aa = max|ω²·u + 2ζω·u̇| (equilibrium).

Sources: [1] §6.6, §6.12, [2].

Nigam-Jennings method (exact solution for linearly interpolated excitation), Newmark and the central difference method

Nigam-Jennings method (a method based on interpolation of the excitation; [1] §5.2, Table 5.2.1, [2]; default): the ground acceleration a_g is assumed to vary linearly within each integration time step Δt, and the equation of motion is solved exactly for this excitation. The state at the end of the step follows from the recurrence formulas u_i+1 = A·u_i + B·u̇_i + C·p_i + D·p_i+1, u̇_i+1 = A′·u_i + B′·u̇_i + C′·p_i + D′·p_i+1 (i - number of the computed time instant; p = -a_g - excitation per unit mass; A, B, C, D, A′, B′, C′, D′ - coefficients depending on ω, ζ and Δt, [1] Table 5.2.1; valid also for ζ = 0). The solution is exact at the computed time instants regardless of the step size, so the default step Δt = Δt_rec/m ≤ T/100 (m - integer) is adopted only to capture the peak response well. If the peak were sought only at the computed instants, the error for purely harmonic motion of period T would be at most 1 - cos(π/100) ≈ 0.05 %; however, the relative velocity of long-period systems (≈ -u̇_g, where u̇_g is the ground velocity) and the absolute acceleration also contain high frequencies of the ground motion, so their peak falls between two computed time instants (error down to -1.6 % for S_v with Δt_rec = 0.02). Therefore, in every step in which the derivative of a quantity changes sign (u̇ for S_d, ü for S_v, dü/dt for S_a, dü_abs/dt for S_aa) the program finds the zero of the derivative, i.e. the instant at which the derivative equals zero, and computes the state at that instant with the same exact solution, so the remaining peak error is negligible (verified: ≤ 2·10^-6 relative).

Newmark (constant average acceleration, γ = 1/2, β = 1/4) and the central difference method (CDM) are offered for comparison; their default integration time step is Δt ≤ T/20, so the CDM stability condition (Δt < T/π) is met (periods for which it is not met with the given step are computed with Newmark, with a note). Both methods change the period of vibration: Newmark elongates it, T̄/T = (ω·Δt/2)/arctan(ω·Δt/2) ≈ 1 + (ω·Δt)²/12, and the CDM shortens it, T̄/T = (ω·Δt/2)/arcsin(ω·Δt/2) ≈ 1 - (ω·Δt)²/24, where T̄ is the numerically obtained period, T the exact natural period, ω = 2π/T the circular frequency of the oscillator and Δt the integration time step; for Δt = T/20 this is about +0.8 % and -0.4 % ([1] §5.5, Fig. 5.5.2; [4] §9.4.3, Fig. 9.6). The period error causes a phase error which, without damping, accumulates from cycle to cycle, so the program prints with the results the measured largest deviation of S_d from the Nigam-Jennings solution with the same integration time step and a recommendation.

The numerical integration step is set by the user in one of two ways; for every period T the whole accelerogram is integrated. (1) Condition Δt ≤ T/n: for every period T the step Δt = Δt_rec/m, m = ⌈n·Δt_rec/T⌉ is adopted, so Δt ≤ T/n and Δt ≤ Δt_rec (empty field: the recommended value n = 100 for Nigam-Jennings, n = 20 for Newmark and CDM; n ≥ 1 is allowed - n = 1 gives Δt = Δt_rec for all periods T ≥ Δt_rec - but a small n is not recommended). (2) A given step Δt, the same for all periods: for Δt ≤ Δt_rec it is aligned with the record step, m = ⌈Δt_rec/Δt⌉, Δt = Δt_rec/m; for Δt > Δt_rec the record values between the analysis instants are skipped, so characteristic recorded values may be omitted (warning with the results). With an aligned step every recorded value is an analysis instant, and a_g is interpolated linearly between them (m - integer, ⌈x⌉ - the smallest integer ≥ x). A smaller n increases the period error of the Newmark and CDM methods, and for periods with Δt ≥ T/π CDM is unstable, so those periods are computed with Newmark (note with the results).

After the end of the record a_g decreases linearly to zero over one step Δt_rec (the record is padded with one zero, so the excitation does not depend on the integration time step), and then one more period T of free vibration (a_g = 0) is computed, because the peak may occur after the end of the record. Very short periods (T < Δt_rec/100) are treated as a rigid system: the relative motion is zero and S_aa = S_pa = max|a_g|. The condition Δt ≤ T/n is honoured up to the computational limit of 10⁷ integration steps for one period (the whole accelerogram and one period of free vibration); when this limit reduces the number of sub-steps m, the requested and the applied step are printed with the results.

Sources: [1] §5.2-5.4, [2], [3].

Ustaljeni odgovor sistema sa TMD-om

X = [X₁, X₂]^T - kompleksne amplitude primarne mase i TMD-a, i - imaginarna jedinica, p - kružna frekvencija pobude; F̂ = [F₀, 0]^T za silu, a F̂ = [-m·A_g, -m_d·A_g]^T za ubrzanje tla (odgovor u odnosu na tlo); hod TMD-a (pomeranje mase TMD-a u odnosu na primarnu masu) je |X₂ - X₁|. Pretpostavke: linearan sistem, samo ustaljeno harmonijsko stanje.

Odgovor se dobija tačno, rešavanjem ovog kompleksnog sistema u frekventnom domenu - bez numeričke integracije jednačina kretanja. Vremenski grafici u(t) crtaju se iz kompleksnih amplituda u trenucima t_i = (i - 1)·Δt, pa polje Δt određuje samo gustinu tačaka grafika, a ne tačnost rezultata.

Izvori: [1] čl. 24, [2].

Optimizacija i opcije 1-8

• Optimizacija 1-7: najmanja amplituda primarne mase A_m pri zadatoj kružnoj frekvenciji pobude p; 8 (minimaks): najmanji maksimum A_m(ω) u opsegu ω = α₈·ω₀ (N tačaka), uz F₀ = const, F₀(ω) = m_e·e·ω² (m_e - ekscentrična masa, e - ekscentricitet) ili F₀(ω) iz txt fajla (promenljiva amplituda samo za pobudu silom).
• Granice: m_d, k_d, c_d (mogu se popuniti iz μ = m_d/m, α = ω_d/ω₀ i ζ_d,min-ζ_d,max), ograničenje ζ_d,min-ζ_d,max; početna tačka x₀.
• Kontrole: (opt. 1-7) SDOF rešenje za m_d, c_d, k_d → 0; (opt. 8, F₀ = const) poređenje sa rešenjem Den Hartoga (sila), odnosno Vorbertona (tlo): f = 1/(1 + μ), ζ_d = √(3μ/(8(1 + μ))), odnosno f = √(1 - μ/2)/(1 + μ), ζ_d = √(3μ/(8(1 + μ)(1 - μ/2))) (f = ω_d/ω₀; ω_d, ω₀ - kružne frekvencije TMD-a i primarnog sistema; μ = m_d/m; m, m_d - njihove mase). Kontrola 1 (čekboks) samo isključuje prigušenje primarnog sistema (c = 0, ζ_d,min = 0) radi poređenja; nezavisno od nje, „Kontrola tačnosti modela” izvodi se uvek (v. ispod).

Postupci: funkcija cilja je A_m (opcije 1-7) ili max A_m(ω) u opsegu (opcija 8); gradijent se ne koristi. Naziv opcije navodi postupak u programu, a u zagradi MATLAB ekvivalent:

• 1-4 (isti postupak): pretraga po mreži od N³ tačaka (N vrednosti svakog od parametara m_d, c_d, k_d u njihovim granicama), pa Nelder-Mead iz najbolje tačke mreže i iz x₀ + završno poboljšanje rešenja;
• 5 i 6 (isti postupak): više nasumičnih početnih tačaka (prva je x₀) + Nelder-Mead + završno poboljšanje najboljeg rešenja (tabela pokretanja prikazuje rešenja pre završnog poboljšanja, a red ⭐ konačni optimum);
• 7 (hibridna): Nelder-Mead sa malim brojem iteracija iz nasumičnih početnih tačaka, pa iz najboljeg rešenja Nelder-Mead + završno poboljšanje + još jedan Nelder-Mead; izbor 1 / 2 (oznake za poređenje sa MATLAB-om: ga / simulannealbnd) ne menja postupak - genetski algoritam i simulirano kaljenje se ne izvršavaju;
• 8 (minimaks): pretraga po mreži od N³ tačaka, pri čemu se max A_m(ω) računa na ređem skupu frekvencija; zatim Nelder-Mead iz 8 najboljih tačaka mreže i iz x₀, Pattern Search i dodatne početne tačke (među njima rešenje Den Hartoga, odnosno Vorbertona).

Završno poboljšanje rešenja: Nelder-Mead sa ponovnim pokretanjem iz dobijenog rešenja + Pattern Search. Nelder-Mead pomera simpleks od n + 1 tačaka (refleksija, širenje, skupljanje) bez izvoda funkcije; Pattern Search ispituje tačke u koordinatnim pravcima ±e_j i smanjuje korak kad nijedna nije bolja. Van granica m_d, c_d, k_d funkcija cilja dobija vrlo veliku vrednost, a ograničenje ζ_d,min-ζ_d,max ulazi kaznenim članom J = A_m + A_ref·10¹²·Δζ_d² (Δζ_d - prekoračenje granice, A_ref - referentna amplituda).

Pri jednoj frekvenciji (opcije 1-7) optimum je po pravilu na donjoj granici ζ_d,min: neprigušen TMD podešen na frekvenciju pobude poništava odgovor primarne mase (antirezonanca, k_d/m_d = p²), pa je ζ_d,min projektna odluka, a ne rezultat.

Kontrola tačnosti modela (uvek, opcije 1-7): A_m iz 2-DOF rešenja sa m_d, c_d, k_d → 0 poredi se sa SDOF rešenjem A_m = λ·F₀/k (sila), odnosno A_m = λ·m·A_g/k (ubrzanje tla, pomeranje u odnosu na tlo); λ = 1/√((1 - η²)² + (2ζη)²) - dinamički koeficijent, η = p/ω₀, ω₀ = √(k/m), ζ = c/(2·m·ω₀) - relativno prigušenje primarnog sistema; m, k, c - masa, krutost i koeficijent prigušenja primarnog sistema, F₀ - amplituda sile, A_g - amplituda ubrzanja tla, p - kružna frekvencija pobude; relativna razlika mora biti manja od 10^-9. Kontrola 1 (isključen čekboks „Uključeno prigušenje primarnog sistema”) daje c = 0 i ζ_d,min = 0; granice c_d i ζ_d,max ostaju kako su upisane. Podešavanje za proveru: c = 0, ζ_d,min = 0, ζ_d,max = 0,05.

Den Hartogovo rešenje u izvornom obliku [1, str. 130-131]: (c_d/c_c)² = 3μ/(8(1 + μ)³), c_c = 2·m_d·ω₀ - srednja vrednost dva uslova horizontalne tangente krive A_m(ω) u fiksnim tačkama P i Q; uz optimalno podešavanje ω_d = ω₀/(1 + μ) to daje ζ_d = c_d/(2·m_d·ω_d) = √(3μ/(8(1 + μ))). Den Hartogovo i Vorbertonovo rešenje važe samo za konstantnu amplitudu pobude (F₀, odnosno A_g = const).

Izvori: [1] čl. 24, [2], [3], [4], [5].

Steady-state response of the system with a TMD

X = [X₁, X₂]^T - complex amplitudes of the primary mass and the TMD, i - imaginary unit, p - forcing circular frequency; F̂ = [F₀, 0]^T for a force and F̂ = [-m·A_g, -m_d·A_g]^T for a ground acceleration (response relative to the ground); the TMD stroke (displacement of the TMD mass relative to the primary mass) is |X₂ - X₁|. Assumptions: linear system, steady harmonic state only.

The response is obtained exactly, by solving this complex system in the frequency domain - without numerical integration of the equations of motion. The time plots u(t) are drawn from the complex amplitudes at the instants t_i = (i - 1)·Δt, so the Δt field only sets the density of plot points, not the accuracy of the results.

Sources: [1] Art. 24, [2].

Optimisation and options 1-8

• Optimisation 1-7: minimum primary-mass amplitude A_m at the given p; 8 (minimax): minimum of the peak A_m(ω) over the band ω = α₈·ω₀ (N points), with F₀ = const, F₀(ω) = m_e·e·ω² (m_e - eccentric mass, e - eccentricity) or F₀(ω) from a txt file (variable amplitude for force excitation only).
• Bounds: m_d, k_d, c_d (can be filled from μ = m_d/m, α = ω_d/ω₀ and ζ_d,min-ζ_d,max), constraint ζ_d,min-ζ_d,max; starting point x₀.
• Checks: (opt. 1-7) SDOF solution for m_d, c_d, k_d → 0; (opt. 8, F₀ = const) comparison with Den Hartog (force) and Warburton (ground): f = 1/(1 + μ), ζ_d = √(3μ/(8(1 + μ))), and f = √(1 - μ/2)/(1 + μ), ζ_d = √(3μ/(8(1 + μ)(1 - μ/2))) (f = ω_d/ω₀; ω_d, ω₀ - circular frequencies of the TMD and the primary system; μ = m_d/m; m, m_d - their masses). Check 1 (checkbox) only removes the damping of the primary system (c = 0, ζ_d,min = 0) for comparison; independently of it, the “Model accuracy check” is always performed (see below).

Procedures: the objective is A_m (options 1-7) or max A_m(ω) over the band (option 8); no gradient is used. Each option name gives the procedure run by the program, with the MATLAB equivalent in brackets:

• 1-4 (same procedure): N³ grid + Nelder-Mead from the best grid point and from x₀ + refinement;
• 5 and 6 (same procedure): multiple random starts (the first from x₀) + Nelder-Mead + refinement of the best (the table of runs shows the solutions before the final refinement, and the ⭐ row the final optimum);
• 7 (hybrid): random starts + quick Nelder-Mead, then from the best Nelder-Mead + refinement + final Nelder-Mead; the genetic algorithm / simulated annealing choice does not change the procedure;
• 8 (minimax): grid search over N³ points, with max A_m(ω) evaluated on a sparser set of frequencies; then Nelder-Mead from the 8 best grid points and from x₀, pattern search and additional starting points (among them the Den Hartog or Warburton solution).

Final refinement: Nelder-Mead restarted from the solution obtained + pattern search. Nelder-Mead moves a simplex of n + 1 points (reflection, expansion, contraction) without derivatives; pattern search polls points in the coordinate directions ±e_j and reduces the step when none is better. Outside the bounds of m_d, c_d, k_d the objective takes a very large value, and the constraint ζ_d,min-ζ_d,max enters by the penalty J = A_m + A_ref·10¹²·Δζ_d² (Δζ_d - violation of the bound, A_ref - reference amplitude).

At a single frequency (options 1-7) the optimum normally lies at the lower bound ζ_d,min: an undamped tuned TMD cancels the primary response (antiresonance, k_d/m_d = p²), so ζ_d,min is a design decision, not a result.

Model accuracy check (always, options 1-7): A_m from the 2-DOF solution with m_d, c_d, k_d → 0 is compared with the SDOF solution A_m = λ·F₀/k (force) or A_m = λ·m·A_g/k (ground acceleration, displacement relative to the ground); λ = 1/√((1 - η²)² + (2ζη)²) - dynamic magnification factor, η = p/ω₀, ω₀ = √(k/m), ζ = c/(2·m·ω₀) - damping ratio of the primary system; m, k, c - mass, stiffness and damping coefficient of the primary system, F₀ - force amplitude, A_g - ground acceleration amplitude, p - excitation angular frequency. Check 1 (checkbox “Primary system damping enabled” cleared) gives c = 0 and ζ_d,min = 0; the bounds of c_d and ζ_d,max stay as entered. Settings for the check: c = 0, ζ_d,min = 0, ζ_d,max = 0.05.

Den Hartog solution in its original form [1, pp. 130-131]: (c_d/c_c)² = 3μ/(8(1 + μ)³), c_c = 2·m_d·ω₀ - the average of the two conditions of a horizontal tangent of A_m(ω) at the fixed points P and Q; with the optimum tuning ω_d = ω₀/(1 + μ) this gives ζ_d = c_d/(2·m_d·ω_d) = √(3μ/(8(1 + μ))). The Den Hartog and Warburton solutions apply only to a constant excitation amplitude (F₀ or A_g = const).

Sources: [1] Art. 24, [2], [3], [4], [5].

Sistem sa TMD-ovima i odgovor u vremenu

Svaki aktivni TMD dodaje jedan stepen slobode (masa m_d) vezan oprugom k_d i prigušivačem c_d za izabrani stepen slobode konstrukcije; TMD-ovi su međusobno nezavisni. Prošireni sistem M·ü + C·u̇ + K·u = F rešava se direktnom numeričkom integracijom (Newmark ili metoda centralnih razlika, kao u tabu Dinamička analiza odgovora), sa korakom numeričke integracije Δt koji zadaje korisnik. Za metodu centralnih razlika program pre proračuna proverava uslov stabilnosti Δt < T_min/π = 2/ω_max (T_min - najkraći sopstveni period, ω_max - najviša kružna frekvencija spregnutog sistema - konstrukcija + TMD-ovi - za najkruće kombinacije parametara iz granica pretrage) i ne izvodi proračun ako uslov nije ispunjen; posle integracije proverava da li odgovor divergira. Newmark (γ = 1/2, β = 1/4) je bezuslovno stabilan, ali produžava period oscilovanja (v. Pomoć, Dinamička analiza odgovora; [1] §5.5). Modalna superpozicija nije ponuđena: matrica prigušenja C = C_s + C_d (C_s - prigušenje konstrukcije, C_d - prigušivači TMD-ova) u opštem slučaju nije klasična (ne važi C·M^-1·K = K·M^-1·C), pa se jednačine tonova ne razdvajaju.

Izvori: [1] §5.3, §5.4, §16.2, [2].

Funkcija cilja (20 opcija)

č - član funkcije cilja, a_č - koeficijent člana (1, β ili α i 1 - α kod kombinacija), w_i - težina stepena slobode i, x_i - pomeranje u, brzina u̇, ubrzanje ü ili međuspratno pomeranje (drift) Δu_i = u_i - u_i-1; mera je efektivna vrednost RMS = √(Σx²/N) (N - broj trenutaka integracije u intervalu), vršna vrednost max|x|, ∫|x|dt ili ∫x²dt na celom intervalu ili na t₁-t₂. Bez TMD-a je J = J₀ - referentna vrednost koja se dobija kad je u svakom uključenom članu odnos sa/bez TMD-a jednak 1: J₀ je zbir koeficijenata članova - J₀ = 1 za funkcije sa jednim članom i za kombinacije sa α i 1 - α (FC-18 do FC-20), a J₀ = 1 + β za FC-5, FC-12 i FC-17 (npr. J₀ = 2 pri β = 1). J < J₀ znači da TMD smanjuje izabranu zbirnu meru; to ne znači da je odgovor manji u svakom trenutku i u svakom stepenu slobode.

Optimizatori (9 opcija)

Particle Swarm Optimization - roj čestica; brzina čestice zavisi od njenog najboljeg i globalno najboljeg položaja (inercija opada od 0,9 do 0,4, c₁ = c₂ = 2). Differential Evolution - mutacija razlikom vektora populacije i binomno ukrštanje (rand/1/bin, F = 0,8, CR = 0,9). Simulated Annealing - prihvata i lošije rešenje sa verovatnoćom exp(-ΔJ/T), uz eksponencijalno hlađenje T. Pattern Search - direktna pretraga po koordinatnim pravcima sa promenljivim korakom. Nelder-Mead - lokalni simpleks postupak; u kombinacijama „… + NM” dorađuje rešenje globalne metode. Globalni postupci su stohastički - ponovljen proračun može dati malo drugačiji optimum. Sekvencijalnom optimizacijom TMD-ovi se optimizuju jedan po jedan, a simultanom svi zajedno.

Sekvencijalna optimizacija: TMD-ovi se dodaju jedan po jedan u zadatom redosledu (permutacija aktivnih DOF-ova); prethodno optimizovani TMD-ovi ostaju u sistemu sa fiksiranim parametrima dok se optimizuje sledeći. U svakom koraku globalni postupak izvodi N_iter,k = round(0,7·N_iter) iteracija (N_iter - upisani „Maks. broj iteracija”; najmanje 1). Obim pretrage: PSw i DE - N_iter generacija populacije veličine N_pop („Veličina populacije”); SA i PSe - N_iter·N_pop koraka. Pattern Search polazi od najbolje od 6 nasumičnih tačaka, korak posle uspeha udvostručuje, a bez uspeha polovi - po prirodi je lokalni postupak. NM kao samostalan postupak polazi iz središta oblasti pretrage, tj. sredine opsega svakog parametra (najviše 3000 iteracija), a dorada „… + NM” iz rešenja globalne metode (najviše 1000 iteracija).

Izvori: [3], [4], [5], [6], [7].

Granice i referentno rešenje

Granice: m_d = μ·M* (M* = φ^T·M·φ/φ_j² - modalna masa izabranog tona normiranog na 1 u stepenu slobode TMD-a, φ_j - vrednost oblika u tom stepenu slobode, μ - odnos masa), ω_d = f·ω (ω - kružna frekvencija tog tona, f - odnos frekvencija), ζ_d. Referentne vrednosti f i ζ_d daju se po Den Hartogu (sila) ili Vorbertonu (akcelerogram) za μ* = m_d/M* (izrazi su u podtabu 2-DOF TMD).

Izvori: [8], [9].

Pobuda i početni uslovi

Pobuda: sile iz txt fajlova, akcelerogram (uz množilac) F = -M·ι·a_g (ι - vektor uticajnih koeficijenata: bez Modela ι = 1 za sve stepene slobode, uz Model 1 za translacije u pravcu akcelerograma koji izaberete, a 0 za ostale; a_g - ubrzanje tla; za dijagonalnu M: F_i = -m_ii·ι_i·a_g) ili dejstva i kombinacije iz Modela; početni uslovi primarnog sistema u₀, u̇₀ i početno stanje TMD-a po izboru: u_TMD(0) = 0, u̇_TMD(0) = 0 (masa TMD-a miruje, hod s = u_TMD - u ima s(0) = -u₀, ṡ(0) = -u̇₀ čvora, opruga TMD-a je na početku izdužena) ili s(0) = 0, ṡ(0) = 0 (u_TMD(0) = u₀, u̇_TMD(0) = u̇₀ - konstrukcija pomerena zajedno sa TMD-om); izbor menja optimum kod slobodnih oscilacija.

System with TMDs and time-history response

Each active TMD adds one DOF (mass m_d) connected by a spring k_d and a damper c_d to the chosen DOF of the structure; the TMDs are independent of each other. The extended system M·ü + C·u̇ + K·u = F is solved by direct numerical integration (Newmark or the central difference method, as in the Dynamic Response Analysis tab) with the time step Δt set by the user. For the central difference method the program checks the stability condition Δt < T_min/π = 2/ω_max before the calculation (T_min - shortest natural period, ω_max - highest circular frequency of the coupled system - structure + TMDs - for the stiffest parameter combinations within the search bounds) and rejects the calculation if it is not satisfied; after the integration it checks whether the response diverges. The Newmark method (constant average acceleration, γ = 1/2, β = 1/4) is unconditionally stable but elongates the period: ΔT/T ≈ (ω·Δt)²/12 (ω = 2π/T - circular frequency of the mode, T - its period); for Δt = T/20 this is about 0.8 % [1, §5.5]. Modal superposition is not offered: the damping matrix C = C_s + C_d (C_s - structural damping, C_d - TMD dampers) is in general not classical (C·M^-1·K = K·M^-1·C does not hold), so the modal equations do not decouple.

Sources: [1] §5.3, §5.4, §16.2, [2].

Objective function (20 options)

c - term of the objective, a_c - term coefficient (1, β, or α and 1 - α in combinations), w_i - weight of DOF i, x_i - displacement u, velocity u̇, acceleration ü or inter-storey drift Δu_i = u_i - u_i-1; the measure is the RMS value √(Σx²/N) (N - number of integration instants in the window), the peak value max|x|, ∫|x|dt or ∫x²dt over the whole record or over t₁-t₂. Without TMD, J = J₀ - the reference value obtained when the with/without-TMD ratio equals 1 in every included term. J₀ is the sum of the term coefficients - J₀ = 1 for single-term objectives and for the combinations with α and 1 - α (FC-18 to FC-20), and J₀ = 1 + β for FC-5, FC-12 and FC-17 (e.g. J₀ = 2 for β = 1). J < J₀ means that the TMD reduces the selected aggregate measure; it does not mean that every response time history is smaller.

Optimisers (9 options)

Particle Swarm Optimization - a swarm of particles; the particle velocity depends on its own best and the global best position (inertia decreasing from 0.9 to 0.4, c₁ = c₂ = 2). Differential Evolution - mutation by a difference of population vectors and binomial crossover (rand/1/bin, F = 0.8, CR = 0.9). Simulated Annealing - also accepts a worse solution with probability exp(-ΔJ/T), with exponential cooling of T. Pattern Search - direct search in the coordinate directions with a variable step. Nelder-Mead - a local simplex method; in the “… + NM” combinations it refines the solution of the global method. The global methods are stochastic - a repeated run may give a slightly different optimum. Sequential optimisation optimises the TMDs one at a time, simultaneous optimisation all together.

Sequential optimisation: the TMDs are added one by one in the given order (a permutation of the active DOFs); previously optimised TMDs remain in the system with fixed parameters while the next one is optimised. In each step the global method uses N_iter,k = round(0.7·N_iter) iterations (N_iter - the entered “Max. number of iterations”; at least 1). Search effort: PSw and DE - N_iter generations of a population of size N_pop (“Population size”); SA and PSe - N_iter·N_pop steps. Pattern Search starts from the best of 6 random points, doubles the step after a success and halves it otherwise - it is essentially a local method. NM alone starts from the centre of the search domain - the midpoint of each parameter range (at most 3000 iterations), and the “… + NM” refinement from the solution of the global method (at most 1000 iterations).

Sources: [3], [4], [5], [6], [7].

Bounds and reference solution

Bounds: m_d = μ·M* (M* = φ^T·M·φ/φ_j² - modal mass of the selected mode normalised to 1 at the TMD degree of freedom, φ_j - mode value at that DOF, μ - mass ratio), ω_d = f·ω (ω - circular frequency of that mode, f - frequency ratio), ζ_d. Reference values of f and ζ_d are given by Den Hartog (force) or Warburton (accelerogram) for μ* = m_d/M* (expressions in the 2-DOF TMD sub-tab).

Sources: [8], [9].

Excitation and initial conditions

Excitation: forces from txt files, accelerogram (with a scale factor) F = -M·ι·a_g (ι - influence vector: without the Model ι = 1 for every DOF, with the Model 1 for translations in the accelerogram direction you choose and 0 otherwise; a_g - ground acceleration; for a diagonal M: F_i = -m_ii·ι_i·a_g), or Model load cases and combinations; initial conditions of the primary system u₀, u̇₀ and the initial state of the TMD by choice: u_TMD(0) = 0, u̇_TMD(0) = 0 (the TMD mass is at rest, the stroke s = u_TMD - u has s(0) = -u₀, ṡ(0) = -u̇₀ of the node, the TMD spring is initially stretched) or s(0) = 0, ṡ(0) = 0 (u_TMD(0) = u₀, u̇_TMD(0) = u̇₀ - the structure displaced together with the TMD); the choice changes the optimum for free vibration.

Prethodna obrada i FFT

Prethodna obrada signala (opcije, podrazumevano uključene): uklanjanje linearnog trenda i srednje vrednosti, Hanov prozor (smanjuje curenje spektra; bez njega se primenjuje pravougaoni prozor).

Linearni trend je prava a + b·i dobijena metodom najmanjih kvadrata (i - redni broj odbirka, a i b - koeficijenti prave); pošto prava ima slobodni član, posle njenog uklanjanja je Σs_i = 0 (s_i - obrađeni signal), pa posebno uklanjanje srednje vrednosti tada nema uticaja. Bez uklanjanja trenda linearni trend se u spektru javlja kao lažni sadržaj na najnižim frekvencijama (curenje spektra [8, str. 137]). Hanov prozor: glavni luk spektra je dvostruko širi nego kod pravougaonog prozora (vrh je širi), a bočni lukovi, tj. curenje u udaljene linije, znatno su niži. Pošto je amplitudni spektar normiran sa Σw_i, vrh sinusoide čija je frekvencija tačno na spektralnoj liniji ima istu visinu sa prozorom i bez njega; kad frekvencija pada između linija, vrh je sa Hanovim prozorom bliži stvarnoj amplitudi (greška do ≈ 15 % prema ≈ 36 % bez prozora, za signal bez dopune nulama). Kod PSD je vrh na liniji sa Hanovim prozorom 1,5 puta niži (ekvivalentna širina šuma prozora je 1,5 linija), a površina ispod vrha se približno čuva.

Spektri: jednostrani amplitudni spektar, kvadrat amplitudnog spektra |A(f)|² i spektralna gustina snage (PSD) periodogramom i Welch metodom.

Spektri se računaju brzom Furijeovom transformacijom (FFT), Cooley-Tukey postupkom sa osnovom 2 (radix-2), uz dopunjavanje signala nulama do najmanjeg stepena broja 2 koji nije manji od N (N_FFT odbiraka) [8, str. 132, pogl. 6]. Normiranje (jednostrani spektri, bez faktora 2 na f = 0 i na Nikvistovoj liniji):

• amplitudni spektar |A(f_k)| = 2|X_k|/Σw_i (vrh ≈ vršna amplituda sinusoide);
• periodogram G(f_k) = 2|X_k|²/(F_s·Σw_i²); Welch G(f_k) = 2⟨|X_k|²⟩/(F_s·Σw_i²), ⟨⟩ - srednja vrednost po segmentima (w_i - Hanov prozor dužine segmenta, N_FFT → nfft);
• Parsevalova teorema [8, str. 132]: ∫G df = σ² (tačno bez prozora, približno sa Hanovim prozorom), σ² - srednji kvadrat obrađenog signala (varijansa kad je srednja vrednost uklonjena);
• X_k = Σ_i=0^N-1 x_i·w_i·e^{-j2πki/N_FFT} - diskretna Furijeova transformacija (DFT) signala pomnoženog prozorom i dopunjenog nulama; x_i - signal posle prethodne obrade; w_i - težine prozora (Hanov prozor w_i = ½(1 - cos(2πi/(N - 1))), bez prozora w_i = 1); Σw_i - zbir težina; Σw_i² - energija prozora; F_s = 1/Δt.

Pretpostavke: stalna perioda odabiranja Δt - vremenski razmak susednih odbiraka (fajl čija kolona vremena nema stalan korak ili se negde vraća ne učitava se); najviša frekvencija u spektru je Nikvistova frekvencija f_N = 1/(2·Δt) - komponente signala viših frekvencija pojavljuju se u spektru na pogrešnim, nižim frekvencijama (preklapanje spektra, engl. aliasing [8, pogl. 4.2]); razmak susednih spektralnih linija (frekvencijska rezolucija) je Δf = 1/(N_FFT·Δt).

Izvori: [1], [3], [4], [5].

Welch metoda

Dužina segmenta N/K, preklapanje 50 %, ≈ 2K - 1 osrednjenih segmenata, Hanov prozor (N - broj tačaka, K - broj segmenata, ograničen na K ≤ 64 i N/K ≥ 8). Osrednjavanje smanjuje rasipanje (varijansu) procene PSD, po cenu grublje frekvencijske rezolucije, jer je segment K puta kraći od signala.

Izvori: [2], [4], [5].

Dominantne frekvencije

Na sva četiri spektra vrhovi se biraju po istaknutosti (u % od maksimuma) i najmanjem rastojanju Δf; frekvencija vrha određuje se tačnije interpolacijom parabolom kroz logaritme vrednosti tri susedne linije. Tabela daje f i ω dominantnog vrha i do 4 sledeća najveća.

Vrh je linija veća od obe susedne (krajnje linije spektra su kandidati ako su veće od susedne). Prikazuje se ako su mu i visina i istaknutost bar p % od najveće vrednosti spektra; istaknutost je visina vrha iznad veće od dve najniže vrednosti - najniže levo i najniže desno od njega, do prvog višeg vrha ili do kraja spektra. Od dva vrha čije je rastojanje manje od Δf zadržava se veći; Δf treba da bude znatno manje od Nikvistove frekvencije f_N = 1/(2·Δt), inače ostaje najviše jedan vrh.

Arijasov intenzitet

Računa se ako je signal označen kao akcelerogram, iz neobrađenog signala pomnoženog množiocem:

g - ubrzanje Zemljine teže u istom sistemu kao a(t), polje „g”; a(t) - ubrzanje, x(t) - učitani signal bez obrade, c - množilac ubrzanja, T - trajanje zapisa; integral trapeznim pravilom; I_A je u jedinici brzine tog sistema. t₅ i t₉₅ su trenuci u kojima I_A(t) dostiže 5 % i 95 % ukupne vrednosti, a D_5-95 je značajno trajanje - trajanje jakog dela zemljotresnog kretanja tla [7].

Izvori: [6], [7].

Pre-processing and FFT

Pre-processing (options, on by default): linear detrending, mean removal, Hann window (reduces spectral leakage; otherwise rectangular).

The linear trend is the least-squares line a + b·i (i - sample index, a and b - line coefficients); since the line has a constant term, after its removal Σs_i = 0 (s_i - processed signal), so a separate mean removal then has no effect. Without detrending, the linear trend appears in the spectrum as spurious content at the lowest frequencies (spectral leakage). Hann window: the main lobe of the spectrum is twice as wide as with the rectangular window (the peak is wider), and the side lobes, i.e. the leakage into distant lines, are much lower. Since the amplitude spectrum is scaled by Σw_i, the peak of a sinusoid whose frequency lies exactly on a spectral line has the same height with and without the window; when the frequency falls between lines, the Hann peak is closer to the true amplitude (error up to ≈ 15 % versus ≈ 36 % without a window, for a signal without zero padding). For the PSD the on-line peak is 1.5 times lower with the Hann window (equivalent noise bandwidth of 1.5 lines), while the area under the peak is approximately preserved.

Spectra: single-sided amplitude spectrum, squared amplitude spectrum |A(f)|² and power spectral density (PSD) by periodogram and by Welch's method.

Cooley-Tukey radix-2 FFT with zero padding to the next power of two (N_FFT). Scaling (single-sided spectra, without the factor 2 at f = 0 and at the Nyquist line):

• amplitude spectrum |A(f_k)| = 2|X_k|/Σw_i (peak ≈ peak amplitude of a sinusoid);
• periodogram G(f_k) = 2|X_k|²/(F_s·Σw_i²); Welch G(f_k) = 2⟨|X_k|²⟩/(F_s·Σw_i²), ⟨⟩ - mean over segments (w_i - Hann window of segment length, N_FFT → nfft);
• Parseval's theorem: ∫G df = σ² (exact without a window, approximate with the Hann window), σ² - mean square of the processed signal (the variance when the mean is removed);
• X_k = Σ_i=0^N-1 x_i·w_i·e^{-j2πki/N_FFT} - DFT of the signal multiplied by the window and zero-padded; x_i - signal after pre-processing; w_i - window weights (Hann window w_i = ½(1 - cos(2πi/(N - 1))), no window w_i = 1); Σw_i - sum of the weights; Σw_i² - window energy; F_s = 1/Δt.

Assumptions: constant sampling period Δt - the time between neighbouring samples (a file whose time column has no constant step or goes back somewhere is not loaded); the highest frequency in the spectrum is the Nyquist frequency f_N = 1/(2·Δt) - signal components of higher frequencies appear in the spectrum at wrong, lower frequencies (aliasing); the spacing of neighbouring spectral lines (frequency resolution) is Δf = 1/(N_FFT·Δt).

Sources: [1], [3], [4], [5].

Welch's method

Segment length N/K, 50 % overlap, ≈ 2K - 1 averaged segments, Hann window (N - number of points, K - number of segments, limited to K ≤ 64 and N/K ≥ 8). Averaging reduces the scatter (variance) of the PSD estimate at the cost of a coarser frequency resolution, because a segment is K times shorter than the signal.

Sources: [2], [4], [5].

Dominant frequencies

On all four spectra the peaks are picked by prominence (in % of the maximum) and minimum spacing Δf; the peak frequency is refined by a parabola through the logarithms of three adjacent bins. The table gives f and ω of the dominant peak and of up to 4 next largest.

A peak is a line larger than both neighbours (the end lines of the spectrum are candidates if larger than their neighbour). It is shown if both its height and its prominence are at least p % of the largest value of the spectrum; prominence is the height of the peak above the higher of the two lowest values to its left and right, up to the first higher peak or the end of the spectrum. Of peaks closer than Δf the larger one is kept; Δf should be much smaller than the Nyquist frequency f_N = 1/(2·Δt), otherwise at most one peak remains.

Arias intensity

Computed if the signal is flagged as an accelerogram, from the raw signal times the scale factor:

g - gravitational acceleration in the same system as a(t), field “g”; a(t) - acceleration, x(t) - the loaded signal without pre-processing, c - acceleration multiplier, T - record length; trapezoidal rule; I_A is in the velocity unit of that system. t₅ and t₉₅ are the times at which I_A(t) reaches 5 % and 95 % of the total, and D_5-95 is the duration of the strong part of the record.

Sources: [6], [7].