IMPLEMENTACIJA VLAKNASTOG "STUB-GREDA" ELEMENTA U AKADEMSKI CAD SOFTVER-MATRIX 3D FIBER BEAM-COLUMN ELEMENT IMPLEMENTATION IN ACADEMIC CAD SOFTWARE MATRIX 3D

Razvoj nelinearnih analiza konstrukcija, kao i napreci u aseizmičkom projektovanju, doveli su do potrebe za sofisticiranijim konačnim elementima (KE). Iako su 2D i 3D KE nesumnjivo tačniji kada je reč o preciznim nelinearnim analizama, postoje mnoge poteškoće s njihovom upotrebom. Neke od najvažnijih su komplikovana interpretacija rezultata, glomazne numeričke procedure koje zahtevaju dosta vremena, te komplikovana izrada geometrije i mreže samog modela. Kod okvirnih konstrukcija, sve ove prepreke mogu se otkloniti upotrebom vlaknastih linijskih KE. Vlaknasti KE mogu biti podeljeni u one koncipirane na metodi sila (FB) i na metodi deformacija. Iako je velika većina kompjuterskih programa za analizu konstrukcija bazirana na metodi deformacija, razvijeni su i efikasni algoritmi za implemen-


FORMULACIJA VLAKNASTOG STUB-GREDA ELEMENTA
Postoje dve mogućnosti za numeričko opisivanje materijalno-nelinearnog ponašanja u gredama i stubovima: koncentrisana i raspodeljena plastičnost. Ukratko, koncentrisana plastičnost, najčešće modelirana nelinearnim zglobovima ( hindges), podrazumeva dva velika ograničenja. Prvo, precizna pozicija o č ekivanih plastičnih zglobova mora biti poznata a priori, i drugo, nelinearna relacija je definisana kao zavisnost između presečnih sila i deformacija. Ovakvi tipovi zavisnosti su poznati iz eksperimenata, zatim uprosečeni i uprošćeni i potom specificirani u različitim tehničkim propisima kao led by Fillipou [18,[22][23][24] and is afterwards improved by same group to include some other aspects such as shear deformations [8,11], effects of bond deterioration [6], and composite sections [3]. This work established theoretical foundation of OpenSees, software developed on University of California, Berkley, and one of the most popular for academic research in the field of nonlinear seismic analyses. Firstly, formulation is developed for linear geometry/nonlinear materials [8,18], and shortly after Neuenhofer and Filippou [19] presented procedure for linear material/nonlinear geometry. Recently, De Souza presented in [9] extension of their work introducing corotational formulation, thus enabling analysis with large displacements and large rotations. His work can be viewed as a generalisation of the state determination procedures established in [8] and [19]. Different approaches led to other improvements on including shear deformations [16,17], thus enabling to capture interaction between axial, flexural and shear responses with great accuracy. All of the mentioned fibre elements are used primarily for dynamic analysis of concrete structures. Excellent ability of these elements to predict behaviour of composite structures is validated in [13,20]. Thai and Kim [25,26] included stability functions and nonlinear springs, and used FB elements for analysis of steel structure.
In this paper, formulation very similar to one used in [26] will be presented, discussed and tested, also on steel structure examples. This formulation is implemented in academic software for structural analysis -Matrix 3D [28,29], developed on the Faculty of technical sciences in Novi Sad, which is primarily developed for nonlinear analyses. Software is based on finite element method and includes analysis of structures modelled by line, surface and special elements. In this paper only line elements are used -Bernoulli-Euler beam, point supports as elastic springs and point loads. System stiffness matrix Ksys and equivalent node load vector P are assembled by code number method, with transformation to global coordinate system. Solution of unknown displacements for system stiffness matrix and node load vector is found using Gaussian elimination. Matrix 3D has graphical user interface developed, which will be shown in chapter 3.1.

FIBRE BEAM-COLUMN ELEMENT FORMULATION
There are two options for describing nonlinear material behaviour in beams and columns: lumped and distributed plasticity. Briefly, lumped plasticity, most often modelled with hinges, has two major shortcomings. Firstly, precise position of expected plastic hinge must be known a priori, and secondly nonlinear relation is specified as relation between force and deformation. These types of relations are known from experiments, further simplified and specified in various codes such as [21], textbooks etc. More realistic way of performing materially-nonlinear-only analysis (MNO) is through the što je npr. [21], priručnicima i tako dalje. Realističniji način sprovođenja samo-materijalno-nelinearne analize (MNO) jeste kroz konstituisanje zakona materijalne nelinearnosti na nivou napon/deformacija, gde vlaknasti elementi imaju prednost.
First fibre elements were used in the similar way as hinges, positioned in the area of expected plastic hinge, only to represent more realistic force/deformation dependence through integration of stresses. Later the possibility of using only one FE per element became possible. Hence, pre-processing became less timeconsuming and, if the goal is to use one element per member, FB elements are more suitable. Still, FB elements may be scarcely found in commercial CAD software. Main obstacle in the widespread use of this type of elements is the difficulty of integrating nonlinear element state determination in an analysis program based on direct stiffness method [24]. This obstacle is overcome by a state determination procedure, proposed by [22] that iteratively determines the element resisting forces and stiffness matrix, while strictly satisfying element equilibrium and compatibility in each iteration. With this, analysis becomes incremental-iterative, which is discussed in 2.4. Since a pure flexibility formulation is not easily feasible in multi-degree freedom structure model [18], most of the CAD software is displacementbased. This means that linear elements are based on appropriate interpolation functions for member displacements. Conventional frame elements are based on Hermitian polynomials for transverse displacement fields, and Lagrange shape functions for axial displacements. On the contrary, formulation of FB elements is based on interpolation functions for the internal forces. The main difference is that in classical stiffness-based 2node element, curvature distribution is assumed linear (and in real member it becomes highly nonlinear as material softening occurs), while in FB element in the absence of distributed element loads, moment distribution is assumed linear, which is true. As result, only one element per member can be used, with accuracy depending on adopted tolerance of numerical error. On the other hand, displacement-based elements, in addition to numerical error, suffer from error due to the approximate nature of interpolation functions. This is the basic advantage of the FB elements, while others are presented in [15,18].
In software Matrix 3D FB fibre element is implemented, and it is tested in this paper. Basic formulation of this element is based on the mixed method, as it is presented in detail in [24], which for a specific choice of section constitutive relation simplifies to the classical flexibility method. Fundamental assumptions of the element formulation are that deformations are small and plane sections remain plane during the loading history (Bernoulli assumption). With the proposed formulation, equilibrium along the element is always strictly satisfied, and compatibility is satisfied not only at the element ends, but also along the element [24]. Pure derivation of the beam-column element formulation will be omitted in this paper, and can be found in [24], while procedure of calculation implemented in Matrix 3D will be presented.
Process of calculation for a single load increment is achieved on three different levels: structure, element and section level. After calculating structure displacements using global stiffness matrix (1), another set of element degrees of freedom is adopted, and deformations are recalculated using transformation matrix L given in [23].
In a flexibility-based formulation the reference system has to be free of rigid body modes. Therefore, on the element level, there is only 5 degrees of freedom, derived by eliminating rigid body modes, and presented in Figure 1. These element deformations are denoted by q, and corresponding internal forces by Q (3). Generalised displacements given in equation (4) are axial extension and two rotations relative to the chord at each end node. The end rotations and corresponding moments, contained in Q, refer to two arbitrary, orthogonal axes y and z.
Lele matrix is a transformation matrix for the element and it is combination of two transformations: the first is the transformation from the global coordinate system to the local system of element and the other represents the transformation of the obtained displacements to the deformation of element in degrees of freedom obtained by eliminating rigid body displacement.
Kele is found by inverting element flexibility matrix F and afterwards it is further altered by stability functions, which will be explained in the next section. Ignoring effects of axial forces for now, element stiffness matrix is final result obtained from the element level of calculation, and then it is assembled in global stiffness matrix of structure by usual procedure. In order to calculate F by (9), flexibility matrix of each section f(x) is needed, which means that the section level of analyses is initiated. It is not needed to store F in memory between increments of load, but it is so for Kele.
In the first iteration of the first increment, Ei is taken as nominal Young's modulus of material, and in every other case, it is calculated in the fibre level and dependant on the stress state of the considered fibre. In (11), nfib is the number of fibres the section is divided in, z and y are coordinates of fibre from the centroid of section, and A is the area, as represented in Figure 2.
Assumptions of internal force's distribution in element are represented in matrix b. For absence of element loads b is presented in (12). It is obvious that linear distribution of moment throughout the element and a constant axial force are assumed with (12), but other assumptions are possible too.  (12) Na nivou preseka, vektori sila u preseku i deformacija preseka su određeni primenom principa virtualnog rada u preseku, respektivno: On the section level, section force and section deformation vectors, by the principle of virtual work, are respectively: (13) Vektor sila u preseku je kumulativan i računa se prema (14): Section force vector is cumulative and is calculated through (14): Sada je moguće odrediti deformacije preseka prema: Now, section deformations are possible to obtain by: (16) U ovom trenutku proračun se premešta na nivo vlakna, ali nakon završetka petlje po vlaknima, kriterijum konvergencije postavlja se za vektor sila u preseku.
By using selected stress-strain relation, stresses of each fibre are found, after which the calculation of section resisting forces DR is possible (17). As said before, element state determination ( j loop) is completed, and global residual forces PR are determined by: (23) Vektor neuravnoteženih sila u štapu PU računa se kao: Structure unbalanced forces can be determined as: (24) Rezidualne sile elementa se koriste u samouravnotežujućem mehanizmu (7) kroz rezidualne deformacije elementa s (25), koje se dobijaju primenom principa virtuelnih sila.
Vectors denoted with δ are calculated in every iteration of element state determination, while vectors denoted only with ∆ are calculated once per NR iteration.
Hence, algorithm consists of three loops: load increment loop (k loop), NR structure (i) loop, and element state determination one ( j ). Inside the last one there are section and fibre loops, but they are unlikely to bring unbalanced quantities since they depend on values passed by the main 3 loops.
One of the most important aspects of implementation of presented FE is convergence criterion enforced on both the element and structure loop. Element loop is terminated when all of the sections achieve convergence, or when prescribed number of iterations is performed (in Matrix 3D [28,29] it is adopted as 25). If latter is the case, outer (structure) loop is also aborted, and smaller load increment is applied. Structure loop continues NR iterations until the tolerance is achieved. Convergence criteria may be displacement-based, forcebased or energy-based. Cries field in [7] explains that energy criterion can be misleading, while Fillipou in [22] proposes this type of criterion. In Matrix3D, for inner loopj, force criterion is defined, while for structure loop, both energy and force criterion have to be satisfied to proceed. Adopted tolerances are 10 -6 for force, and 10 -12 for energy criterion.

Numerical integration
Since the goal of using one fibre FE per member is to accurately represent true behaviour of a member, two numerical choices are crucial. First one is type of numerical integration, and second one i s number of sections per element. Obviously, the greater number of sections, greater is the accuracy, but this leads to greater cost of analysis. Usually, three to five sections per element are used, giving more than satisfactory accuracy. Even more important is the proper type of integration. This choice dictates positions of the considered sections, as shown in Figure 2. Regarding this particular issue, advantages of Gauss-Lobatto integration rule (26) surpasses the ones of Gauss integration. Even though the latter is derived optimising both the weights and the integration points (stations), former has fixed stations at the ends of interval where the highest bending moments/curvatures are expected. Hence, there are minimum three stations, two endpoints being fixed and arbitrary number of free stations in the middle, symmetrical about the midpoint of the element. The weights of free abscissas are given in (27),and of the endpoints in (28). where n is number of total integration points (stations), and P(x) is a Legendre polynomial.

Material nonlinearity -material models
There are few material models extensively used in numerical researches of steel structures lately. To perform dynamic, or cyclic analysis, hysteretic properties need to be specified. Usually isotropic, kinematic or mixed strain hardening is chosen. All of them are implemented in aforementioned software. Besides, any multi linear material rule can be specified. Special advantage of the flexibility approach is that negative slopes (material failure) can be specified stresses are obtained from strains and not the other way around. However, multi linear type of stress-strain relation may lead to numerical obstacles. Since ksec is calculated with tangent stiffness of each fibre, during collapse analysis (such as pushover, especially close to buckling), if the difference between two adjoining slopes is considerable, convergence will be difficult to achieve. Thus, material models described with smooth curves are preferred. This condition is fulfilled with either of Ramberg-Osgood, Monti-Nuti or Menegotto-Pinto model. Among these models Ramberg-Osgood has an important disadvantage. Since described algorithm requires calculation of stresses from strains, and this model is described with objective function of stress, it means that equation of material law must be solved for each strain. Because there are many fibres, sections, iteration and increments, it becomes vastly time-consuming. Menegotto-Pinto model is improved by Fillipou et al [12] to include isotropic strain hardening. For the analysis performed in this paper, regular model is u sed with the following values: Normalizovani napon σ* dobija se jednačinom (29) u odnosu na normalizovanu dilataciju ε* (30). Oblik krive zavisi od odnosa ojačanja b i od parametra R, koji simulira Baushinger-ov efekat, i dat je u (31). Efekat ova dva parametra na oblik krive napon/dilatacija tokom neizmeničnog opterećenja prikazan je na slici 3. Parametar R kontroliše radijus prelaza između dve grane na slici 3, dok ξ menja R nakon promene smera opterećenja.
Normalised stress σ* is found by equation (29) in relation to normalised strain ε* (30). The curve is controlled by hardening ratio b and a parameter R, which simulates Baushinger effect, given in (31). The effect of those parameters to stress strain curve during load reversals is shown in Figure 3. Parameter R controls smoothness of transition between two branches on Figure 3, while ξ changes R throughout load reversals. where a1 and a2 are material constants and ξ is the absolute value of the plastic strain of the last excursion, normalised with plastic strain εy, as shown in Figure 3. Besides determination of stress for a given strain of fibre, tangent modulus (25) ought to be calculated, and for Menegotto-Pinto model, Ei is given by: Tokom j petlje u procesu određivanja stanja elementa nailazi se na promenu znaka inkrementalne dilatacije, posebno u blizini izvijanja. U takvim situacijama, vrednosti ε0 i εr menjaju vrednost. Ali, da bi se sačuvala glatka kriva istorije napona/dilatacija jednog vlakna, prethodne vrednosti (ε0 ,εr) čuvaju se u memoriji, u slučaju da se znak dva puta promeni tokom iteriranja.

Geometrical nonlinearity
Two effects regarding geometrical nonlinearity are represented through implemented FB formulation: P-∆ and P-δ. P-Δ, also known as geometric nonlinearity, is the nonlinear effect that occurs in each structure with elements subjected to axial forces. This effect introduces secondary forces as a result of the relative displacement of the member ends, ∆. P-δ is the effect that models the secondary forces caused by a local deformation of the member, i.e. elastic line deviation of the chord of the member. Former ( P-Δ) is introduced complementing geometric matrix Kg to the element stiffness matrix in global degrees of freedom (34).It is denoted by Kg, and calculated by (8), where Z0 is zero matrix 3x3. This way of including P-∆ effect in analysis is more accurate, and also more suitable for presented formulation in comparison with calculating only shear due to second navanje samo transverzalnih sila usled efekata drugog reda, i potom ponavljanje proračuna sa ovom silom kao čvornim opterećenjem (kao što je npr. slučaj u programu SAP2000). order effects, and then reiterate with this force as node load (as is used in SAP2000 e.g. Efekat sile koja deluje u pravcu tetive štapa (P-δ) predstavljen je korišćenjem funkcija stabilnosti u matrici krutosti elementa (5x5 -u prirodnim stepenima slobode). Izvođenje ovih funkcija (12) može se naći u [4,5], a prvi ih je uveo James još 1935. godine u radu koji se bavio metodom distribucije momenta. One su izvedene za moment savijanja na jednom kraju poluuklještene grede, uvodeći poprečna pomeranja usled rotacije jednog kraja stuba-grede, i stoga se mogu koristiti samo za kruto vezane elemente. Za proste štapove (npr. elemente sprega), potrebno je koristiti drugačiji set funkcija. Ove funkcije su prikazane u [4].
Effect of axial force acting through the lateral displacement of the beam-columns element (usually called P-small delta), is captured using stability functions in the element stiffness matrix (5x5 -natural degrees of freedom). Derivation of these functions (36) can be found in [4,5] and they were introduced by James in 1935 in a work dealing with moment-distribution method. They are derived for a propped cantilever, introducing lateral displacement due to rotation on one side of beam-column, and hence can be used only for continuous beam-columns. For different boundary conditions (i.e. brace member) different set of functions should be used. These functions are shown in [4].  functions, and hence denoted as sc. Parameter c is termed the carry-over factor, same as in slope-deflection method. The diagrams of s, c and sc are shown in Figure 4, as functions of ρ=P/PE, where PE is Euler's critical force for a simple supported column with length L. With increase of axial compression, s decreases and carry-over factor increases. For ρ=2.046, s becomes 0, and c→∞, with sc finite. The zero value of s represents infinite flexibility at P=2.046PE, a value already recognized to be the first buckling load for a fixed-hinged beam. At values ρ>2.046 moment at hinged end becomes restraining moment, and this explains why the stiffness becomes negative.  -δy(x) umesto 0, i b3,1 postaje -δz(x) [11,20], gde je x prostorna koordinata duž elementa.
not reached, such as full or modified Newton-Raphson (NR), bisection method or combination with NR, etc. In Matrix 3D, full NR is used since the biggest difficulty of the method is calculation of tangent stiffness in each iteration. But since the element stiffness matrices are calculated either way, pure assembling of global matrix is much more beneficial than using modified NR method. In explained algorithm, limit of 20 iterations is imposed as maximum, before load increment is decreased.
latter being more problematic. Namely, the hyper-sphere with its centre in the last converged solution intersects equilibrium path in two points, one being the desired solution, and the other on the equilibrium path portion already determined. Ramm's method is based on finding intersection between desired field of resistance with updated hyper plane. For all mentioned methods load increment is variable, which requires another equation to be solved. Their main differences are in that equation. Both Ramm's and Criesfield methods are implemented in Matrix 3D in conjunction with full NR method, but only Ramm's equation (39) will be shown, since none of the advantages of Criesfield's method is uncovered in present use by the authors. It however should be noted that authors' experience with its use led to conclusion that β factor, which is a value that weights the importance of the contribution of displacements and loads, preferably be taken nonzero, contrary to what many authors suggests.
where: δλ iterative change of load-level parameter (see

VERIFICATION OF FOMULATION
Implementation of above-mentioned formulation is verified for accuracy and efficiency by the comparison of predictions with accessible results and those generated by commercial finite element packages through two examples. First one investigates a column with all end restrains fixed except for longitudinal displacement at the upper end, as shown in Figure 6, with the aim to verify the accuracy and efficiency of the proposed element in capturing the buckling loads of columns with imperfections. Second example is t he Vogel's portal frame, which is usually adopted for calibrating advanced second-order inelastic analysis. Since the frame collapse is determined by the column's inelastic buckling, this example is a good benchmark test for any inelastic formulation [6]. Figure 6. Column clamped on both ends (left), Vogel's frame (centre), residual stress pattern (right)
Both models include initial out-of-plumbness imperfections, modelled by equivalent load in first example, and rearrangement of nodes in the second. Both models include already mentioned out-ofstraightness initial imperfection throughout their length. Even though residual stresses incorporation is not so convenient for both flexibility approach and Menegotto-Pinto model, it may be overcome with few coding improvisations. Residual stress pattern of cross section used for this research is proposed by European Convention for Constructional Steelwork -ECCS as illustrated in Figure 6. Residual stresses are assigned directly to fibres as initial stresses at first increment.
The load-deflection curves of the column obtained by the proposed element ( Figure 7) are compared with those generated by SAP2000 and SeismoStruct. Material m odel used in Matrix 3D (PW) and SeismoStruct (SS)is Menegotto-Pinto with initial modulus of 210GPa and post-elastic hardening of 1%, while bilinear material with the same hardening is used for SAP2000 model. Fibre element of SAP2000 overpredicts buckling load by about 18% if two elements are used, and 6% if 4 elements per member are used. Critical force obtained by Seismo Struct and proposed element are almost identical, while the difference in displacement can be partly explained by the absence of residual stresses in SS model. It should be noted that Seismo Struct has FB element implemented while for SAP2000 model, displacement-based fibre element is used. Design buckling resistance of a column determined by EC3 rules of design (labelled as Nb,Rd) is marked on the graph for comparison. Increasing the number of displacement-based elements in SAP2000 package would lead to better agreement with the results obtained in the SS and PW models. Accuracy depending on the number of elements can be found in [25], but it should be noted that the increase of the element number makes analysis more time-consuming. Figure 7. Load-displacement curve of column Na slici 7 su prikazana poređenja krivih silapomeranje za Vogelov ram, dobijenih predloženim elementom s rezultatima nekoliko autora. Prvo ga je analizirao Vogel 1985. [27], čiji je model plastične zone (PZ) nepoznat autorima, a kasnije Klark 1994. [6], koji je koristio 50 PZ elemenata po stubu i 20 za gredu. Kim i Li [14] korišćenjem komercijalnog softvera Abaqus, modelirali su ram sa 8952 2D-shell elementa. Alvarenga i Silveira [1] 2009. godine koristili su vlaknasti model sa osam elemenata po stubu i šest u gredi. Poređenje sile nosivosti (granične tačke) i maksimalnih horizontalnih pomeranja gornjeg desnog čvora prezentovano je u Tabeli 1.  [27], whose plasticzone (PZ) model is unknown to authors, and later by Clarke in 1994. [6], who used 50 PZ elements per column and 20 for the beam. Kim and Lee [14], using the Abaqus commercial software, modelled the frame with 8952 2D-shell elements. Alvarenga and Silveira [1] in 2009. used fibre model with 8 elements per column and 6 for the beam. Comparison of ultimate load (limit point) and maximum horizontal displacement of upper right joint is presented in Table 1 Kao što se može videti iz Tabele 1, predloženom FB modelu je potreban samo jedan element po štapu kako bi postigao rezultate koji su dobijeni programom Abaqus sa dvodimenzionalnim konačnim elementima korišćenim u [14] ili osam vlaknastih elemenata kao u [1]. Ipak, preciznost u određivanju granične sile je značajno veća nego u određivanju pomeranja. U određenoj meri razlike mogu biti uslovljene različitim post-elastičnim modulom. U svim citiranim radovima korišćen je bilinearni model elastično-savršeno plastično, dok je u ovom radu usvojeno post-elastično ojačanje od 2.05 GPa kako bi se izbegle numeričke poteškoće i zbog čega se javlja donekle kruće ponašanje. Konačno, prezentovani podaci prikazuju zadovoljavajuće poklapanje rezultata dobijenih korišćenjem samo jednog predloženog FB elementa po štapu sa analizama s većim nivoom diskretizacije.

Slika 7. Krive sila-pomeranje stuba
As it can be recognized from Table 1, proposed FB need only one element per member to match the results generated by Abaqus with two-dimensional finite elements used in [14] or 8 fibre elements as in [1]. Again, precision in determination of limit load is considerably greater than for a displacement. Some divergence of results can be attributed to different post-elastic modulus. In all cited analyses elastic-perfectly plastic bilinear model was used, while in present study postelastic hardening of 2.05 GPa is adopted in order to avoid numerical difficulties and hence somewhat stiffer behaviour. Altogether, presented data illustrate considerable agreement between results obtained using only 1 proposed FB element per member and analyses with higher order of discretization.

Result presentation in matrix 3D
Beside mentioned commercial software ( SAP2000 and SiesmoStruct) there are several academic packages that use similar or more advanced (corotational) formulation of fibre element, than the one presented in this paper. Two most widely used are Open Sees, developed on Berkley, which is open source software, and Engineer's Studio, made in Tokyo University. These programs are mainly used for research purposes and rarely for design, since they do not possess graphical interface. Matrix 3D, even though primarily made as academic software, is meant to be practical engineers tool for more refined analyses. Therefore, presentation of results is possible in various manners. Figure 9 depicts displays of stresses in whole structure for a selected load step: full range of stresses in a structure (left), only stresses in observed sections (middle) and spread of plastic strains through the structure (right). The last one is actually a smoothened presentation (between sections) of a stresses higher and lower than two set stress limits, but if the limit is set to yield stress, the output corresponds to plastic strain distribution. Slika 9. Različiti tipovi prikaza rezultata Figure 9. Various displays of results Pored globalnog pregleda rezultata, naponi i dilatacije, kao i radni dijagrami, mogu se sagledati za svako pojedinačno vlakno bilo kog preseka. Ovo je predstavljeno na slici 10 za HEA profil sa usvojenim početnim naponima kao u prikazanom primeru.
Besides the overall inspection of results, stresses and strains, along with stress strain curve, can be acquired for each fibre in any section. This is shown in Figure 10 for a HEA profile with initial stresses distribution as used in aforementioned examples. Figure 10. Results for a single fibre

CONCLUSION AND DISSCUSION
Fibre elements have the essential advantage of representing very accurate behaviour of structures, while requiring simple pre-processing and is capable of delivering results common for engineers, i.e. internal force diagrams. This actually means that no extra time is needed for structural modelling with fibre elements than with commonly used line elements. Beside this, many advantages of more sophisticated and accurate analysis are on the reach with introduction of fibre elements. Such as, for example, introducing real material behaviour, finding a limit load and critical configuration, or even exploring the structural behaviour beyond limit point. And in the case of nonlinear dynamic analysis, elements with distributed plasticity need no presumptive positioning of hinges, and are able to show "true" spreading of plasticity. With the use of proposed formulation, there is a possibility to analyse structure behaviour through the use of a very sophisticated material models, and especially investigate collapse of structure, or real strength of a structure, since material damage fail to induce any numerical problems.
Considering steel structures, complete geometrically and materially nonlinear analysis with initial imperfection (GMNIA) can be performed with this type of element and setting in Matrix 3D. Including advanced heuristic algorithm Matrix 3D uses for load stepping, limit load of the structure can be very precisely determined. This fact aspire the new concept of safety coefficients. In Eurocode 3, design of compressed elements includes many approximations that had to be done in the past. Considering present concept of combining loads and safety coefficients, it is unbiased to say that computer aid is unavoidable. When such, analysis presented can determine safety coefficients more accurately which is the aim of design in general.

Zahvalnost
Ovaj rad je proistekao iz dela istraživačkog projekta TR 36043, koji je podržalo Ministarstvo prosvete i nauke Republike Srbije. user interface. This is the basic advantage of the presented academic software compared to OpenSees, which is considered one of the leading packages for numerical scientific research. Future work that is already in progress will include introducing bond-slip model (for reinforcement and composite sections), including shear in FB formulation, warping, adaptive pushover etc. The present work is the second in a series of papers that will confirm and verify the accuracy and possibilities of newly academic software for testing of structures.