Large Deformations of Framed Structures Under Static and Dynamic Loads?

Download as pdf or txt
Download as pdf or txt
You are on page 1of 9

Compders & Sln!crums. Vol. 6. pp. 53947. Peqamon Press 1976.

hinted in Great Britain

LARGE DEFORMATIONS OF FRAMED STRUCTURES


UNDER STATIC AND DYNAMIC LOADS?
CENAP ORAN+ and ASLAM KASSIMALIO
Civil Engineering Department, University of Missouri-Columbia,Columbia, MO 65201,U.S.A.

(Received 30 Nouember 1975)

Abstract-With reference to the problem of large deformations and stability of elastic framed structures, this paper
explores the computational capabilities of a general beam-column type method which was recently developed by the
senior author. The method is flexible in that the coordinate system used may be either Eulerian or Lagmngian. In
addition, various types and levels of consistent approximations can be introduced into the analysis in a rather routine
fashion.
In an effort to evaluate the merits of the method in the static case, extensive numerical studies were carried out on
a groupof specially selected and relatively simple st~ctural systems. In formulating the dynamic case, use was made
of a well-known numerical integration technique, namely, the so-called “Newmark’s /3-method”. Again, numerical
studies were carried out, although on a smaller scale than in the static case. These studies clearly indicate that the
suggested method is both practical and highly accurate.

INTRODUCIION equilibrium equations of the structure in terms of joint


The methods that have been suggested in the literature for displacements, {F} must be expressed in terms of {u} for
the analysis of large deformations and stability of elastic each individual member. For this purpose, it is convenient
framed structures are either based on the finite element to introduce two alternate sets of member end forces and
concept, or else they make use of the so-called end displacements as shown in Fig. l(b) and (c). From
“beam-column’* approach. In the former group, the basic geometry,
idea is relatively simple in theory, but a large number of
elements are often needed for satisfactory accuracy. In {F) = [R]{P}; {Ag}= [RIT{Au] (1)
the latter group, member force-deformation relations are
in which
determined from a more complex, but also more refined,
analysis so that fewer elements are generally needed than
in the case of a finite element method.
The study presented herein deals mainly with methods
of the beam-column type. Some dilbculties associated
with these methods, including a number of inconsistencies with
encountered in the literature, were examined by Oran[ l- m =cosu; n = sin ff (3)
31 in the course of a theoretical investigation that
ultimately led to “consistent” tangent stiffness matrices and the superscript T denotes transpose. Note that the
for plane and space frames with prismatic and non- angle u refers to the deformed configuration of the
prismatic members. No attempt was then made, however, member.
to evaluate the merits of these theoretical developments Similarly,
from a practical computational viewpoint. Such an
exploratory study constitutes the main objective of the {I’, = [ri]{s}; {An} = [@‘{ha} (4)
present investigation. Also considered briefly herein is the
dynamic analysis of framed structures in the large in which
deformation range.
In what follows, the analysis of plane frames under .?,=M, St=& s,=QL (5)
static loads is first examined in detail. The ideas thus
developed are then extended to more complex problems c,=e, i&=9* a,=; (6)
(i.e. space frames and/or dynamic loads) without,
however, giving detailed derivations. and

STATlC ANALYSIS OFPLANEFRAMES - 0 0 1


Member end efects. Consider an arbitrary prismatic 1 1+so
1
member of a plane frame, and let {F} and {v} denote 1+6
member end forces and end displacements in global [B]= 1 0 0 (7)
coordinates as shown in Fig. l(a). In formulating the 0 0 -1
-- 1 1
1ts
tPresented at the Second National Symposium on Computer-
ized Structural Analysis and Design at the School of Engineering -- 1 0 1
O
and Applied Science, George Washington University, Washing- with
ton, D.C., 29-31 March 1976.
SProfessor.
OGraduate Student.
(8)

539
C. ORAN and A. KASSIMALI

in which xi = generalized coordinates (essentially, transla-


tions and rotations of the joints); P, = load terms; and
J = nonlinear functions representing internal forces cor-
responding to x, (obtained by combining member end
forces acting at the respective joints). When the problem
is analyzed by the beam-column approach, (f} cannot be
expressed explicitly in terms of {x}, because, as previ-
ously indicated, member axial forces do remain in the
picture as auxiliary parameters.
(a) Member forces in global coordinates Equations (12) can be rewritten in terms of differentials
as

{AP} = [T]{AX} (13)

where {Ax} and {AP} represent displacement and load


increments respectively, and

[T] = [$I
I
(14)
(b) Member forces in member coordinates

is the “system tangent stiffness matrix”. As will become


apparent shortly, eqn (13) plays a central role in
practically all computational methods that have been
suggested for the solution of eqns (12).
Member tangent stiffness matrices. In determining [T],
the first step consists in establishing incremental relation-
ships between member end forces and end displacements
for each member of the structure. For this purpose, let
(c) Member forces and relative member deformations
IAF1= [WV) (15)
Fig. I. Member end effects.
in which TT1=member taneent stiffness matrix. As
Note that L( 1+ 8) represents the chord length in the shown in dk&il elsewhere [ 1, i],
deformed configuration of the member.
-What is needed at this stage is a relationship between V’l= UU[~l[tl[~l’W +$, S,[RlW’I[W (16)
{S} and {zi}. Within the limitations of the beam-column s
theory [ I],
in which [t] is defined through the equation
{A.?}= [t]{Ati} (17)
and is given by
M2 = ~(~~0, t c,e,)

(18)

in which A = area of cross section; I = moment of inertia:


E = modulus of elasticity; c,, c2 = stability functions; and
with
cb= b,(e,
t e,)*tbde,- B?)* (IO)

where b,, b2 = bowing functions. GI=cle,tc;e2


The chain of relationships defined in the preceeding G2 = de, t c;e2
paragraphs does implicitly express {F} in terms of {u}. (19)
n*
However, the situation is complicated by the fact that c,,
ct, b, and b2 are functions of a dimensionless axial force
parameter,
In these equations,

4=$g (11)
(20)
which, in turn, is an implicit function of {u}.
Equilibrium equations. The equilibrium equations of
the structure as a whole can be written symbolically as is the usual “slenderness ratio” for a column pinned at
both ends, and a prime superscript on c, or bi denotes a
fihx*,..., x,)-P,=0 for i=l,2 ,..., n (12) differentiation with respect to q. Also,
Largedeformationsof framedstructures
understaticanddynamicloads 541

(21) [T]=LBl~~l~~l~+~,S[~“‘l (31)


_

(22) in which
WI = [RIWI, t81 (32)
in which f’) and p” are 3 x 3 matrices defined by
[G’“‘]= [B][R],[g’“‘](R],‘[~Ir (33)
ff’ = 0 except fli = J+$$i
=h (23) can be evaluated directly, without actually performing
any matrix multiplication on the computer. Note that, for
-1 consistency, the assumption of smah p would have to also
A’=0 except fg =- (24)
(1+@ be incorporated into the system equilibrium equations,
eqn (12), via eqn (1).
Equation (16), while it is somewhat different in form, is The approximate theory outlined in this section may be
essentially equivalent to eqn (45) of Ref. [ 11,except that S interpreted as a transition from Eulerian to Lagrangian
was considered to be small in comparison with unity in coordinates, although the relationships that have been
eqn (45). derived may or may not be consistent with those given by
Eulerion us Lugmngian coordinates. The relationships other writers who use Lagrangian coordinates.
derived in the preceeding sections are valid for arbitrarily It should also be mentioned at this time that the
large member rotations provided, of course, that relative practical advantages of the Lagrangian approach are
member deformations are small enough to justify the use better appreciated in space than in plane problems.
of the beam-column theory. If the rotations are also small, Computational techniques and approximations. Two
then it would be possible to introduce a useful simplifica- techniques are commonly used in solving nonlinear
tion into the analysis, with only minor effects on the problems. The first consists in applying the loads in small
accuracy of the method. increments, and determining the corresponding changes in
For this purpose, let the configuration of the structure from a sequence of
linearized analyses (Fig. 2a). This technique does not
a =d+p (25) involve iteration, but requires that the tangent stiffness
matrix be updated at each load level. The second
in which ~5refers to the initial undeformed orientation of technique is an iterative one (Fig. 2b). For any desired
the member, and p represents the angle of rotation of the load level, an approximate solution is first assumed or
chord. The transformation matrix [RI can now be calculated; it is then improved step by step via a
expressed as Newton-Raphson type of iteration until the joint equilib-

WI = mR1, (26)

in which [I?] and [RI, are to be obtained from [R] in eqn


(2) by replacing a by C?and p, respectively.
If p is small, then

cos p =s1; sinp-p

1
(27)
so that, from eqn (2),

1
-p 0
bl= p 1 0 w
[ 0 0 1

Note that, in Lagrangian coordinates,

{I?}= [B]‘(u); {E} = [R]‘(F) (29) Lwd


AFisumd
sohltii
I -
Substituting eqn (26) into eqn (16),

ITI = ~~I~Rl,~~l~tl~~l’~Rl~~~lr

+$ S,[R][R],[~“)][R],T[~]~. (30)
k-1
Judging from appearances, eqn (30) would seem to be
even more complicated than eqn (16). If the actual
computational work is performed according to the
detailed sequence indicated in eqn (30), this observation
may indeed be valid. For computational purposes, (b) Ileratii
however, eqn (30) can be rewritten in alternate and more
compact forms, such as Fig. 2. Commoncomputationaltechniques.
542 C. ORAN and A. KASSEMALI

rium equations (eqn 121,are satisfied within a prescribed indicating that the imperfect system will not buckle (at
tolerance. Various combinations of iterative and finear- least in a theoretics sense).
ized incremental techniques are also commonly used, Numerical sohtions. In an effort to assess the compu-
In analyzing plane frames by the beam-column ap- tational merits of the proposed general method of
proach outlined in the preceding sections, one may analysis, six special problems were selected for detailed
consider the possibility of reducing the required computa- numerical studies. These problems, which are summar-
tional work by introducing some simplifying approxima- ized schematically in Fig. 4, are relatively small and
tions. For example, one could: simple, but have certain interesting and representative
(a) Use Lagrangian coordinates. features:
(b) Set Gi = Gz = H = 0 in eqn (18). Then, [t] would be (a) No instability. Exact solution available.
a function of q only. (b) Limit point when shallow, b~ur~ation point when
(c) Use undeformed member Iength in writing member steep. Exact solution available.
equilibrium equations, i.e. neglect S in comparison with fc) Column (perfect or imperfect). Exact solution
unity in eqns (7, 23, 24). available for perfect case.
These and other conceivable approximations can be (d)
(e) Previously examined by other investigators. Selected
used individually or in various combinations.
mainly for comparison purposes.
Detection ofinstability. It may be of interest to note at (f)
this point that many writers view a lack of convergence in
the iterative process as an indication of instability. This
reasoning, while it is clearly justified from a theoretical
viewpoint, appears to be unreliable as a practical buckling
criterion. Better results are obtained by noting[4] that the
determinant of the system tangent stiffness matrix
vanishes at all critical points,

det [T] = 0. (34)


I(a) Cantilevw bwm
(b) Twa+v truss

If the value of det [T] is plotted as function of a load


parameter, P, a limit point (i.e. snap-though bucking) is
characterized by a vertical tangent as illustrated in Fig.
3(a), and a bifurcation point by a finite angle of
intersection of the type shown in Fig. 3(b). In the presence
of a small imperfection, the bifurcation point may
I(cl Cantilever calunn :d)One-stow
bent

Lz
transform into a limit point indicating that the system is
imperfection-sensitive, or it may disappear completely
F-’ 9 /: l-=$=-j
e (

G-----j
(ff Fixed beam

Fig. 4. Planar structures selected for numerical studies.

Each problem was solved approximately 6-8 times, by


using a slightly different numerical technique in each case,
*~ y----- with particular reference to the following considerations:
Eulerian vs Lagrangian coordinates, linearized incremen-
tal vs iterative procedures, small vs large load increments,
effect of neglecting 6 and/or G,, Gz, H. Clearly, it would
(a) Typical Ikit PM be neither possible nor desirable to present extensive
numerical data within the space limitations of the present
report. The following brief study is accordingly limited to
a relatively small number of representative solutions.
The results summarized in Fig. 5 are for a cantilever
beam with a concentrated load acting at its tip. The solid
line represents the exact solution[S], in terms of elliptical
integrals. Also incuded in the figure are results obtained
by using some of the approximations that were suggested
in a previous section of this report. It is seen that the
iterative solutions (both Eulerian and Lagrangian) are in
excellent agreement with the exact solution. The linear-
ized incremental solutions do, as would be expected,
show a deviation toward the linear solution. It is
interesting to note that the error becomes more substan-
tial when GI, G2 and H are neglected.
(b) Typical axon paint The results obtained for a simple two-bar truss are
summarized in Fig. 6. Of the two iterative solutions
Fig. 3. Detection of instability. shown, the Eulerian does happen to coincide with the
Largedeformations of framed structures under static and dynamic loads 543

l hra?im (Et&ion)

0 IO 20 30 40 50 m
verbl~~ofpoint8,h

Fig. S. Load deflection curvesfor cantilever beam.

0 4 6 12 16 20 24 26
'&ricddeflcchand,h

Fig. 6. Load deflection curves for shallow truss.

exact solution, and the Lagrangian, while it is not exact, is converge at subtantially lower load levels. A lack of
seen to be highly accurate. As in the case of the cantilever convergence in the iterative process was viewed by those
beam of Fig. 5, the linearized incremental solutions do authors as an indication of instability. In light of the
appear to have a tendency to drift away from the exact additional numerical data presented herein, such a
solution. However, their accuracy seems to remain conclusion would not necessarily be justified. It should
satisfactory for most practical purposes. also be noted that the bent shown in Fig. 9 is not
As for the det [r] vs P relationship, it is seen, from Fig. imperfection-sensitive, so that, in the presence of a small
7, that all the solutions that are presented do exhibit the imperfection, the structure would not even buckle.
typical behavior that characterizes a limit point. The A recent comparative study by Ebner and Ucciferro [7]
results obtained for a steeper arch are summarized in Fig. gives numerical results obtained by several different finite
8, where the existence of a bifurcation is quite obvious. element techniques for a group of specially selected
Note also that, in the presence of a small imperfection, the problems. The group includes, in particular, the cantilever
critical point transforms into a limit point, thus confnming beam of Fig. 4(e), with L, = 52.30 in., .L2= 50.72 in.. E =
the imperfection-sensitive character of the system. 30x lVpsi, I= 1/6OOOin.‘, A = 1/5in.*, P,=0.85lb.,
The results obtained for a simple one-story bent are P2 = 1.35lb. The results reported for this relatively simple
shown in Fig. 9. It is particularly important to note that the problem are reproduced in Table 1, along with results
iterative procedures developed as part of the present obtained by the linearized incremental version of the
investigation did function properly (i.e. with no significant present method in Eulerian coordinates. It is seen that the
convergence difficulties) up to, even beyond, the critical present solution, although it uses only two elements, is at
load of the limiting case of the perfect bent. It may be of least as accurate as, if not even more so, than the
interest to note, in this connection, that an alternate solutions obtained by various versions of the finite ele-
iterative technique suggested by James et a!.[61 failed to ment approach.
C. ORAN and A. KASSIMALI

STATIC ANALYSIS OF SPACE FRAMES


Static analysis of space frames [2] is essentially similar
to that of plane frames except, of course, that all matrices
are now substantially larger. It may be recalled, in this
connection, that a joint has six degrees of freedom in
space (as opposed to three ip plane), and a member has
twelve end displacements in space (as opposed to six in
plane). From a theoretical viewpoint, the most important
difference between the two classes of problems stems
from the fact that large rotations cannot be treated as
vectors in space problems, thus implying that the vector
{u} is no longer defined, while {Au} still is. For book-
keeping purposes, the deformed orientation of each joint
is represented by a “joint orientation matrix” the columns
of which are the direction cosines of three mutually
0 Iteration
Eulerd
+ lkmtian (LlqmngM)
perpendicular lines that are assumed to be rigidly con-
nected to the joint. For a detailed examination of these
ideas, the reader is referred to Oran[2].
When the rotations are small, the analysis can be
simplified substantially by switching to Lagrangian coor-
dinates. By analogy to eqn (26), [R] would then be
expressed as a product, except that [RI,, would be a
function of not one but three rotation components p,, pr.
pr. More importantly, the Lag-rangian approach would
I make it possible to use the joint rotations as generalized
coordinates, thus eliminating the need for joint orientation
0 4ccl 800 1200 1603 3 matrices.
Lwd P, kips Numerical solutions were obtained for two simple
space problems, one of which is the cantilever beam
Fig. 7. Load-determinant curves for shallow truss. shown in Fig. 4(a) and previously treated as a plane
I , r I I # 1 -

3- 0 lherution (Eulerion)
+ Iteratiw (Eulerian with 8 *O)

Fig.8. Load-determinantcurvesfor steep truss

Table I. Comparison of solutions for the cantilever beam shown in Fig. 4(e)
Largedeformationsof framedstructuresunderstaticand dynamic toads 545

t.lcecddeflsclwnofpdnr8,n

Fii. 9. haddefkction curvesfor one storybent.

problem. In an effort to check on the adequacy of known at time tK, and it is desired to determine the
book-keeping process suggested herein for rotations (i.e. configuration at fK+,= fK+ At. For this purpose, let eqn
the use of joint orientation matrices) solutions were (35) be rewritten in an incremental form as
obtained by con~de~g three different initial o~en~tions
in space for the fzintilever beam. AII three solutions thus MM + t~l&] - tbp] = 0 (37)
obtained were found to agree with the results of the plane
case withii the tolerance specified for the iterative pro- in which {AP} represents (known) load increments, and
cess. {ai} and {Ax}denote (yet unknown) increments of accel-
erations and displacements. It is now necessary to take
DYNAMIC ANALYSIS advantage of the fact that {a} is essentially the second
In extending the method to the dynamic case, the derivative of {x}. Newmark’s #Lmethod[l3] appears to
masses are assumed to be lumped at the joints oniy. In be particularly convenient for this purpose. Using, for
addition, the effects of rotary inertia and damping, which example, a value of /3 = l/4, Newmark’s equations reduce
lead to complications in the Eulerian formulation, are to
neglected.
Plane fruntes. Under the assumptions stated, the equa- {Ax}= (AtKi”)+f (At)V’“)} +$(At)2{?K+1)} (38)
tions of motion of a planar framed structure can be
obtained from eqn 12 by simply inserting the necessary and
in& terms. Thus
{A~]=~(At)ii(K))+~(ntXi(“+l)) (3%
f~]~~x’)+VP-~~~=O (35)

in which (M] = mass matrix, and a dot on xi indicates a in which & = velocity. Substituting {_?x’}+{ti~ for
differentiation with respect to time. With rotary inertia {PK+‘)}in eqn (38), and solving for {AZ},
effects out of the picture, [A4Jis a diagonal matrix,
IAil’ ---+Ax}
(At)2 -&{P}-2li’9 (40)
M,
M*
Substituting this expression into eqn (37). and reananging,

Ml= (36) 4
~‘]~~}~~A~~ [M]{Axf = {AP)

‘& [M](.F~ f 2[M]{i’“)}. (41)

Two different computation procedures can now be


0
devised, one with iteration and one without. The latter
technique, which is herein referred to as the the linearized
in which Mm+,= Mn+2= ***= MB= 0 correspond to incremental procedure, involves the following steps:
translations of massless joints, as well as rotations of all
joints. 1. Evaluate [T] at t = tfi
Assume now that the configuration of the structure is 2. Consider a time increment Af, and calculate {AP).
546 C. ORANand A. KASSIMALI

Time t, set

Fig. 10. Response curves for fixedbeam.

3. Solve eqn (41) for {AX},then calculate {x(~+‘)}. The same problem was also solved by a three-
4. From eqn (40), determine {Ai}, then {,PK+‘)}. dimensional version of the method described above, and
5. From eqn (39). determine {ti}, then {i(K+‘)}. the results (not presented herein) were found to be in
close agreement with those shown in Fig. 10.
In the iterative technique, an approximate solution is
first obtained at f = fK+I by the linearized incremental
procedure outlined in the preceeding paragraph. This SUMMARY AND CONCLUSIONS
solution is then improved step by step, until the equations Some of the more important conclusions to be drawn
of motion are satisfied within a prescribed tolerance. In from this investigation are:
what follows, the subscript j is used to refer to the jth I. The proposed method of analysis proves to be highly
iteration cycle, while t is kept constant at I = fK+,. The accurate even when the deformations of the structure are
computational work that the jth cycle involves can be quite large, provided the Eulerian approach is used.
summarized as follows: 2. Judging from the limited numerical data presented,
1. From eqn (35). calculate the unbalanced joint forces: the Lagrangian approach can be expected to yield
generally satisfactory results for small or moderately
{AQX = IV - [Ml{~]i -U]i (42) large deformations.
3. Neglecting G,, G2, H in [f] does adversely affect the
2. Let {AAx},,{AAih, {Ufi denote the corrections behavior of the tangent stiffness matrix. By retaining
that {AQX implies for {Ax}, {&}, {Ai}. Then from eqn these terms, one may be able to reduce the total
VO), computational work through improved accuracy, although
the work required per load increment and/or per iteration
cycle may show a slight increase.
4. In most problems of practical interest, 6 is indeed
small in comparison with unity, and can be neglected.
3. Noting, from eqn (37), that However, it may be preferable to retain S in the
formulation, because (a) the required additional computa-
4 tional work is often negligibly small, and (b) in some
bliIAAxh +m [MlWxI, = IAQX WI problems (such as the steep truss of Fig. 8) S may be
somewhat large.
determine {AAx}~ 5. Critical points (whether of bifurcation or limit type)
4. If ,{AAx}~ obtained in the preceding step is suffi- can be detected conveniently and accurately by monitor-
ciently small in comparison with {x} according to a ing the determinant of the system tangent stiffness matrix.
prescribed tolerance, then stop iteration, increase t by At, 6. The observations summarized in items 14 are also
obtain a new linearized incremental solution, and return to applicable to the dynamic case. It may be added that the
step 1. Otherwise: advantages of the Lagrangian approach become much
5. From eqn (43), determine {AAi}, then {i}, and return more apparent in the dynamic case, especially when
to step 1. analyzing space problems.
Numerical solutions were obtained for the problem
shown in Fig. 10 by assuming that the system is initially at REFERENCES
rest and the force P is applied suddenly. Of the three 1. C. Oran, Tangent stiffness in plane frames. J. Strucr. Dick
solutions presented, the iterative one agrees closely with ASCE 99(ST6), 973-985 (1973).
the corresponding results reported by Weeks[14]. Note 2. C. Oran, Tangent stiffness in space frames. I. Sfruct. Div..
also that the linearized incremental solutions do tend to ASCE 99(ST6), 987-1001 (1973).
drift away, much like some of the less accurate solutions 3. C. Oran, Geometric nonlinearity in nonprismatic members. _I.
included in Week’s comparative study. SPUC~. Div., ASCE lOOfST7),1473-1488 (1974).
Large deformatins of framed structures under static and dynamic loads 541

4. S. Mau and R. H. Gallagher, A finite element procedure for structural analysis. Prog. Aeronnutical Sciences, 4. Macmil-
nonlinear prebuckling and initial postbuckling analysis, lan, New York (1964).
NASA Cr-1936 (1972). 10. A. Jennings, Frame analysis including change of geometry. J.
5. K. E. Bisshopp and D. C. Drucker, Large deflection of Sfrrtcf. Diu., ASCE 94(ST3),627-644 (196lQ.
cantilever beams. Quart. 1. Appl. Math. III(3). 272-275(1945). 11. G. H. Powell, Theory of nonlinear elastic structures. J. Srruct.
6. M. E. James, T. J. Kozik .and J. E. Martinez, Effect of h’u., ASCE 95(STl2), 2687-2701(1%9).
curvature on nonlinear frame analvsis. 1. Strucr. Dia., ASCE 12. R. Frisch-Fay, A new approach to the analysis of the
loo(sT7,, 1451-1457(1974). . deflection of thin cantilevers. .I.Appl. Mech. 28,87-90 (l%l).
7. A. M. Ebner and J. J. Ucciferro, A theoretical and numerical 13. T. P. Tung and N. M. Newmark, A review of numerical
comparison of elastic nonlinear finite element methods. integration methods for dynamic response of structures.
Comput. Struct. 2(5/6), 1043-1061(1972). Univ. Illinois Civil Eng. Studies, Srrucr. Rex Series, No. 69
8. H. C. Martin, On the derivation of stiffness matrices for the (1954).
analysis of large deflection and stability problems. Proc. Con!. 14. G. ’ Weeks, Temporal operators for nonlinear structural
Matrix Merhods in Srr~crurul Mechanics, Wright-Patterson dynamics problems. 1. E. M. Div., ASCE !Xt(EMS),1087-l104
Air Force Base, Ohio (1%5). (1972).
9. J. H. Argyris. Recent advances in matrix methods of

You might also like