Open access

Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic Load

Written By

Tomasz Kubiak

Submitted: March 11th, 2012 Published: October 24th, 2012

DOI: 10.5772/48961

Chapter metrics overview

4,258 Chapter Downloads

View Full Metrics

1. Introduction

A thin plate or thin-walled constructions are used in the sports industry, automotive, aerospace and civil engineering. As an example of such structural elements snowboard, skis, poles may be mentioned, as well as all kinds of crane girders, structural components of automobiles (car body sheathing or all longitudinal members), aircraft fuselages and wings, supporting structures of the walls and roofs of large halls and warehouses. All the above structures, as well as many others which can be regarded as a thin, exhaust carrying capacity not by exceeding the allowable stresses but by the stability loss. Therefore, not only critical load but also the postbuckling behaviour of thin-walled structures subjected to static and dynamic load is essential knowledge for designers. The use of more accurate mathematical models allows to explore the phenomena occurring after the loss of stability and to describe more precisely their behaviour. Engineers and designers need guidelines to construct as well as quick and easy software to use for analyse the behaviour of thin-walled structures. Therefore, the author of this chapter decided to explore this issue, propose a mathematical model and the method of analysis of orthotropic thin-walled structures subjected to static and dynamic load.

1.1. Static buckling

The buckling and postbuckling of thin-walled structures subjected to static load have been investigated by many authors for more than one hundred years. To the group of precursors of the investigation on the stability of thin-walled structures problem should be included following scientists: Euler [1], Timoshenko [2] and Volmir [3].

This chapter considered the thin plate or thin-walled structures composed of flat plates. Such structures have a various of buckling modes which can differ from one another both in quantitative (e.g., by the number of half-waves) and in qualitative (e.g., by global and local buckling) way.

The stability loss or buckling is a system transition from one equilibrium to another (the bifurcation point), or jump from the stable to the unstable equilibrium path (the limit point). Load resulting in the loss of stability is called the critical load. The behaviour of the structure subjected to load higher than the critical one can be described by a stable (the grow of displacement is caused by increased load) or unstable (displacements grow with decreasing load) postbuckling equilibrium path.

The postbuckling behaviour of the structures depends on their type. For example, the cylindrical shells subjected to axial compression change their equilibrium stage (buckling) by unstable bifurcation point or limit point. Long rods or columns subjected to axial compression have usually a sudden global buckling (bifurcation point of passage to the unstable postbuckling equilibrium path). Thin plates supported on all edges lose their stability having the local buckling mode and the stable postbuckling equilibrium path. Mentioned above type of buckling and postbuckling behaviour for given thin-walled structures are the same for ideal structures as for structures with geometrical imperfection. Columns made ​​of thin prismatic plates could have the local buckling mode, global (flexural, torsional or distorsional) one or coupled. The structures after local buckling are able to sustain further load, because increasing the displacement is only possible by increasing the load value (stable postbuckling equilibrium path), further increasing the load leads to plasticity or reaching the new, this time unstable bifurcation point (global buckling). The dangerous form of stability loss is the interactive buckling (coupled buckling), which usually causes the structure transition to the unstable equilibrium path what leads to the destruction of the structure with load lower than the critical load corresponding to each mode separately. The interaction of different buckling modes occurs when the critical loads corresponding to the different buckling modes are close to each other.

A more comprehensive review of the literature concerning the interactive buckling analysis of an isotropic structure can be found for example in Ali and Sridharan [3], Benito and Sridharan [5], Byskov [6], Koiter and Pignataro [7], Kolakowski [8–10], Manevich [11], Moellmann and Goltermann [12], Pignataro et al. [13], Pignataro and Luongo [14, 15], Sridharan and Ali [16, 17]. The interactive buckling of orthotropic structures can be found for example in [18, 19].

1.2. Dynamic buckling

In literature a quantity of ‘‘pulse intensity’’ [20] or ‘‘pulse velocity’’ [21] is introduced. The analysis of dynamic stability of plates under in-plane pulse loading can be divided into three categories depending on pulse duration and magnitude of its amplitude. For pulses of high intensity the impact phenomenon is observed whereas for pulses of low intensity the problem becomes quasi-static. The phenomenon of dynamic stability and dynamic buckling are often confused with each other. In this chapter the dynamic buckling phenomenon is examined but the concept of dynamic stability is broader and applies also to the stability of motion, which for thin-walled structures can be found for example in [22, 23]. The dynamic buckling occurs when the loading process is of intermediate amplitude and the pulse duration is close to the period of fundamental natural flexural vibrations (in range of milliseconds). In such case the effects of dumping are neglected [24]. Damping neglecting is only possible for problem solved in elastic range [25].

It should be noted that dynamic stability loss may occur only for structures with initial geometric imperfections; therefore the dynamic bifurcation load does not exist. For the ideal structures (without geometrical imperfection) the critical buckling amplitude of pulse loading tends to infinity [26]. The dynamic buckling load should be defined on the basis of the assumed buckling criterion.

The precise mathematical criteria were formulated for structures having unstable postcritical equilibrium path or having limit point [26, 27]. But for the structures having stable postbuckling equilibrium path (thin plate, thin-walled beam-columns with minimal critical load corresponding to local buckling) the precise mathematical criterion have not been defined till now.

Therefore Simitses [27] suggested not to define the dynamic buckling for the structures with stable postbuckling behaviour, but rather it should be defined as a dynamic response to pulse loads.

It is a reason why in world literature a lot of criteria can be found. In the sixties of the twenty century Volmir [28] proposed a criterion for plates subjected to in-plane pulse loading. The Volmir criterion - considered the easiest to use - states that the dynamic critical load corresponds to the amplitude of pulse force (of constant duration) at which the maximum plate deflection is equal to some constant value k (k - one half or one plate thickness)[28].

In many publications the dynamic buckling load is determined on the basis of stability criterion of Budiansky and Hutchinson [26, 29, 30]. However, this criterion was formulated for shell structures but also it can be used for the plate structures [31-34]. Budiansky and Hutchinson noticed that in some range of the amplitude value, the deflection of structures grows more rapidly than in other. Budiansky and Hutchinson formulated the following criterion: Dynamic stability loss occurs when the maximum deflection grows rapidly with the small variation of the load amplitude[26].

In the end of 90’s Ari-Gur and Simonetta [20] analysed laminated plates behaviour under impulse loading and formulated four own criteria of dynamic buckling, two of them of collapse-type conditions. One of them states: Dynamic buckling occurs when a small increase in the pulse intensity causes a decrease in the peak lateral deflection[20].

The failure criterion was proposed by Petry and Fahlbush [34], who suggest that for structures with stable postbuckling equilibrium path the Budiansky-Hutchinson criterion is conservative because it does not take into account load carrying- capacity of the structure.

Based on examples [35] it was noticed that for the thin-walled structures subjected to pulse loading, which lose their stability according to Budiansky-Hutchinson criterion or Volmir criterion, the maximal radius rmax calculated from characteristic root χ=a+jb (where j=1) of Jacoby matrix is equal or greater than unity in complex plane.

Therefore the criterion for thin-walled structures proposed by author [35] can be formulated as follows: Thin-walled structures subjected to pulse loading of finite duration lose their stability even if one characteristic root χ=a+jb of Jacoby matrix find for every time moment from 0 to 1.5Tplies in the complex plane outside the circle with radius equal to unity.

Teter [36, 37] in his works analysed the long columns with longitudinal stiffeners and basing on the phase portrait for dynamic response of these structures defined the following criterion: The dynamic buckling load for the tracing time of solutions has been defined as the minimum value of the pulse load such that phase portrait is an open curve.

The dynamic buckling problem has been well known in the literature for over 50 years and was the subject of numerous works [20, 24, 26-34]. The extensive list of work dealing with dynamic buckling can be found for example in the book edited by Kowal-Michalska [38] or written by Simitses [39] or Grybos [40]. It seems that the analysis of dynamic buckling of thin-walled structures, especially structures with flat walls is not sufficiently investigated. There is a lack of both single-and multimodal analysis of dynamic buckling of columns with complex cross-sections made ​​of thin flat walls. The author of this paper decided to fill this gap presenting a method for the analysis of the local (single mode) and interactive (coupled mode – local and global) buckling of thin-walled structures subjected to pulse loading.

It should be mentioned that the presented method can be used only if the structures are in the elastic range. The case of dynamic buckling in elasto-plastic range including the viscoplastic effect has been investigated by Mania and Kowal-Michalska [25, 41-43]. In world literature it is also possible to find the paper dealing with the dynamic buckling of thin-walled structures subjected to combined load [44]. Czechowski [45] modelled the girder subjected to twist and bending considering only one plate subjected to shear and compression. The general summary showing which parameters have an influence on dynamic buckling of plated structures can be found in [46, 47].


2. Thin orthotropic plate theory

The thin isotropic or orthotropic plates with constant or widthwise variable material properties are considered. The thin-walled beam-columns or girders composed of mentioned above plates are also analysed. In order to taken into account all buckling modes (global, local and their interaction) the plate model was adopted to the analysed structures.

2.1. Basic assumptions

The basic assumption for thin plate are given by Kirchhoff for linear and by von Kármán and Marquerre for nonlinear thin plate theory. They made their assumption for isotropic material; lots of authors extended these assumptions for orthotropic or even for orthotropic multilayer thin plate [18, 48]. The assumptions are as follows:

  • the plate is homogeneous (for example orthotropic homogenisation is made for fibre composite – resin matrix and fibre-reinforcement)

  • the plate is thin – other dimensions (length and width) are at least 10 times higher than plate thickness;

  • the material of the plate subjects to Hooke’s law;

  • the plane stress state is considered for the plate – stress acting in the plate plane dominates the plate behaviour, stress acting in normal to plate plane direction are assumed to be zero;

  • all strains (normal and shear) in plate plane are small compared to unity and they are linear;

  • the strains of the plate to its normal direction are neglected (thickness of the plate do not change after deformation) – this assumption are made according to the Kirchhoff-Love hypothesis;

  • straight lines normal to the mid-surface of the plate remain straight and normal to the mid-surface after deformation

  • there is no interaction in normal direction between layers parallel to middle surface;

  • deflections of the plate can be considerable in terms of nonlinear geometrical relations;

Additionally, it is assumed that principal axes of orthotropy are parallel to the edges of analysed structures (plate, beam, column, beam-column or girder).

2.2. Geometrical equations for thin plate

A plate model has been assumed for a thin plates and thin-walled beam-columns or girders. For easier explanation the plate (Figure 1a) or each i-th strip (Figure 1b) of the plate (or wall of the girder) or each i-th wall of the girder (Figure 1c) are called plate.

To describe the middle surface strains for each plate the following strain tensor have been assumed:

εixm=ui,x+12(wi,x2+ui,x2+vi,x2),  εiym=vi,y+12(wi,y2+ui,y2+vi,y2),γixym=ui,y+vi,x+wi,xwi,y+ui,xui,y+vi,xvi,y,E1

where: ui, vi, wi - displacements parallel to the respective axes xi, yi, zi of the local Cartesian system of co-ordinates, whose plane xiyi coincides with the middle surface of the i-th plate before its buckling (Figure 1).

In the majority of publications devoted to structure stability, the terms(ui,x2+vi,x2), (ui,y2+vi,y2)and (ui,xui,y+vi,xvi,y)are in general neglected for εixm, εiym,  γixymcorrespondingly, in (1) in the strain tensor components.

The change of the bending and twisting curvatures of the middle surface are assumed according to [48, 49] as follows:


The geometrical relationship given by equations (1) and (2) allow to consider both out-of-plane and in-plane bending of the plate.

Figure 1.

Possible models: plates, strips or wall with assumed dimension, coordinate systems and direction of deflections

2.3. Constitutive equations for orthotropy

Let’s consider orthotropic plate with principal axes of ortothropy 1and 2parallel to plate edges (Figure 2).

As same as in previous paragraph let’s consider i-th plate or strip of structures under analysis. The stress – strain relationship for orthotropic plate can be written in following form:




and Ei1, Ei2are Young modulus in longitudinal 1and transverse 2direction respectively, νi12is a Poisson ratio for which strains are in longitudinal direction 1and stress in transverse direction 2, Gi12is a shear modulus (Kirchhoff modulus) in 12plane.

Figure 2.

Plates or walls with principal axes of orthotropy

Young modulus and Poisson ratio occurring in (5) according to Betty-Maxwell theorem or according to symmetry condition of stress tensor should fulfil following relation:


For isotropic plate (wall of beam-columns) the constitutive equations are as follows:


2.4. Generalized sectional forces

Substituting stress-strain relation from previous subchapter, the sectional moments and forces:

  • for i-th isotropic plate or wall of beam-column are expressed by:



  • for i-th orthotropic strip or wall are:



Dix=Eixhi312(1νixyνiyx),   Diy=Eiyhi312(1νixyνiyx),   Dixy=Gixyhi36.E11

2.5. Dynamic equations of stability for thin plate

Differential equations of motion of the plate were derived basing on Hamilton’s principle. It states that the dynamics of a physical system is determined by a variation problem for a functional based on a single function, the Lagrangian, which contains all physical information concerning the system and the forces acting on it. In dynamic buckling problem the motion should be understand as the time dependent deflection.

The Hamilton’s principles for conservative systems states that the true evolution (compatible with constrains) of the system between two specific states in specific time range (t0, t1) is a stationary point (a point where the variation is zero) of the action functional Ψ. Action functional Ψ for i-th plate is described by following equation:


where Λ is the Lagrangian function for the system, Kis a kinetic energy of the system and Π is a total potential energy of the system.

The subscript idenoting i-th plate or strip in all equations in this subchapter is omitted – all equations are presented for one plate, which could be i-th plate, wall or strip of considered plate, beam-columns or girder (Figure 1).

Taking the action functional Ψ in form (9) the Hamilton’s principle can be written as:


The total potential energy variation δΠ for i-th thin plate (or strip) can be written in form:


where δQ is a variation of internal elastic strain energy:


and Ω is the volume of the plate and Sis its area, the volume can be expressed as lbhor Sh.

The variation of internal elastic strain energy for i-th plate or strip could be expressed by strain and sectional forces and moments in a following way:


The work Wof external forces (neglecting the out-of plane load) done on i-th plate can be written as follows:


where: p0(x), p0(y), τ0xy(x), τ0xy(y) are the prebuckling load applied to the middle surface of the considered plate (wall or strip)

For thin plates, it is assumed that the displacements uand vdo not depend on rotation w,xand w,yand therefore do not depend on the coordinate z. This approach results in exclusion of rotational inertia [50] in the equation for kinetic energy, which for the i-th thin plate (strip) can be written as:


The Hamilton’s principle, it is the variation of the action functional δΨ (10) for i-th thin plate (strip or wall) which after taking into consideration equations from (11) to (15) can be written as:


The Lagrangian function for the whole system is equal to the sum of the Lagrangian functions of all nplates of which the system was composed. To determine the variation of action δΨ for i-th plate, the following identity:

X δY=δ(XY)YδXE20

was used.

In the obtained equation, terms with the same variations were grouped, and then each of the obtained groups of terms (due to the mutual independence of variations) were equated to zero, giving:

  • equilibrium equations:

  • boundary conditions for lateral edges of the plate (x= const):

  • boundary conditions for longitudinal edges of the plate (y= const):

  • boundary condition for the plate corners (x= const and y= const):

  • initial conditions for t= const:


Above conditions are fulfilled for the entire structure, so if one apply the restrictions in moment of the initial t0 and in moment of the final t1 that the displacement variations are zero at all points of the structure. Then the system of equations (23) vanishes.

  • already used the relationship between deformations and internal forces and moments (8) or (9):


3. Solution method

To determine the critical loads, natural frequencies and the coefficients of the equation describing the postbuckling equilibrium path, the analytical-numerical method has been employed. The proposed method also allows analysing dynamic response of the structure subjected to pulse loading. Taking the time courses of deflections and applying the relevant dynamic buckling criteria it is possible to determine the dynamic critical load.

3.1. Equilibrium equations

The differential equations of equilibrium for orthotropic plate or strip directly from the equations (19) can be derived and have the form:


Above equilibrium equations after omitting the inertia forceshρu.., hρv..and hρw..becomes the equilibrium equations for thin plates allowing analysis of both local and global buckling mode.

3.2. Boundary and initial condition

As the wave propagation effects have been neglected, the boundary conditions referring to the simply supported columns at their both ends, i.e. x= 0 and x= l, according to (20), are assumed to be:


The condition written as a first of equations of (27) is satisfied for the prebuckling state and first-order approximation, the condition for deflection v(27) is satisfied for the first and second order of approximations, while the other two conditions are met for prebuckling state as well as for the first and second order of approximation. The condition of displacement in the ydirection in the prebuckling state can be found for example in [51]. This approach allows to take into account the Poisson effect on the edges of the walls of the column. The boundary conditions described by equations (27) assume the lack of displacement possibility of points lying at the loaded edges in the transverse vand normal wdirections to the surface in a wall or column. Furthermore, it is assumed that the moments Mix (as a vector parallel to the edge of the plate or end edge of the column walls) are zero.

For structures with material properties varying widthwise the strip model was adopted what forces the boundary conditions modification in the second order approximations [51]. Modification consists of changing the first condition of (27) onto the following form:


Summation is performed only for the Jnumber of the strips, between which the angle φi,i+1 (Figure 3) is equal to zero.

To determine the boundary conditions on the longitudinal edges of plates or free edges of columns with open cross-sections the equations (20) were used. Whereas, directly from equations (23) result the following initial conditions:

ui.(xi,yi,t=t0)=ui~(xi,yi)   and   ui(xi,yi,t=t0)=ui¯(xi,yi),vi.(xi,yi,t=t0)=vi~(xi,yi)   and   vi(xi,yi,t=t0)=vi¯(xi,yi),wi.(xi,yi,t=t0)=wi~(xi,yi)   and   wi(xi,yi,t=t0)=wi¯(xi,yi),E31

where the following functions ui¯, vi¯, wi¯, ui~ ,vi~wi~are given for the initial moment t= t0.

3.3. Interaction condition between adjacent plates

Static and kinematic junction conditions on the longitudinal edges of adjacent plates (Figure 3), according to (21), can be written as:




Figure 3.

The geometrical dimensions and local coordinate systems adjacent plates

3.4. Buckling and postbuckling equilibrium paths

A non-linear stability problem has been solved by means of the Koiter’s asymptotic theory. The displacement fieldU¯, and sectional force field N¯have been expanded into the power series with respect to the parameter ξ, - the buckling linear eigenvector amplitude (normalised with the equality condition between the maximum deflection and the thickness of the first plate h1).


It was assumed that the dimensionless amplitude of the initial deflections (imperfections) correspond to the considered buckling mode (for s-th buckling mode) is:


By substituting expansions (32) into equations of equilibrium (26) with neglected inertia terms (static buckling problem), junction conditions (30) and boundary conditions (27), the boundary problem of the zero (superscript (0) in Equations (32) and further), first (superscript (i)) and second (superscript (ij)) order has been obtained [18, 50, 52, 53]. The zero approximation describes the prebuckling state, whereas the first order approximation allows for determination of critical loads and the buckling modes corresponding to them, taking into account minimisation with respect to the number of half-waves min the lengthwise direction. The second order approximation is reduced to a linear system of differential heterogeneous equations, which right-hand sides depend on the force field and the first order displacements only.

The most important advantage of this method is that it enables us to describe a complete range of behaviour of thin-walled structures from all global (i.e. flexural, flexural–torsional, lateral, distortional buckling and their combinations) to the local dynamic stability. In the solution obtained, the shear lag phenomenon, the effect of cross-sectional distortions and also the interaction between all the walls of structures are included.

Having found the solutions to the first and second order of the boundary problem, the coefficients aijs, bijks have been determined [18, 50, 52, 53]:


where: λs – is the critical load corresponding to the s-th mode, L11 is the bilinear operator, L2 is the quadratic operator and σ(i), σ(ij) are the stress field tensors in the first and second order.

The postbuckling static equilibrium paths for coupled buckling can be described by the equation:

(1λλs)ξs+aijsξiξj+bijksξiξjξk=λλsξs*;(= 1,,N),E37

which for the uncoupled problem have the form:


where λcr is the critical load value.

In a special case, i.e. for the so-called ideal structure without initial imperfections (ξ*=0) and when the equilibrium path (a111) is symmetrical, the postbuckling equilibrium path is defined by the equation:


3.5. Natural frequencies

Determination of the natural frequencies is similar to the determination of critical buckling load and the natural frequencies are found by solving the eigenvalue problem.

Natural frequencies of thin-walled structures were determined by solving a dynamic problem, which uses the approach proposed by Koiter in his asymptotic stability theory of conservative systems in the first-order approximation [52].

To determine the natural frequencies [55] of the structure the adopted equilibrium equations (26) contain cross-sectional inertia forces acting in the direction normal to the middle surface of the plate (column wall) and in the middle plane of plate (i.e. hρu..≠ 0 and hρv..≠ 0).

3.6. Lagrange equations

In the dynamic analysis (while finding the frequency of natural vibrations [55]), the independent non-dimensional displacement ξ and the load factor λ become a function dependent on time, and dynamic terms were added to equations describing postbuckling equilibrium path. Neglecting the forces associated with the inertia terms of prebuckling state and the second-order approximations, and taking into account the orthogonality conditions for the displacement field in the first U¯(i)and second-order approximationU¯(ij), the Lagrange equations can be written as [56]:

1ωs2ξ..s+(1λλs)ξs+aijsξiξj+bijksξiξjξk=ξsλλs;(s=1,2, ,N)E40

where ωs is a natural frequency with mode corresponding to buckling mode; aijsand bijksare the coefficients (34) describing the postbuckling behaviour of the structure (independent of time); however the parameters of load λ and the displacement ξ are the functions of time t.

For the uncoupled buckling, i.e. the single-mode buckling (where index s= N= 1), the equations of motion may be written in the form:


It is assumed that in the initial moment of time t= 0 the non-dimensional displacement ξ, as well as the velocity of displacement are equal to zero, i.e.:

 ξ(t=0)=0 and ξ.(t=0)=0E42

The Runge-Kutta method [57] for solving the equation (39) requires the following substitutions:


which lead to the system of two differential equations. ‘‘Complete’’ equations of motion (41) are solved with the numerical Runge–Kutta method of order 8 (5,3), thanks to Dormand and Price (with step-size control and density output).


4. Exemplary results of calculations

The exemplary results of numerical calculation are presented in this sub-chapter. All results are obtained using explained above proposed analytical-numerical method (ANM) based on the nonlinear orthotropic plate theory.

The material properties (E – Young modulus, ν – Poisson ratio, G=E/[2(1+ν)] – Kirchhoff modulus; ρ – density) for materials taken into account are presented in Table 1.

material type:E
epoxy resin3.50.331249
glass fibre710.222450

Table 1.

Assumed material properties

The fibre composite material was modelled as orthotropic but for components (resin and fibre) the isotropic material properties (Table 1) was assumed. Necessary equations for material properties homogenization based on theory of mixture [57, 58] are as follows:


where Em and Ef are the Young’s modulus of elasticity for matrix and fibre, respectively, Gm and Gf are the shear modulus for matrix (subscript m) and fibre (subscript f), νm and νf are the Poisson’s ratios for matrix and fibre and f= Vf /(Vm + Vf) is the fibre volume fraction.

For static buckling the critical buckling load and corresponding modes are presented as well as the postbuckling equilibrium paths.

For dynamic buckling the proposed by Budiansky and Hutchinson parameter called Dynamic Load Factor DLF is introduced. The DLF is defined as a ratio of pulse loading amplitude to static buckling load. The results are presented of nondimensional deflection ξ versus DLF. The critical dynamic load factor DLFcr corresponding to dynamic buckling has been estimated using different criteria – the obtained results were compared.

For the proposed method the validation of the results was made by comparison with the other Authors [34] calculations (Figure 4) or with the results obtained with FEM [38]. The results presented in Figure 4 were obtained for thin (ratio length to thickness equals 200) aluminium square plate simply supported at all edges and subjected to sinusoidal pulse load. The time of pulse duration was equal to the period of natural vibration of the plate. The considered plate has a geometrical imperfection corresponding to buckling mode with amplitude equal to 0.05 of the plate thickness.

Figure 4.

The results of different calculation comparison

4.1. Plates

The rectangular thin plates simply supported on loaded edges with different boundary conditions along the unloaded ones were considered (Figure 5). On the longitudinal edges five different boundary condition cases were taken into account. Following notations is used in Figure 5: s – simply supported edge, c – clamped edge, e – free edge.

Figure 5.

Analysed plates with different boundary conditions

material:boundary conditionPcr[kN]ω [rad/s]
cc15.6 (m= 1)
13.9 (m= 2)
4423 (m= 1)
8363 (m= 2)
f= 0.5

Table 2.

Critical load Pcrand natural frequencies ω for analysed plates

Exemplary results were calculated for steel and epoxy glass composite (fibre volume factor f= 0.5) square plates subjected to rectangular compressive pulse loading. The buckling load for plate under analysis is presented in Table 2. The pulse duration Tp was equal to the period of natural vibration with mode corresponding to the buckling mode.

The dimensions of analysed plates were assumed as follows: the length (width) a= b= 100 mm and thickness h= 1 mm.

The geometrical imperfection was assumed in the shape corresponding to the buckling mode with amplitude ξ* = 0.01, where ξ* is an amplitude of deflection divided by the plate thickness. The Figure 6 presents postbuckling equilibrium paths for ideal flat composite plate (Figure 6a) and for plate with geometrical imperfection with amplitude ξ* =0.01 (Figure 6b).

Figure 6.

Postbuckling equilibrium paths for square ideal plates (a) and plates with imperfection (b) with different boundary conditions on non-loaded edges

In the dynamic buckling case the results are shown as graphs presenting nondimensional deflection ξ or radius rcalculated from real and imaginary part of maximal characteristic root of Jacoby matrix as a function of dynamic load factor DLF. The graphs mentioned above allow to find critical amplitude of pulse loading using the proposed criterion (PC) [35] and to compare the obtained results with Budiansky-Hutchinson (B-H) or Volmir (V) criteria. In brackets the notation used in Figures and Tables is given. The critical deflection according to Volmir criterion was assumed as ξcr= 1.

Figure 7.

Nondimensional deflection ξ (a) and maximum radiusrmax (b) vs. DLF for square plates with different boundary conditions on non-loaded edges [35]

Basing on curves presented in Figure 7 the critical value of dynamic load factor can be found. The comparison of obtained critical DLF values using different criteria is presented in Table 3. All critical DLF values except the case denoted as seobtained from the proposed criterion (PC) are in the range obtained from Budiansky-Hutchinson criterion (B-H). For the case denoted as se the greatest differences between critical dynamic load factors from the proposed and Budiansky-Hutchinson criterion were obtained but these differences are less than 10%.

boundary conditions for unloaded edgesmodeB-HVPC

Table 3.

Comparison of DLFcr obtained from different criteria for compressed plate

Figure 8.

Nondimensional deflection ξ vs. DLF for simply supported square plate

The comparison for postbuckling behaviour of the rectangular plate simply supported at all edges subjected to static and dynamic load are obtained using proposed analytical-numerical method and presented in Figure 8.

4.2. Segments of the girders

As a next example the static and dynamic buckling of composite (epoxy glass composite with different volume fibre fraction f) girders with open cross-section (Figure 9) is presented. The assumed boundary conditions on loaded edges correspond to simply support. The calculation was carried out for short segment of girder with length to web width ratio l/b1= 1 and for the following dimensions of the cross-section: b1/h= 50, b2/h= 25 and b3/h= 12.5.

Figure 9.

Cross-sections of analysed segment of the girders

The geometrical imperfection was assumed in the shape corresponding to the buckling mode with amplitude ξ* = 0.01. The static buckling load and fundamental flexural natural frequency obtained with analytical-numerical method for girders made of composite with different fibre fraction are presented in Tables 4. The static critical buckling loads are presented in Table 4 and postbuckling equilibrium paths are presented in Figure 10.

volume fibre fraction f:
critical load
natural frequency
ω [rad/s]
channel (Figure 9a)1526228110761076
channel with inner stiffeners (Figure 9c)2821421713081308
omega (Figure 9b)2819421413081308

Table 4.

Critical load and natural frequencies for analysed girder’s segments

Figure 10.

Postbuckling equilibrium paths for segment of girders

Buckling load and natural frequency for girder segment with omega and stiffened channel cross-section are similar – it is true only for local buckling case. For girder with channel cross-section the buckling was caused by flanges – this is a reason why for this cross-section the buckling load and natural frequency are smaller than for two others analysed cross-sections. Looking at obtained results (Table 4) it can be seen that increasing the volume fibre fraction fleads to an increasing the buckling loads as well as the natural frequencies. The postbuckling equilibrium paths for stiffened cross-section (channel with inner stiffeners and omega) overlap. The postbuckling path for channel cross-section lies below the equilibrium paths of girders with stiffened cross-section – it is obvious because the girders with stiffened cross-section have similar stiffness and the girder with channel cross-section is more flexible.

In dynamic buckling case the time of pulse duration Tpwas assumed as the period of natural fundamental flexural vibration corresponding to the local buckling mode. Considered shapes of pulse loading are presented in Figure 11.

Figure 11.

Shapes of pulse loading: three triangular impulses T1 (a), T2 (b), T3 (c), rectangular R (d) and sinusoidal S (e)

Figure 12.

Dimensionless deflection vs.DLFfor channel girder – results comparison obtained with analytical-numerical method ANM and finite element method FEM [56]

The results presented in Figure 12 were obtained with the proposed analytical-numerical method (ANM) and compared with FEM computations. They are similar for both assumed shapes of pulse loading (Figure 11): rectangular (R) and sinusoidal (S). Nevertheless, for rectangular shape pulse loading some small differences in deflection are visible for DLFgreater than 2. It should be noted that obtained curves (Figure 12) from both methods allow to find the same critical dynamic load factor DLFcrusing Budiansky-Hutchinson or Volmir criterion – for rectangular pulse loading DLFcr≈1.4 (both criteria) and for sinusoidal pulse loading DLFcr≈2.1 (Budiansky-Hutchinson) or 2.3 (Volmir).

In Figure 13 the dynamic response comparison of girders made of composite material (f= 0.5) with different cross-section subjected to triangular T3 pulse was presented. The curves for omega cross-section and channel cross-section with inner stiffeners cover each other’s.

Figure 13.

Nondimensional deflection ξ as a DLF function for girder with different channel cross-section subjected to T3 pulse loading [56]

Dynamic responses for girder with channel cross-section with inner stiffeners for different pulse loading are presented in Figure 14. The curves denoted by S=Rwere obtained for sinusoidal pulse loading with the same area as rectangular pulse loading (for the same pulse duration the amplitude was higher for sinusoidal pulse). The highest deflection was obtained for pulse denoted by S=Rbecause this pulse has the highest amplitude. The rest of compared pulses have the same amplitude and the same duration. Analysing the curves ξ(DLF) for rectangular, sinusoidal and three triangular impulses it can be seen that the highest increment of deflection for the smallest DLFtakes place for rectangular pulse loading.

Figure 14.

Dimensionless deflection vs. dynamic load factor for different shapes of pulse loading - channel with inner stiffeners, composite materialf= 0.7 [56]

The comparison of DLFcrobtained using Budiansky-Hutchinson criterion for girders with different cross-section made of composite materials (f= 0.5) are presented in Table 5. Only the average values from the obtained critical ranges are presented. The dynamic load factors for different cross-sections are in the same relation as buckling loads (Table 4) – the same or similar DLFcrfor cross-sections with stiffeners (omega and channel with inner stiffeners).

4.3. Columns

Next, the exemplary results for the dynamic interactive buckling of channel cross-section the columns are presented. The columns subjected to rectangular compressive pulse loading were analysed. The calculation was carried out for columns with various length to web width ratio l/ b1 = 4; 6 and 8 and for the following dimensions of the cross-section: width of the web to its thickness b1/h= 100, width of the flange to its thickness b2/h= 50.

analysed cross-section:

type of pulse
channelomegachannel with inner stiffeners

Table 5.

The DLFcrvalue for different shape of applied pulses

Figure 15.

Nondimensional deflection ξ (a) and maximum radiusrmax (b) as a function of DLF for channel beam-columns with different length ratiol/b1 and pulse durationTp[35]

The interaction between global buckling mode m= 1 and the first local buckling mode m> 1 was considered. The geometrical imperfection were assumed in the shape corresponding to the buckling mode with amplitude ξ* ≡ ξ2*equals 1/100 wall thickness for local mode and ξ1* equal to length to one thousand wall thickness (l/1000h) for global mode. Time of pulse duration Tpwas assumed as T1equal to the period of natural fundamental flexural vibration or Tm(where mis a number of half waves of local buckling mode) equal to the period of natural vibration with mode corresponding to the local buckling mode.

Figure 15 presents dimensionless deflection ξ as a function of dynamic load factor and maximal radius rmax calculated for maximal characteristic root of Jacoby matrix as a function of DLF.

From curves presented in Figure 15a the critical value of dynamic load factor based on Volmir (V) or Budiansky-Hutchinson (B-H) criterion can be found. The curves presented in Figure 15b help to find critical DLF value based on proposed criterion (PC) [35]. The obtained critical DLF’s according to mentioned above criteria are presented in Table 6.

columns length ratiopulse duration
buckling modesB-HVPC
l/b1= 4T1= 1.6m=1; 3 local1.0÷
l/b1= 6T1= 1.9m=1; 5 local0.95÷
l/b1= 8T1= 2.7m=1; 6 local1.1÷
l/b1= 4T3= 0.8m=1; 3 local1.6÷1.751.401.58
l/b1= 6T5= 0.7m=1; 5 local1.6÷1.751.391.59
l/b1= 8T6= 0.8m=1; 6 local1.6÷2.051.431.51

Table 6.

Comparison of DLFcr obtained from different criteria for interactive buckling

The comparison of the obtained results shows that they are in good agreement. In all cases with pulse duration equal to the period of natural vibration of the form corresponding to local buckling mode the results obtained from the proposed criterion (PC) are between the results obtained from Volmir (V) and Budiansky-Hutchinson (B-H) criteria. For loading with time of pulse duration equal to the period of natural fundamental vibration T1 the critical DLF values obtained using the proposed criterion (PC) are equal or a bit greater (about 6%) than the critical dynamic load factors from Budiansky-Hutchinson (B-H) or Volmir (V) criteria.

Should be pointed out that in the dynamic buckling problem also for the short columns the multimodal buckling analysis should be carried out. It has been proven on exemplary channel columns with following dimensions: b1/h= 100, b2/h= 50, b3/h= 25 and l/b1= 4.

The problem has been calculated with the analytical-numerical method and the finite element method [56].

The dimensionless deflection ξ as a function of dimensionless time (time divided by pulse duration) for channel column is presented in Figure 16. The characteristic points are located in the middle cross-section of the columns and in the middle of the web (point 1 – Figure 16) and on the edge between the web and the flange (point 2 – Figure 16).

Figure 16.

Deflection in time for channel columns subjected to rectangular pulse loading with the DLF = 1.6

Analysing the results presented in Figure 16, it can be said that the global buckling appears for channel cross-section columns. The column edge deflections are greater than deflections of the middle part of the web. The FEM results of calculations presented in Figure 16 have initiated the need for a multimodal buckling analysis also for short columns subjected to pulse loading.

Results for linear buckling and modal analyses obtained with proposed analytical-numerical method are presented in Table 7. As it will be presented below (see Figure 17) the four modes should take into consideration in ANM to obtain similar results to this obtained with FEM. The finite element method gives results (global mode) even in the case when only one buckling mode as the initial imperfection (for example, the local buckling mode m=3) has been taken into account [56].

local modem= 353614
primary local modem= 1123312
secondary local modem= 1’972880
global modem= 1’’51222001

Table 7.

Buckling stress and natural vibration for the channel column

A comparison between the results obtained with the analytical-numerical method and the finite element method on plots presenting a dimensionless deflection vs. a dynamic load factor for columns with channel cross-sections are shown in Figure 17. Some differences appear because the analytical-numerical model has only a few degrees of freedom in contrary to FE model, which has thousands DoF. However the curves presented in Figure 17 are different the critical DLF values estimated using the Budiasky-Hutchinson criterion is similar. From the ANM, the DLFcr= 2.7 and from the FEM, the DLFcr= 2.6.

Good agreement between the results obtained with ANM and FEM is possible because the interactive dynamic buckling problem has been solved in the analytical-numerical method. Four modes have been taken into account. The buckling stress and the natural frequency obtained with the analytical-numerical method for all the modes assumed in the multimodal analysis are listed in Table 7. The buckling modes taken into consideration in the multimodal analysis are presented in Figure 18, correspondingly.

Figure 17.

Dimensionless edge deflection ξ vs. theDLFfor channel columns subjected to rectangular pulse loadingTp=T3= 1.6 ms [56]

Figure 18.

Buckling modes for the channel cross-section column


5. Conclusion

Taken into consideration the nonlinear thin plate theory for orthotropic material allows, as is presented in exemplary results of calculation, to analyse thin-walled structures composed of flat plates and subjected to static and dynamic load. The nonlinear orthotropic plate theory is the base for the proposed analytical-numerical method which allows to find buckling load with corresponding buckling mode, natural frequencies with corresponding modes and to analyse the postbuckling behaviour – drawing the postbuckling equilibrium paths for plate, segment of girders or columns made of isotropic, orthotropic or even composite materials. As it was shown not only static load can be considered but also dynamic load with intermediate velocity – the dynamic buckling can be analysed using assumed plate theory and the proposed method of solution.

The proposed analytical–numerical method gives almost the same results for eigenvalue problem (buckling loads, natural frequencies with corresponding modes) and similar results for dynamic buckling as the finite element method. However the dimensionless deflection versus dynamic load factor relation obtained with both (proposed and FEM) methods are not identical (especially for higher DLF value). These relations allow to find similar critical value of DLF taking into consideration one of the well-known criterion. The differences in the dimensionless deflection ξ appear because the numerical model in the FEM has more degrees of freedom than the model in the analytical–numerical method, but the results from the ANM are obtained in a significantly faster way than those from the finite element method.


This publication is a result of the research work carried out within the project subsidized over the years 2009-2012 from the state funds designated for scientific research (MNiSW - N N501 113636).


  1. 1. BernoulliJ.EulerL.1910Abhandlungen uber das Gleichegewicht und die Schwingungen der Ebenen Elastischen Kurven, Wilhelm Engelmann, Nr 175, Leipzig.
  2. 2. TimoshenkoS. P.GereJ. M.1961Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadMcGraw-Hill Book Company, Inc. NewYork, Toronto, London.
  3. 3. Volmir A.S.1967Stability of Deformation Systems, Moscow, Nauka, Fizmatlit /in Russian/.
  4. 4. AliM. A.SridharanS.1988Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadInternational Journal of Solids and Structures245481486
  5. 5. BenitoR.SridharanS.1985Mode Interaction in Thin-Walled Structural Members, Journal of Structural Mechanics,124517542
  6. 6. ByskovE.1988Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadMechanics and Structures Machines,154413435
  7. 7. KoiterW. T.PignataroM.1974An Alternative Approach to the Interaction Between Local and Overall Buckling in Stiffened Panels, Buckling of structures Proceedings of IUTAM Symposium, Cambridge:133148
  8. 8. KolakowskiZ.1989Some Thoughts on Mode Interaction in Thin-Walled Columns Under Uniform Compression, Thin Wall Structures,712335
  9. 9. KolakowskiZ.1993Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadWall Structures,153159183
  10. 10. KolakowskiZ.1993Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadInternational Journal of Solids and Structures301925972609
  11. 11. Manevich A.I.1988Interactive Buckling of Stiffened Plate under Compression. Mekhanika Tverdogo Tela,5152159in Russian/.
  12. 12. MoellmannH.GoltermannP.1989Interactive Buckling in Thin-Walled Beams, Part I: Theory, Part II: Applications. Internatioan Journal of Solids and Structures, 25(7):715728and 729-749.
  13. 13. PignataroM.LuongoA.RizziN.1985Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadThin Wall Structures, 3(4):283321
  14. 14. PignataroM.LuongoA.1987Asymmetric Interactive Buckling of Thinwalled Columns with Initial Imperfection, Thin Wall Structures,55365386
  15. 15. PignataroM.LuongoA.1987Multiple Interactive Buckling of Thin-Walled Members in Compression, Proceedings of the International Colloqium on Stability of Plate and Shell Structures, Ghent: University Ghent:235240
  16. 16. SridharanS.AliM. A.1985Interactive Buckling in Thin-Walled Beam-Columns, Journal of Engineering Mechanics ASCE,1111214701486
  17. 17. SridharanS.AliM. A.1986Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadInternational Journal of Solids Structures,224429443
  18. 18. KolakowskiZ.Kowal-MichalskaK. .Eds1999Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic Loadity of Lodz Press, A Series of Monographes, Lodz, 1999.
  19. 19. KolakowskiZ.KubiakT.2007Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadComposite Structures222232
  20. 20. Ari-GurJ.SimonettaS. R.1997Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadComposites Part B28301308
  21. 21. CuiS.HaiH.CheongH. K.2001Numerical Analysis of Dynamicbuckling of Rectangular Plates Subjected to Intermediate Velocity Impact, International Journal of Impact Engineering,252147167
  22. 22. AwrejcewiczJ.KryskoV. A.(2003) Nonclassical Thermoelastic Problems in Nonlinear Dynamics of Shells, Springer-Verlag, Berlin.
  23. 23. AwrejcewiczJ.AndrianovI. V.ManevitchL. I.2004Asymtotical Mechanics of Thin-Walled Structures. A Handbook, Springer-Verlag, Berlin.
  24. 24. KounadisA. N.GantesC.SimitsesG.1997Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadInternational Journal of Impact Engineering1916380
  25. 25. Mania R.J.2011Viscoplastic Thin-Walled Columns Response to Pulse Load, Proc. ICTWS, Thin-Walled Structures, (ed.) Dubina D.:415430
  26. 26. BudianskyB.1965Dynamic Buckling of Elastic Structures: Criteria and Estimates, Report SM-7, NASA CR-66072.
  27. 27. Simitses G.J.1990Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadSpringer Verlag, New York.
  28. 28. VolmirS. A.1972Nonlinear Dynamics of Plates and Shells, Science, Moscow /in Russian/.
  29. 29. BudianskyB.HutchinsonJ. W.1966Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadProceedings of the Eleventh International Congress of Applied Mechanics, Goetler H. (ed.), Munich:636651
  30. 30. HutchinsonJ. W.BudianskyB.1966Dynamic Buckling Estimates, AIAA Journal,43525530
  31. 31. BisagniC.2005Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadThin-Walled Structures433499514
  32. 32. GilatR.AboudiJ.2002The Lyapunov Exponents as a Quantitive Criterion for the Dynamic Buckling of Composite Plates, International Journal of Solid and Structures392467481
  33. 33. KubiakT.2005Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadInternational Journal of Solid and Structures422055555567
  34. 34. PetryD.FahlbuschG.2000Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadThin-Walled Structures38267283
  35. 35. KubiakT.2007Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadThin-Walled Structures888 EOF892 EOF
  36. 36. TeterA.2010Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadMechanics and Mechanical Engineering142165176
  37. 37. TeterA.2011Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadThin-Walled Structures495589595
  38. 38. Kowal-MichalskaK.(ed2007Dynamic Stability of Composite Plate Structures, WNT, Warsaw.
  39. 39. SimitsesG. J.1987Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadApplied Mechanics Revision,401014031408
  40. 40. GrybosR.1980Statecznosc Konstrukcji pod Obciazeniem Uderzeniowym, PWN /in Polish/.
  41. 41. Mania R.J.2011Dynamic Buckling of Orthotropic Viscoplastic Column, Thin-Walled Structures,495591588
  42. 42. ManiaR.Kowal-MichalskaK.2006Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadThin-Walled Structures
  43. 43. ManiaR. J.Kowal-MichalskaK.2010Elasto-Plastic Dynamic Response of Thin-Walled Columns Subjected to Pulse Compression, Proc. of SSTA 2009, CRC Press:183186
  44. 44. KrolakM.ManiaJ. R.(eds.(2011) Stability of Thin-Walled Plate Structures,1of Statics, Dynamics and Stability of Structures, A Series of Monographs, TUL Press, Lodz.
  45. 45. CzechowskiL.2008Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadMechanics and Mechanical Engineering124309321
  46. 46. Kowal-MichalskaK.2010Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadMechanics and Mechanical Engineering14269279
  47. 47. Kowal-MichalskaK.ManiaJ. R.2008Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadMechanics and Mechanical Engineering122135146
  48. 48. KubiakT.2007Dynamic Coupling Buckling of Thin-Walled Columns, 998 Scientific Bulletin of TUL, Lodz /in Polish/.
  49. 49. PietraszkiewiczW.1989Geometrically Nonlinear Theories of Thin Elastic Shells. Advances in Mechanics,12151130
  50. 50. WozniakCz.(ed.2001Mechanika Sprezysta Plyt i Powlok. Wydawnictwo Naukowe PWN, Warszawa /in Polish/.
  51. 51. KubiakT.2001Postbuckling Behavior of Thin-Walled Girders with Orthotropy Varying Widthwise, International Journal of Solid and Structures, 38 (28-29): 4839-4856.
  52. 52. Koiter W.T.1963Elastic Stability and Post-Buckling Behaviour, Proceedings of the Symposium on Nonlinear Problems, Univ. of Wisconsin Press, Wisconsic:257275
  53. 53. KrolakM.(ed.1990Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic Loadin Polish/.
  54. 54. KrolakM.(ed.1990Statecznosc, Stany Zakrytyczne i Nosnosc Cienkosciennych Konstrukcji o Ortotropowych Scianach Plaskich, Monographs, TUL Press, Lodz, /in Polish/.
  55. 55. TeterA.KolakowskiZ.2003Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadThin-Walled Structures414291316
  56. 56. KubiakT.2011Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadComputers and Structures2001 EOF2009 EOF
  57. 57. PrinceP. J.DormandJ. R.1981Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadJournal of Computational and Applied Mathematics716775
  58. 58. JonesR. M.1975Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadInternational Student Edition, McGraw-Hill Kogakusha, Ltd., Tokyo.
  59. 59. KellyA. .Ed.1989Nonlinear Plate Theory for Postbuckling Behaviour of Thin-Walled Structures Under Static and Dynamic LoadPergamon Press.

Written By

Tomasz Kubiak

Submitted: March 11th, 2012 Published: October 24th, 2012