- Email: [email protected]
beam components of arbitrary cross-section II. Shear stresses
Computers and Structures 80 (2002) 391–401 www.elsevier.com/locate/compstruc
Sensitivity analysis and optimization of truss/beam components of arbitrary cross-section II. Shear stresses V. Apostol, J.L.T. Santos *, M. Paiva Dep. of Mech. Eng., LEMAC, Instituto Superior T ecnico, Av. Rovisco Pais 1, 1096 Lisbon Codex, Portugal Received 10 January 2001; accepted 20 January 2002
Abstract This paper presents a general approach for detailed analysis and design optimization of arbitrary cross-sections of truss/beam built up structures. The approach allows arbitrary shape parametrization of 2-D cross-sections, as long as the coordinates of the contour vertices and their velocities are available, and is well suited for integration with existing CAD modelers and FEM analyzers. It leads to an inexpensive 2-D size/shape optimization in an alternative to costly 3-D shape optimizations, virtually impossible for real-life built up structures. Any composite multi-contour crosssection is first discretized with elementary triangles. Direct integration on the surface, using closed-form formulas, allows computation of the cross-section axial properties. Numerical integration on the boundary, along the line segments used to describe the contour, allows the computation of the shear properties. The power-series method is used to obtain the equilibrium equations and their governing linear warping system. The design sensitivities are calculated by the direct differentiation method requiring only backward substitutions on the triangular stiffness matrix. Numerical tests extensively verify the accuracy and the practical use of the formulation and implementation. 2002 Elsevier Science Ltd. All rights reserved. Keywords: Cross-section optimization; Design sensitivity; Structural optimization; Shear stress; Power-series
1. Introduction Although truss/beam like components are extensively used in civil and mechanical engineering structural applications, only very few authors have addressed the topic of design optimization of arbitrary cross-sections. Apostol and Santos [1,2] described an arbitrary parametric cross-section through geometrical shape variables, and expressed each vertex coordinate as a linear combination of these shape variables. The cross-section axial properties were computed and differentiated analytically with respect to the selected shape variables. These properties, representing finite el-
*
Corresponding author. E-mail address: [email protected] (J.L.T. Santos).
ement input quantities, were used as performance measures in the optimization model. However, previous research on design optimization of arbitrary cross-sections, has been focused on axial properties and stresses. Still, it is well known that the shear stresses induced by transverse shear loads and especially by torsion can reach significant values for space frames and therefore they can neither be neglected, nor poorly approximated. The analysis of the torsion problem may be performed by the semi-inverse method proposed by SaintVenant by the mid of the last century. The method is exact, provided that a warping function satisfying a differential equation and certain boundary conditions is known. A comprehensive description of the semi-inverse method can be found in Timoshenko [3]. One of the first papers addressing cross-section optimization of elastic bars under torsion seems to belong
0045-7949/02/$ - see front matter 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 2 ) 0 0 0 1 0 - X
392
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
to Banichuk [4], although shear stresses were not considered. Pilkey and Liu [5] avoided finite or boundary elements for the computation of the warping function by using the direct integration of the boundary integral equation instead. Their approach to the solution of the torsion problem was used by Schramm and Pilkey [6] for structural optimization in conjunction with shape description via B-splines. The optimization contained the torsion moment of inertia and the cross-section area only, whereas the transverse shear, the shear center and the shear stresses were not taken into consideration. The same authors extended the previous research to thinwalled beam theory, Schramm et al. [7], describing the shape with non-uniform rational B-splines (NURBS), Schramm and Pilkey [8,9]. A variational approach based on power-series, was introduced by Mindlin [10] for the computation of the warping function. This approach was further improved by Kosmatka [11,12], and it was selected in this paper as the basis for the cross-section analysis and shape optimization. Kosmatka used exact Gauss-quadrature integration for higher order polynomials, for numerical evaluation of the integrals on a surface, according to Dunavant [13]. In this paper, the power-series approach used by Kosmatka is used for numerical evaluation of the warping function. However, the Green’s formula is used to convert the surface integrals into contour integrals, see Press et al. [14]. The direct differentiation method is used to differentiate the equilibrium equations with respect to crosssection shape design parameters. The design sensitivities are finally used to demonstrate the validity of the approach for numerical shape optimization of truss/beam arbitrary cross-sections. Although two design parametrization approaches are developed and compared in this paper, the current method allows the use of any general parametrization like those provided by parametric and associative CAD systems.
2. Shape parametrization The aim of design parametrization is to relate the higher level shape variables with the lower level variables, vertex and control points, used to define the crosssection contours. Once the design parametrization is known, the design sensitivities of the vertex coordinates with respect to the shape variables may be easily computed. This sensitivity information is referred as the cross-section velocity field. Apostol and Santos [2] used a direct linear parametrization, assuming each vertex coordinate as being a linear combination of shape variables. The approach led to an important simplification, since the constant coef-
ficients of the linear combination represent in fact the required velocity field. Direct non-linear parametrizations, allowing nonlinear dependencies between the shape variables and the vertex coordinates is also sometimes used in practice. Recently, Langelaan and Livne [15] parametrized the distance along the axis and the radii in transversal directions of slender curved hollow tubes, providing analytically the cross-section velocity field. A numerical example is used in this paper to show the disadvantages of this method, mainly related to the non-smoothness of the final shape. Indirect parametrization, using interpolation curves for the contour definition, assures the smoothness of the cross-section shape. Braibant and Fleury [16] introduced the use of B-splines for shape structural optimization, whereas Schramm and Pilkey [8,9] used NURBS for shape parametrization. However, interpolation curves may also make it difficult to achieve desired regular shapes. If proper control dimensions (radii, angles, etc.) are selected as shape variables, the final shape remains regular and smooth as desired. Such direct parametrization strategy in connection with advanced CAD systems seems to be attractive for structural optimization, see Lindby and Santos [17].
3. Axial properties To calculate the axial properties and their design sensitivities with respect to a set of shape parameters b, Ref. [2] describes a complex cross-section as the sum of closed polygonal contours. In the initial reference system, the X-axis is along the beam, and the beam transversal plane ZOY is used for locating all vertices used to define the closed contour making-up the crosssection. The cross-section centroid ðzG ; yG Þ and the axial moments of inertia related to the reference system ZOY allow the computation of the principal axes of e 2GX e 3 , with the X e 1 -axis remaininertia denoted by X 1 ing along the beam. Consequently, I2 and I3 become the maximum and the minimum axial moments of inertia. Fig. 1 shows the dependency between the two reference systems. The principal angle of inertia ap is conventionally defined as the angle made by the arbitrary Z-direction with the principal maximum axis of e 2 . All vertex coordinates are translated and inertia X rotated to the principal axes of inertia.
1 Notation denotes a measure related to the principal system of inertia.
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
393
Express the warping functions as the double infinite power-series 3 wk ð~x2 ; ~x3 Þ ¼
1 X 1 X i¼0
ðkÞ cij ~xi2~xj3
ð2Þ
j¼0
which are truncated in reality to a reasonable high order n, with i þ j 6 n. In vector notation the components of the warping function are defined as wk ð~x2 ; ~x3 Þ ¼ N cðkÞ
Fig. 1. Cross-section subject to shear loads and torque.
The design sensitivities of all axial properties explicitly known in terms of the vertex coordinates, can be obtained via analytical differentiation, see Apostol and Santos [2]. The axial stress is expressed in terms of a constant stretching force, and two linear bending moments, and depends upon the coordinates (~x2 ; ~x3 ) of the cross-section point under analysis. The design sensitivities of the axial stress also need their corresponding velocity field (~x02b ; ~x03b ). 2 4. Shear properties
A significant simplification in the formulation proposed by Kosmatka [12] is achieved introducing two vectors of size m, namely g, l, containing the corresponding exponents of ~x2 , ~x3 in N n N ¼ ~xg21 ~xl3 1 ; ~xg22 ~xl3 2 ; ~xg23 ~xl3 3 ; ~xg24 ~xl34 ; ~xg25 ~xl35 ; . . . ; ~xg2m 2 ~xl3 m 2 ; o ~xg2m 1 ~xl3m 1 ; ~xg2m ~xl3 m : ð4Þ For the shear properties under investigation, only the first derivatives of the shape functions with respect to coordinates are required oNi ¼ gi~x2gi 1~xl3i ; o~x2
ð1Þ
ð2Þ
s1j ¼ s1j h þ s1j
4.1. Power-series method
s12
wð~x2 ; ~x3 Þ ¼ w1 ð~x2 ; ~x3 Þh þ w2 ð~x2 ; ~x3 Þ
ð1Þ
i ¼ 1; . . . ; m
ð5Þ
P2 ð3Þ P3 þ s1j ; EI3 EI2
j ¼ 2; 3
ð6Þ
where ð1Þ
ow1 ow1 ð1Þ
~x3 ; s13 ¼ G þ ~x2 o~x2 o~x3 ow2 m 2 ow2 ð2Þ s13 ¼ G
~x2 ~x23 ;
m~x2~x3 ¼G o~x2 2 o~x3 ow3 ow3 m 2 ð3Þ
m~x2~x3 ; s13 ¼ G
~x3 ~x22 ¼G o~x2 o~x3 2
s12 ¼ G ð2Þ
s12
ð3Þ
P2 P3 þ w3 ð~x2 ; ~x3 Þ EI3 EI2
oNi ¼ li~xg2i ~xl3 i 1 ; o~x3
Shear stress formulas may be expressed in terms of the three components of the warping function, corresponding to the torsion and the two transverse shear loads, Kosmatka [12],
The calculation of the warping function w via powerseries leads to the shear properties of the cross-section, namely the torsion moment of inertia I1 , the cross-section twist rate h, and the shear center ð~xs2 ; ~xs3 Þ. The cross-section in-plane shear stresses s12 , s13 acting along the principal directions may be determined next. Finally, the design sensitivities of the shear properties, required for design optimization, are obtained by analytical differentiation.
The warping function describes the warping out of the cross-section plane due to transverse loads P2 , P3 applied at the centroid and along the principal directions, and to a torque M1 applied along the beam axis, see Fig. 1. Define the warping function as, Kosmatka [12],
ð3Þ
ð7Þ The three sets of unknown warping coefficients deðkÞ fined in Eq. (3), ci ; i ¼ 1; . . . ; m, can be determined solving the linear system resultant from the virtual work principle, dU ¼ dW . According to Kosmatka [12], after replacing the strains and stresses by their values, into the virtual work principle, the coefficients cðkÞ are obtained from the following linear system of equations
where the twist rate h is a function of M1 , P2 , P3 . 2
Design sensitivities are written as ð Þ0b .
3 ð Þk , or ð ÞðkÞ , k ¼ 1, 2, 3, refers to one component of the warping function.
394
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
KcðkÞ ¼ FðkÞ
ð8Þ
ðkÞ ðkÞ FðkÞ FðkÞ w Fc ¼ ðFw Fc ÞQ
ð9Þ
the stiffness matrix defined as T T ) Z ( oN oN oN oN K ¼ GL þ dA o~x2 o~x2 o~x3 o~x3 A and the force matrices given by Z
Z Fw ¼ EL 0; ~x2 NT dA; ~x3 NT dA
Fc ¼ GL
Z (
2
A
~x2 4
~x3
ð10Þ
ð11Þ
A
oN o~x3
T T ) oN ; o~x2
3
m~x2~x3
2m ð~x23 ~x22 Þ 5dA
2m ð~x22 ~x23 Þ
m~x2~x3
ð12Þ
with Qð1Þ ¼ f1; 0; 0gT ;
nv Z X F2 d~x3 þ F3 d~x2 ¼ F2 d~x3 þ F3 d~x2
C
with the load vector
A
I
Qð2Þ ¼ f0; 1; 0gT ;
Qð3Þ ¼ f0; 0; 1gT ð13Þ
i¼1
ð15Þ
Ci
For a comparison between the surface and contour integration approaches, consider a polynomial of order n ¼ 11. Whereas the exact integration on the contour requires 11 Gauss points along a line, the associated triangle made of the two line end vertices, and the origin, must be integrated using 79 Gauss points for the surface integral. Savings of up to 85% in run-time follow, if contour integrals are used. Regardless of the approach used to carry-out the required integrations, the stiffness matrix and the load vector have to be created by assembling the corresponding values of the base elements, triangles if surface integrals are used, straight lines if contour integrals are used. In the classical finite element method, one finite element depends on a subset of degrees of freedom, and therefore the elementary stiffness matrix adds terms only at particular locations inside the overall stiffness matrix. In the power-series method, any elementary stiffness matrix depends on all warping coefficients cðkÞ , and thus contributes to all locations of the global stiffness matrix. The assembly process degenerates here into a simple fullmatrix summation over all Gauss points, and no sparsity can be exploited.
4.2. Numerical integration on the surface
4.4. Parametric form of a straight line
The integration of the stiffness matrix and force vector is performed on the cross-section domain. The cross-section is viewed as a sum of triangles defined by the corresponding centroid and two adjacent vertices. For complex shapes, with a larger number of vertices, higher order polynomials must be used. Hence, numerical integration on the cross-section domain requires to build the stiffness matrix as proposed by Kosmatka [12]. An effective alternative for the integration of the stiffness matrix and force vector is here proposed. The domain integrals can be transformed into contour integrals, with significant gains in computational time, and allowing for simplified computation of the design velocities required for cross-section shape design sensitivity analysis.
The parametric form of a straight line defined in the e 2GX e 3 system is X
4.3. Numerical integration on the contour Green’s formula states that given two functions F2 and F3 , the following transformation holds I Z oF2 oF3
dA ¼ ð14Þ F2 d~x3 þ F3 d~x2 o~x2 o~x3 A C If the contour C is a polygon with nv vertices, and nv straight sides, the contour integral can be written as
~x2 an ¼ a0 n þ a1 ) d~x2 ¼ a0 dn ~x3 bn ¼ b0 n þ b1 ) d~x3 ¼ b0 dn
ð16Þ
where a0 , a1 , b0 , b1 are constant coefficients, and
1 6 n 6 1 is the parametric variable along the line. Auxiliary variables (an ; bn ) are shortcuts for the longer formulas which follow, and they simply express the coordinates at a certain location n along the line. Considering two adjacent vertices of coordinates ð~xi2 ; ~xi3 Þ and ð~xðiþ1Þ2 ; ~xðiþ1Þ3 Þ, i ¼ 1; . . . ; nv , coefficients ~xðiþ1Þ2 ~xi2 ; 2 ~xðiþ1Þ3 ~xi3 b0 ¼ ; 2
a0 ¼
~xðiþ1Þ2 þ ~xi2 ; 2 ~xðiþ1Þ3 þ ~xi3 b1 ¼ 2
a1 ¼
are determined for each segment of the contour C. For i ¼ nv , the last and the first vertices of the closed contour determine the parametric coefficients. In order to keep formulas easier to read, no script on coefficients a, b indicates a contour segment. The summation introduced in Eq. (15) is implicitly assumed. Design derivatives can be computed straightforward if the velocities ~x02b , ~x03b are known. Hence,
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
a00b ¼ ~x02b
b00b ¼ ~x03b
~x0ðiþ1Þ2b ~x0i2b a0nb
2 ¼ a00b n þ a01b
~x0ðiþ1Þ3b ~x0i3b b0nb
;
;
2 ¼ b00b n þ b01b
a01b ¼
~x0ðiþ1Þ2b þ ~x0i2b 2
;
ð17Þ
h¼
1 P2 P3 M1 a2
a3 EI3 EI2 a1 ðkÞ
b01b ¼
~x0ðiþ1Þ3b þ ~x0i3b 2
;
ð18Þ
395
ð24Þ
ðkÞ
with s12 , s13 given by Eqs. (6) and (7). Replacing the shear stresses into Eq. (23), for k ¼ 1, follows Z Z Z ow1 ow1 2 2 ~x2 a1 ¼ G
~x3 dA þ ~x3 dA þ ~x2 dA o~x3 o~x2 A A A ð25Þ
4.5. Solving the linear warping system An element of the symmetric stiffness matrix Kij , i ¼ 1; . . . ; m, j ¼ i; . . . ; m, introduced in Eq. (10), can be written as Z g þg 2 l þl g þg l þl 2 dA Kij ¼ GL gi gj~x2i j ~x3i j þ li lj~x2i j ~x3 i j ZA oF2 oF3
dA ¼ GL o~x2 o~x3 A ð19Þ Identifying now the two functions of the Green’s formula it follows gi gj g þg 1 l þl ~x i j ~x3 i j ; gi þ gj 1 2 li lj g þg l þl 1 ~x i j ~x3i j F3 ¼
li þ lj 1 2 F2 ¼
ð20Þ
Green’s formula (14), and the parametric form of the line (16) assure that Z 1 Kij ¼ GL ðF2 b0 þ F3 a0 Þdn ð21Þ
1
or Z
ð22Þ
4.6. Computation of shear properties
A
and the twist rate as
5. Design sensitivity analysis The parametric definition of the cross-section assures the connection between the vertex coordinates and the selected design parameters. The design sensitivities P0ab of the axial properties ~2 ; x ~3 ; r11 g I2 ; I3 ; ap ; x
The linear symmetric warping system (8) can be solved now for the three sets of coefficients cðkÞ . Therefore, a two-step based solver, using forward elimination and backward substitution, is more suitable than an iterative solver.
Define the warping constants, Kosmatka [12], Z ðkÞ ðkÞ ~x2 s13 ~x3 s12 dA ak ¼
where the axial moments of inertia are identified with the last two integral terms in Eq. (25). Eq. (26) clearly shows the error involved in approximating I1 ffi I2 þ I3 , since the first integral term on the right hand side is neglected. In fact, this expression is correct only for circular cross-sections. Warping coefficients a2 , a3 can be used to locate the shear center a3 a2 ð~xs2 ; ~xs3 Þ ¼ ð27Þ ;
EI2 EI3
Pa ¼ fzT ; yT ; AT ; SzO ; SyO ; A; IzzO ; IyzO ; IyyO ; IzzG ; IyzG ; IyyG ; zG ; yG ;
1
gi gj g þg 1 l þl a i j bn i j b0 gi þ gj 1 n
1 ! li lj gi þgj li þlj 1 a bn a0 dn
li þ lj 1 n
Kij ¼ GL
The warping constant a1 is nothing but the torsional rigidity GI1 . Hence, the torsion moment of inertia I1 becomes Z ow ow ~x2 1 ~x3 1 dA þ I2 þ I3 ð26Þ I1 ¼ o~x3 o~x2 A
ð23Þ
were presented in Apostol and Santos [2]. In this section the design sensitivities P0sb of the shear properties below are obtained, Ps ¼ fI1 ; h; ~xs2 ; ~xs3 ; s12 ; s13 g: Computation of the shear properties requires the solution of a linear system of equations. The analogy with the finite element method allows a parallel between the design sensitivity methods developed for the linear finite element method and for the power-series method. Whereas the displacement field is the unknown solution in the finite element method, the present method uses the warping coefficients cðkÞ as the solution of the governing ðkÞ0 system of equations. Hence, the design sensitivities cb are the basis for the design sensitivities of the shear properties.
396
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
Since the entire design sensitivity field is required, the direct differentiation method is here preferred against the adjoint variable method. Differentiating wk ¼ N cðkÞ with respect to design, ðkÞ0
w0kb ¼ N0b cðkÞ þ Ncb
ð28Þ
Differentiating the first derivative of wk with respect to design follows ow0kb oN0b ðkÞ oN ðkÞ0 ¼ c þ c ; o~xj b o~xj o~xj
j ¼ 2; 3
ð29Þ
Taking advantage of the closed-form formulas (5), the above equation becomes m 0 X ow0kb ðkÞ g 1 l g 1 l ðkÞ0 ¼ gi~x2i ~x3i ci þ gi~x2i ~x3 i cib b o~x2 i¼1 ð30Þ m 0 X 0 owkb ðkÞ gi li 1 gi li 1 ðkÞ0 ¼ li~x2 ~x3 ci þ li ~x2 ~x3 cib b o~x3 i¼1 or arranged in compact form m nh i X ow0kb g 2 l g 1 l 1 ¼ gi ðgi 1Þ~x2i ~x02b ~x3i þ gi li ~x2i ~x3i ~x03b o~x2 i¼1 o ðkÞ g 1 l 0ðkÞ ci þ gi ~x2i ~x3 i cib m nh i X ow0kb g 1 l 1 g l 1 ¼ li gi ~x2i ~x02b ~x3i þ li ðli 1Þ~x2i ~x3 i ~x03b o~x3 i¼1 o ðkÞ ðkÞ0 ci þ li ~xg2i ~xl3 i 1 cib
ð31Þ These formulas, and Eqs. (6)–(27) show that any design sensitivity of a shear property Ps can be later computed if additionally to the previous computations, the veloc~03b Þ, and the warping derivatives cðkÞ0 ity field ð~ x02b ; x are b available. Recall that the warping coefficients cðkÞ are computed from the linear system KcðkÞ ¼ FðkÞ , by differentiating it with respect to design ðkÞ0
K0b cðkÞ þ Kcb
ðkÞ0
¼ Fb
ð32Þ
Moving the first term from the left side to the right side, results in ðkÞ0
Kcb
ðkÞ0
¼ Fb K0b cðkÞ
ð33Þ
This system has the same stiffness matrix as the original, and only backward substitutions are required to comðkÞ0 pute cb . Derivatives of the stiffness matrix, load vectors, and warping constants ak , require integration on the contour, i.e. the velocity field in any Gauss point along the contour. This information is obtained according to Section 4.4.
If surface integration is used instead of contour integration, as proposed by Kosmatka, no velocities would be available in the interior of the contour. This disadvantage is even more severe than the higher runtime used by the surface integration, and ultimately influenced the use of contour integration. The algorithm used for the design sensitivity calculation of the shear properties is as follows: • Compute the stiffness matrix K, and the loads FðkÞ . The stiffness matrix is scaled and decomposed after the Gauss forward elimination. • Compute design sensitivities of the stiffness matrix ðkÞ0 K0b , and loads Fb . They are used to build the pseudoloads associated with the right side of Eq. (33). • Solve 3 np times the system (33), where np is the number of design variables, through backward substitutions on the unchanged, decomposed stiffness ðkÞ0 matrix, and get the coefficients cb . • Integrate the design sensitivities of the warping constants a0kb . • Compute the design sensitivities of the shear proper0 ties I1b , h0b , ð~x0 x0 s2b ; ~ s3b Þ. • Compute the design sensitivities of the shear stresses ðkÞ0 sb , for any desired point along the contour. It is worth noticing that the theoretical results presented in this paper together with those published in Apostol and Santos [2] enable the computation of the axial and shear design sensitivity contributions of any constraint that can be written as a function of the axial and shear properties included in both vectors Pa and Ps .
6. Optimization of a thin-walled cross-section subject to shear stress constraints This section presents the optimization of a thin-walled cross-section, with respect to shear stress constraints, for two different parametrizations, direct vertex parametrization, and indirect vertex parametrization using NURBS. 6.1. Direct vertex parametrization A thin-walled cross-section is parametrized using direct non-linear relationships between vertex coordinates and design variables, see Fig. 2. The initial crosssection shape is described as a regular 20-sided polygon, with inner and outer vertices expressed in terms of 18 geometrical dimensions and the wall thickness. The cross-section is made of steel with a Young’s modulus of 2.1e þ 05 MPa, and the Poisson’s ratio is 0.3. The loading consists of a torque of Mx ¼ 1e þ 05 N mm, and two shear forces Py ¼ Pz ¼ 1e þ 04 N.
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
397
Table 1 Optimization results of a thin-walled cross-section
Fig. 2. Definition and loading of a thin-walled cross-section.
The optimization problem to be solved is the minimization of the cross-section area subject to shear stress constraints along the cross-section boundary. The shear stresses are defined at nine equidistant points within each boundary segment. A total of 234 stress performances are included in the optimization model. The shear stress constraints are defined as s ¼ ðs2xy þ s2xz Þ1=2 , where smax ¼ 100 MPa. The design parametrization requires the wall thickness to be constant along the wall cross-section contour. The horizontal vertex coordinates must be kept constant, and the vertical vertex coordinates are allowed to vary linearly with respect to the design variables, for the inner vertices, and non-linearly for the outer vertices. The analytical expressions for the vertices may be easily obtained, and therefore they are not presented here. Whenever a new design is investigated, the vertex coordinates must be recomputed as well as the nodal velocity field. The direct differentiation method described in Section 5 is used for nodal velocity computations. The output of these computations is used for the cross-section shear related geometric and performance computations. The velocity field within the contour is interpolated linearly from the nodal velocities. The optimization problem is solved with a sequential quadratic algorithm, see Apostol [18] for more details. Table 1 shows the initial and the final results after 20 iterations. The initial design is feasible since all shear stresses are smaller than 100 MPa. The 15% structural reserve allows the area to decrease from 453.18 to 351.58 mm2 . At the 20th iteration, 16 out of 234 constraints become active, i.e. the corresponding shear stresses are close to 100 MPa. Fig. 3 shows the first and the last shapes, for comparison. The optimizer did not converge yet, but no significant changes occurred after the 20th iteration. Fig. 4 presents the objective function iteration history. However, the optimum shape is not likely to be suitable
DV
Initial (mm)
Final (mm)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
9.270 17.633 24.270 28.531 30.000 28.531 24.270 17.633 9.270
9.270
17.633
24.270
28.531
30.000
28.531
24.270
17.633
9.270
10.057 19.675 30.381 35.018 34.507 29.810 27.302 19.586 15.365
5.303
12.582
32.162
36.486
41.115
39.878
24.541
17.982
5.515
Properties Thickness t (mm) Area A (mm2 ) Tors. mom. I1 (mm4 ) Max. s (MPa)
2.000 453.18 3.842eþ05 85.267
1.395 351.58 3.240eþ05 100.020
for manufacturing since it has a very irregular contour. A new parametrization, typical of many CAD systems, is used next to overcome the above problems and verify the appropriateness of the proposed formulation for indirect vertex parametrization. This example clearly shows the application of the theoretical developments presented in Sections 4 and 5. Despite the highly non-linear relationships between the shear stress constraints and the design variables, still, regular optimization paths and convergence were achieved, demonstrating that the power-series method is a practical tool for structural optimization. 6.2. Indirect vertex parametrization using NURBS The approach of direct vertex definition with shape variables was earlier intensively used. Nowadays, interpolation curves are preferred in modern shape optimization, having at least two advantages: they are controlled by fewer design variables, and they assure the smoothness required for manufacturing. In particular NURBS were used by Schramm and Pilkey [8,9] for design parametrization and integration with structural shape optimization. The thin-walled structure illustrated in Fig. 2 is now parametrized with two NURBS, see Fig. 5. The inner vertices are defined by the two independent sets of control points: p1 , p2 , p3 , p4 , p5 , p5 , and p5 , p6 , p7 , p8 , p1 , p1 , with coordinates ðzj ; yj Þ given under the column ‘Initial’
398
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
Fig. 3. Initial and final shapes of a thin-walled cross-section.
in Table 2. Remark that the last vertex has to be entered twice. This is a characteristic of B-splines. The polynoTable 2 Optimization results using NURBS
Fig. 4. Optimization history of a thin-walled cross-section.
Fig. 5. Definition of the thin-walled cross-section using NURBS.
DV
Initial (mm)
Final (mm)
z1 y1 w1
30.00 0.00 1.00
30.00a 0.00a 1.00a
z2 y2 w2
30.00 19.60 1.00
32.00a 29.68 1.92
z3 y3 w3
0.00 34.00 1.00
4.49 37.46 1.25
z4 y4 w4
30.00 17.60 1.00
31.96 10.09 1.96
z5 y5 w5
30.00 0.00 1.00
30.00a 0.00a 1.00a
z6 y6 w6
30.00
17.60 1.00
32.00a
40.00a 1.19
z7 y7 w7
0.00
34.00 1.00
20.00a
40.00a 1.96
z8 y8 w8
30.00
17.60 1.00
20.00a
5.00a 1.99
2.00 456.40 3.977eþ05
1.21 304.81 2.154eþ05
83.29
100.03
Properties Thickness t (mm) Area A (mm2 ) Tors. mom. I1 (mm4 ) Max. s (MPa) a
Design variable bounds imposed for constructive reasons.
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
399
Fig. 6. Initial and final shapes of inner vertices using NURBS.
mial order is fixed to n ¼ 3, whereas each B-spline is transformed into a polygon with m ¼ 31 vertices. Fig. 6 left shows that, due to the non-uniform knot sequence, the intervals between vertices are not equal. The shape is approximately circular. However, the initial cross-section area defined through NURBS differs by less than 1% relatively to the previous optimization with direct vertex parametrization. The outer vertices are analytically computed such that the thickness t remains constant over the crosssection. Each horizontal flange does not change its dimensions, except for the thickness. Control points p1 and p5 are kept constant during the optimization process. The other six control points define 18 design variables (yj ; zj ; wj ), where wj are the weighting factors. The two curves act independently, and therefore, sharp corners could appear at the connection between the upper and the lower curves. The bounds on the adjacent control points may control this undesirable effect. The thickness t is the 19th design variable. The optimization problem to be solved is the same as that of Section 6.1. It is required to minimize the cross-
section area subject to shear stress constraints, along the cross-section boundary. The shear stresses are defined in this case as three shear stresses at equally spaced points along each segment, on the outer contour. A total of 198 constraints of the type s 6 smax , with smax ¼ 100 MPa, are specified. The convergence of the optimization process was rather slow, requiring 40 iterations. The final results of the model are given under column ‘Final’ in Table 2. Some of the design variables reached their bounds. This fact shows that slightly improved shapes could be found if the bounds were relaxed. However, in practical optimization the objective is to find meaningful shapes rather than a mathematical optimum. For comparison, Figs. 6 and 7 present the initial and the final shapes of the inner contour, and of the entire cross-section, accordingly. Fig. 8 shows the optimization history of the crosssection area. The jump in the objective function value observed at the 21st iteration was caused by the stop of the optimizer after the 20th iteration. After restarting the optimizer, the second order Hessian was lost, and a poorer solution was reached at the 21st iteration.
Fig. 7. Initial and final shapes using NURBS.
400
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
Fig. 8. Optimization history using NURBS.
Although the optimizer did not converge even after 40 iterations, since the variations in design were at least four orders of magnitude smaller than the units in millimeters, the process was stopped. Some comments are due to this point, regarding the optimization of the beam cross-section. As already verified in classical shape optimization, the parametrization using NURBS leads to a final smooth shape, whereas the parametrization using piecewise linear segments leads to an irregular final design, see Fig. 3. Also, the smoother parametrization associated with the NURBS allows convergence to a significantly better design than the piecewise parametrization of the boundary, leading to final cross-section areas of 304.81, and 351.58 mm2 , respectively. The optimal shape associated with the NURBS boundary parametrization illustrates better the optimal material distribution required to reach the highest stiffness under the shear constraints, by aligning the material preferably along two quasi straight-sided segments, almost parallel to each other, and to the 45 oriented shear force. A final comment on the selection of the NURBS weighting factors is appropriate. By defining these weighting factors as design variables it was possible to search for the optimum within a large design space. In Fig. 7a rounded curve at the right-upper side of the cross-section is associated with a weighting factor w3 ¼ 1:25, whereas a sharp corner at the lower-left side corresponds to a weighting factor w7 ¼ 1:96.
7. Conclusions A general approach for detailed analysis and design optimization of cross-sections of arbitrary shape of truss/beam built up structures is proposed. The approach allows for both direct and indirect parametric
definition of an arbitrary cross-section, and may be easily integrated with the parametrization supported by state-of-the-art CAD systems. The truss/beam components are subjected to axial and transversal forces, and to bending and torsion moments. Both shear properties and stresses are considered in both the analysis and the optimization process. The power-series method is used for the solution of the shear problem, and boundary integration is proposed for carrying-out the integrations required for the solution of the warping equations. A method for computing the design sensitivities of the shear properties and stresses with respect to crosssection design parameters is developed. This method uses direct differentiation of the warping equations and only requires backward substitutions on the triangular stiffness matrix, making it a quite attractive alternative for structural shape optimization of truss/beam crosssections of built up structures. A numerical example illustrates the applicability of the method to design optimization of a real-life crosssection.
References [1] Apostol V, Santos JLT. Design sensitivity analysis and optimization of truss/beam structures with arbitrary crosssections. In: Olhoff N, Rozvany G, editors. WCSMO-1 Proceedings of the First World Congress of Structural and Multidisciplinary Optimization, 28 May–2 June, Goslar, Germany. Oxford: Pergamon Press; 1995. p. 353–8. [2] Apostol V, Santos JLT. Sensitivity analysis and optimization of truss/team components of arbitrary cross-section. I––Axial stresses. Comput Struct 1996;58(4):727–37. [3] Timoshenko SP, Goodier JN. Theory of elasticity. 3rd ed. New York: McGraw-Hill; 1982. [4] Banichuk NV. Optimization of elastic bars in torsion. Int J Solids Struct 1976;12:275–86.
V. Apostol et al. / Computers and Structures 80 (2002) 391–401 [5] Pilkey WD, Liu Y. Field theory: A two-dimensional case for not using finite or boundary elements. Finite Elem Anal Des 1993;13(2/3):127–36. [6] Schramm U, Pilkey WD. Structural shape optimization for torsion problem using direct integration and B-splines. Comput Meth Appl Mech Eng 1993;107:251–68. [7] Schramm U, Kitis L, Kang W, Pilkey WD. On the shear deformation coefficient in beam theory. Finite Elem Anal Des 1994;16:141–62. [8] Schramm U, Pilkey WD. The coupling of geometric descriptions and finite elements using NURBs––a study in shape optimization. Finite Elem Anal Des 1993;15: 11–34. [9] Schramm U, Pilkey WD. Optimal shape design for thin walled beam cross-sections. Int J Numer Meth Eng 1994;37(23):4039–58. [10] Mindlin RD. Solution of St. Venant’s torsion problem by power series. In: Solid structures, vol. 11. Oxford: Pergamon Press; 1975. p. 321–8. [11] Kosmatka JB. On the behavior of pretwisted beams with irregular cross-sections. ASME––J Appl Mech 1992;59(1): 146–52.
401
[12] Kosmatka JB. Flexure–torsion behavior of prismatic beams, part I: Section properties via power series. AIAA J 1993;31(1):170–9. [13] Dunavant DA. High degree efficient symmetrical Gaussian quadrature rules for the triangle. Int J Numer Meth Eng 1985;21:1129–48. [14] Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical recipes. Cambridge: Cambridge University Press; 1989. [15] Langelaan JW, Livne E. Analytic sensitivities and design oriented structural analysis for airplane fuselage shape synthesis. Comput Struct 1997;62(3):505–19. [16] Braibant V, Fleury C. Shape optimal design using Bsplines. Comput Meth Appl Mech Eng 1984;44:247–67. [17] Lindby T, Santos JLT. 2-D and 3-D shape optimization using mesh velocities to integrate analytical sensitivities with associative CAD. Struct Optimization 1997;13(4): 213–22. [18] Apostol V. Design sensitivity analysis and optimization of geometrically nonlinear structures in the pre- and postbuckling domains. PhD thesis, Instituto Superior Tecnico Lisbon, Portugal, 1998.
Sensitivity analysis and optimization of truss/beam components of arbitrary cross-section II. Shear stresses V. Apostol, J.L.T. Santos *, M. Paiva Dep. of Mech. Eng., LEMAC, Instituto Superior T ecnico, Av. Rovisco Pais 1, 1096 Lisbon Codex, Portugal Received 10 January 2001; accepted 20 January 2002
Abstract This paper presents a general approach for detailed analysis and design optimization of arbitrary cross-sections of truss/beam built up structures. The approach allows arbitrary shape parametrization of 2-D cross-sections, as long as the coordinates of the contour vertices and their velocities are available, and is well suited for integration with existing CAD modelers and FEM analyzers. It leads to an inexpensive 2-D size/shape optimization in an alternative to costly 3-D shape optimizations, virtually impossible for real-life built up structures. Any composite multi-contour crosssection is first discretized with elementary triangles. Direct integration on the surface, using closed-form formulas, allows computation of the cross-section axial properties. Numerical integration on the boundary, along the line segments used to describe the contour, allows the computation of the shear properties. The power-series method is used to obtain the equilibrium equations and their governing linear warping system. The design sensitivities are calculated by the direct differentiation method requiring only backward substitutions on the triangular stiffness matrix. Numerical tests extensively verify the accuracy and the practical use of the formulation and implementation. 2002 Elsevier Science Ltd. All rights reserved. Keywords: Cross-section optimization; Design sensitivity; Structural optimization; Shear stress; Power-series
1. Introduction Although truss/beam like components are extensively used in civil and mechanical engineering structural applications, only very few authors have addressed the topic of design optimization of arbitrary cross-sections. Apostol and Santos [1,2] described an arbitrary parametric cross-section through geometrical shape variables, and expressed each vertex coordinate as a linear combination of these shape variables. The cross-section axial properties were computed and differentiated analytically with respect to the selected shape variables. These properties, representing finite el-
*
Corresponding author. E-mail address: [email protected] (J.L.T. Santos).
ement input quantities, were used as performance measures in the optimization model. However, previous research on design optimization of arbitrary cross-sections, has been focused on axial properties and stresses. Still, it is well known that the shear stresses induced by transverse shear loads and especially by torsion can reach significant values for space frames and therefore they can neither be neglected, nor poorly approximated. The analysis of the torsion problem may be performed by the semi-inverse method proposed by SaintVenant by the mid of the last century. The method is exact, provided that a warping function satisfying a differential equation and certain boundary conditions is known. A comprehensive description of the semi-inverse method can be found in Timoshenko [3]. One of the first papers addressing cross-section optimization of elastic bars under torsion seems to belong
0045-7949/02/$ - see front matter 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 2 ) 0 0 0 1 0 - X
392
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
to Banichuk [4], although shear stresses were not considered. Pilkey and Liu [5] avoided finite or boundary elements for the computation of the warping function by using the direct integration of the boundary integral equation instead. Their approach to the solution of the torsion problem was used by Schramm and Pilkey [6] for structural optimization in conjunction with shape description via B-splines. The optimization contained the torsion moment of inertia and the cross-section area only, whereas the transverse shear, the shear center and the shear stresses were not taken into consideration. The same authors extended the previous research to thinwalled beam theory, Schramm et al. [7], describing the shape with non-uniform rational B-splines (NURBS), Schramm and Pilkey [8,9]. A variational approach based on power-series, was introduced by Mindlin [10] for the computation of the warping function. This approach was further improved by Kosmatka [11,12], and it was selected in this paper as the basis for the cross-section analysis and shape optimization. Kosmatka used exact Gauss-quadrature integration for higher order polynomials, for numerical evaluation of the integrals on a surface, according to Dunavant [13]. In this paper, the power-series approach used by Kosmatka is used for numerical evaluation of the warping function. However, the Green’s formula is used to convert the surface integrals into contour integrals, see Press et al. [14]. The direct differentiation method is used to differentiate the equilibrium equations with respect to crosssection shape design parameters. The design sensitivities are finally used to demonstrate the validity of the approach for numerical shape optimization of truss/beam arbitrary cross-sections. Although two design parametrization approaches are developed and compared in this paper, the current method allows the use of any general parametrization like those provided by parametric and associative CAD systems.
2. Shape parametrization The aim of design parametrization is to relate the higher level shape variables with the lower level variables, vertex and control points, used to define the crosssection contours. Once the design parametrization is known, the design sensitivities of the vertex coordinates with respect to the shape variables may be easily computed. This sensitivity information is referred as the cross-section velocity field. Apostol and Santos [2] used a direct linear parametrization, assuming each vertex coordinate as being a linear combination of shape variables. The approach led to an important simplification, since the constant coef-
ficients of the linear combination represent in fact the required velocity field. Direct non-linear parametrizations, allowing nonlinear dependencies between the shape variables and the vertex coordinates is also sometimes used in practice. Recently, Langelaan and Livne [15] parametrized the distance along the axis and the radii in transversal directions of slender curved hollow tubes, providing analytically the cross-section velocity field. A numerical example is used in this paper to show the disadvantages of this method, mainly related to the non-smoothness of the final shape. Indirect parametrization, using interpolation curves for the contour definition, assures the smoothness of the cross-section shape. Braibant and Fleury [16] introduced the use of B-splines for shape structural optimization, whereas Schramm and Pilkey [8,9] used NURBS for shape parametrization. However, interpolation curves may also make it difficult to achieve desired regular shapes. If proper control dimensions (radii, angles, etc.) are selected as shape variables, the final shape remains regular and smooth as desired. Such direct parametrization strategy in connection with advanced CAD systems seems to be attractive for structural optimization, see Lindby and Santos [17].
3. Axial properties To calculate the axial properties and their design sensitivities with respect to a set of shape parameters b, Ref. [2] describes a complex cross-section as the sum of closed polygonal contours. In the initial reference system, the X-axis is along the beam, and the beam transversal plane ZOY is used for locating all vertices used to define the closed contour making-up the crosssection. The cross-section centroid ðzG ; yG Þ and the axial moments of inertia related to the reference system ZOY allow the computation of the principal axes of e 2GX e 3 , with the X e 1 -axis remaininertia denoted by X 1 ing along the beam. Consequently, I2 and I3 become the maximum and the minimum axial moments of inertia. Fig. 1 shows the dependency between the two reference systems. The principal angle of inertia ap is conventionally defined as the angle made by the arbitrary Z-direction with the principal maximum axis of e 2 . All vertex coordinates are translated and inertia X rotated to the principal axes of inertia.
1 Notation denotes a measure related to the principal system of inertia.
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
393
Express the warping functions as the double infinite power-series 3 wk ð~x2 ; ~x3 Þ ¼
1 X 1 X i¼0
ðkÞ cij ~xi2~xj3
ð2Þ
j¼0
which are truncated in reality to a reasonable high order n, with i þ j 6 n. In vector notation the components of the warping function are defined as wk ð~x2 ; ~x3 Þ ¼ N cðkÞ
Fig. 1. Cross-section subject to shear loads and torque.
The design sensitivities of all axial properties explicitly known in terms of the vertex coordinates, can be obtained via analytical differentiation, see Apostol and Santos [2]. The axial stress is expressed in terms of a constant stretching force, and two linear bending moments, and depends upon the coordinates (~x2 ; ~x3 ) of the cross-section point under analysis. The design sensitivities of the axial stress also need their corresponding velocity field (~x02b ; ~x03b ). 2 4. Shear properties
A significant simplification in the formulation proposed by Kosmatka [12] is achieved introducing two vectors of size m, namely g, l, containing the corresponding exponents of ~x2 , ~x3 in N n N ¼ ~xg21 ~xl3 1 ; ~xg22 ~xl3 2 ; ~xg23 ~xl3 3 ; ~xg24 ~xl34 ; ~xg25 ~xl35 ; . . . ; ~xg2m 2 ~xl3 m 2 ; o ~xg2m 1 ~xl3m 1 ; ~xg2m ~xl3 m : ð4Þ For the shear properties under investigation, only the first derivatives of the shape functions with respect to coordinates are required oNi ¼ gi~x2gi 1~xl3i ; o~x2
ð1Þ
ð2Þ
s1j ¼ s1j h þ s1j
4.1. Power-series method
s12
wð~x2 ; ~x3 Þ ¼ w1 ð~x2 ; ~x3 Þh þ w2 ð~x2 ; ~x3 Þ
ð1Þ
i ¼ 1; . . . ; m
ð5Þ
P2 ð3Þ P3 þ s1j ; EI3 EI2
j ¼ 2; 3
ð6Þ
where ð1Þ
ow1 ow1 ð1Þ
~x3 ; s13 ¼ G þ ~x2 o~x2 o~x3 ow2 m 2 ow2 ð2Þ s13 ¼ G
~x2 ~x23 ;
m~x2~x3 ¼G o~x2 2 o~x3 ow3 ow3 m 2 ð3Þ
m~x2~x3 ; s13 ¼ G
~x3 ~x22 ¼G o~x2 o~x3 2
s12 ¼ G ð2Þ
s12
ð3Þ
P2 P3 þ w3 ð~x2 ; ~x3 Þ EI3 EI2
oNi ¼ li~xg2i ~xl3 i 1 ; o~x3
Shear stress formulas may be expressed in terms of the three components of the warping function, corresponding to the torsion and the two transverse shear loads, Kosmatka [12],
The calculation of the warping function w via powerseries leads to the shear properties of the cross-section, namely the torsion moment of inertia I1 , the cross-section twist rate h, and the shear center ð~xs2 ; ~xs3 Þ. The cross-section in-plane shear stresses s12 , s13 acting along the principal directions may be determined next. Finally, the design sensitivities of the shear properties, required for design optimization, are obtained by analytical differentiation.
The warping function describes the warping out of the cross-section plane due to transverse loads P2 , P3 applied at the centroid and along the principal directions, and to a torque M1 applied along the beam axis, see Fig. 1. Define the warping function as, Kosmatka [12],
ð3Þ
ð7Þ The three sets of unknown warping coefficients deðkÞ fined in Eq. (3), ci ; i ¼ 1; . . . ; m, can be determined solving the linear system resultant from the virtual work principle, dU ¼ dW . According to Kosmatka [12], after replacing the strains and stresses by their values, into the virtual work principle, the coefficients cðkÞ are obtained from the following linear system of equations
where the twist rate h is a function of M1 , P2 , P3 . 2
Design sensitivities are written as ð Þ0b .
3 ð Þk , or ð ÞðkÞ , k ¼ 1, 2, 3, refers to one component of the warping function.
394
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
KcðkÞ ¼ FðkÞ
ð8Þ
ðkÞ ðkÞ FðkÞ FðkÞ w Fc ¼ ðFw Fc ÞQ
ð9Þ
the stiffness matrix defined as T T ) Z ( oN oN oN oN K ¼ GL þ dA o~x2 o~x2 o~x3 o~x3 A and the force matrices given by Z
Z Fw ¼ EL 0; ~x2 NT dA; ~x3 NT dA
Fc ¼ GL
Z (
2
A
~x2 4
~x3
ð10Þ
ð11Þ
A
oN o~x3
T T ) oN ; o~x2
3
m~x2~x3
2m ð~x23 ~x22 Þ 5dA
2m ð~x22 ~x23 Þ
m~x2~x3
ð12Þ
with Qð1Þ ¼ f1; 0; 0gT ;
nv Z X F2 d~x3 þ F3 d~x2 ¼ F2 d~x3 þ F3 d~x2
C
with the load vector
A
I
Qð2Þ ¼ f0; 1; 0gT ;
Qð3Þ ¼ f0; 0; 1gT ð13Þ
i¼1
ð15Þ
Ci
For a comparison between the surface and contour integration approaches, consider a polynomial of order n ¼ 11. Whereas the exact integration on the contour requires 11 Gauss points along a line, the associated triangle made of the two line end vertices, and the origin, must be integrated using 79 Gauss points for the surface integral. Savings of up to 85% in run-time follow, if contour integrals are used. Regardless of the approach used to carry-out the required integrations, the stiffness matrix and the load vector have to be created by assembling the corresponding values of the base elements, triangles if surface integrals are used, straight lines if contour integrals are used. In the classical finite element method, one finite element depends on a subset of degrees of freedom, and therefore the elementary stiffness matrix adds terms only at particular locations inside the overall stiffness matrix. In the power-series method, any elementary stiffness matrix depends on all warping coefficients cðkÞ , and thus contributes to all locations of the global stiffness matrix. The assembly process degenerates here into a simple fullmatrix summation over all Gauss points, and no sparsity can be exploited.
4.2. Numerical integration on the surface
4.4. Parametric form of a straight line
The integration of the stiffness matrix and force vector is performed on the cross-section domain. The cross-section is viewed as a sum of triangles defined by the corresponding centroid and two adjacent vertices. For complex shapes, with a larger number of vertices, higher order polynomials must be used. Hence, numerical integration on the cross-section domain requires to build the stiffness matrix as proposed by Kosmatka [12]. An effective alternative for the integration of the stiffness matrix and force vector is here proposed. The domain integrals can be transformed into contour integrals, with significant gains in computational time, and allowing for simplified computation of the design velocities required for cross-section shape design sensitivity analysis.
The parametric form of a straight line defined in the e 2GX e 3 system is X
4.3. Numerical integration on the contour Green’s formula states that given two functions F2 and F3 , the following transformation holds I Z oF2 oF3
dA ¼ ð14Þ F2 d~x3 þ F3 d~x2 o~x2 o~x3 A C If the contour C is a polygon with nv vertices, and nv straight sides, the contour integral can be written as
~x2 an ¼ a0 n þ a1 ) d~x2 ¼ a0 dn ~x3 bn ¼ b0 n þ b1 ) d~x3 ¼ b0 dn
ð16Þ
where a0 , a1 , b0 , b1 are constant coefficients, and
1 6 n 6 1 is the parametric variable along the line. Auxiliary variables (an ; bn ) are shortcuts for the longer formulas which follow, and they simply express the coordinates at a certain location n along the line. Considering two adjacent vertices of coordinates ð~xi2 ; ~xi3 Þ and ð~xðiþ1Þ2 ; ~xðiþ1Þ3 Þ, i ¼ 1; . . . ; nv , coefficients ~xðiþ1Þ2 ~xi2 ; 2 ~xðiþ1Þ3 ~xi3 b0 ¼ ; 2
a0 ¼
~xðiþ1Þ2 þ ~xi2 ; 2 ~xðiþ1Þ3 þ ~xi3 b1 ¼ 2
a1 ¼
are determined for each segment of the contour C. For i ¼ nv , the last and the first vertices of the closed contour determine the parametric coefficients. In order to keep formulas easier to read, no script on coefficients a, b indicates a contour segment. The summation introduced in Eq. (15) is implicitly assumed. Design derivatives can be computed straightforward if the velocities ~x02b , ~x03b are known. Hence,
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
a00b ¼ ~x02b
b00b ¼ ~x03b
~x0ðiþ1Þ2b ~x0i2b a0nb
2 ¼ a00b n þ a01b
~x0ðiþ1Þ3b ~x0i3b b0nb
;
;
2 ¼ b00b n þ b01b
a01b ¼
~x0ðiþ1Þ2b þ ~x0i2b 2
;
ð17Þ
h¼
1 P2 P3 M1 a2
a3 EI3 EI2 a1 ðkÞ
b01b ¼
~x0ðiþ1Þ3b þ ~x0i3b 2
;
ð18Þ
395
ð24Þ
ðkÞ
with s12 , s13 given by Eqs. (6) and (7). Replacing the shear stresses into Eq. (23), for k ¼ 1, follows Z Z Z ow1 ow1 2 2 ~x2 a1 ¼ G
~x3 dA þ ~x3 dA þ ~x2 dA o~x3 o~x2 A A A ð25Þ
4.5. Solving the linear warping system An element of the symmetric stiffness matrix Kij , i ¼ 1; . . . ; m, j ¼ i; . . . ; m, introduced in Eq. (10), can be written as Z g þg 2 l þl g þg l þl 2 dA Kij ¼ GL gi gj~x2i j ~x3i j þ li lj~x2i j ~x3 i j ZA oF2 oF3
dA ¼ GL o~x2 o~x3 A ð19Þ Identifying now the two functions of the Green’s formula it follows gi gj g þg 1 l þl ~x i j ~x3 i j ; gi þ gj 1 2 li lj g þg l þl 1 ~x i j ~x3i j F3 ¼
li þ lj 1 2 F2 ¼
ð20Þ
Green’s formula (14), and the parametric form of the line (16) assure that Z 1 Kij ¼ GL ðF2 b0 þ F3 a0 Þdn ð21Þ
1
or Z
ð22Þ
4.6. Computation of shear properties
A
and the twist rate as
5. Design sensitivity analysis The parametric definition of the cross-section assures the connection between the vertex coordinates and the selected design parameters. The design sensitivities P0ab of the axial properties ~2 ; x ~3 ; r11 g I2 ; I3 ; ap ; x
The linear symmetric warping system (8) can be solved now for the three sets of coefficients cðkÞ . Therefore, a two-step based solver, using forward elimination and backward substitution, is more suitable than an iterative solver.
Define the warping constants, Kosmatka [12], Z ðkÞ ðkÞ ~x2 s13 ~x3 s12 dA ak ¼
where the axial moments of inertia are identified with the last two integral terms in Eq. (25). Eq. (26) clearly shows the error involved in approximating I1 ffi I2 þ I3 , since the first integral term on the right hand side is neglected. In fact, this expression is correct only for circular cross-sections. Warping coefficients a2 , a3 can be used to locate the shear center a3 a2 ð~xs2 ; ~xs3 Þ ¼ ð27Þ ;
EI2 EI3
Pa ¼ fzT ; yT ; AT ; SzO ; SyO ; A; IzzO ; IyzO ; IyyO ; IzzG ; IyzG ; IyyG ; zG ; yG ;
1
gi gj g þg 1 l þl a i j bn i j b0 gi þ gj 1 n
1 ! li lj gi þgj li þlj 1 a bn a0 dn
li þ lj 1 n
Kij ¼ GL
The warping constant a1 is nothing but the torsional rigidity GI1 . Hence, the torsion moment of inertia I1 becomes Z ow ow ~x2 1 ~x3 1 dA þ I2 þ I3 ð26Þ I1 ¼ o~x3 o~x2 A
ð23Þ
were presented in Apostol and Santos [2]. In this section the design sensitivities P0sb of the shear properties below are obtained, Ps ¼ fI1 ; h; ~xs2 ; ~xs3 ; s12 ; s13 g: Computation of the shear properties requires the solution of a linear system of equations. The analogy with the finite element method allows a parallel between the design sensitivity methods developed for the linear finite element method and for the power-series method. Whereas the displacement field is the unknown solution in the finite element method, the present method uses the warping coefficients cðkÞ as the solution of the governing ðkÞ0 system of equations. Hence, the design sensitivities cb are the basis for the design sensitivities of the shear properties.
396
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
Since the entire design sensitivity field is required, the direct differentiation method is here preferred against the adjoint variable method. Differentiating wk ¼ N cðkÞ with respect to design, ðkÞ0
w0kb ¼ N0b cðkÞ þ Ncb
ð28Þ
Differentiating the first derivative of wk with respect to design follows ow0kb oN0b ðkÞ oN ðkÞ0 ¼ c þ c ; o~xj b o~xj o~xj
j ¼ 2; 3
ð29Þ
Taking advantage of the closed-form formulas (5), the above equation becomes m 0 X ow0kb ðkÞ g 1 l g 1 l ðkÞ0 ¼ gi~x2i ~x3i ci þ gi~x2i ~x3 i cib b o~x2 i¼1 ð30Þ m 0 X 0 owkb ðkÞ gi li 1 gi li 1 ðkÞ0 ¼ li~x2 ~x3 ci þ li ~x2 ~x3 cib b o~x3 i¼1 or arranged in compact form m nh i X ow0kb g 2 l g 1 l 1 ¼ gi ðgi 1Þ~x2i ~x02b ~x3i þ gi li ~x2i ~x3i ~x03b o~x2 i¼1 o ðkÞ g 1 l 0ðkÞ ci þ gi ~x2i ~x3 i cib m nh i X ow0kb g 1 l 1 g l 1 ¼ li gi ~x2i ~x02b ~x3i þ li ðli 1Þ~x2i ~x3 i ~x03b o~x3 i¼1 o ðkÞ ðkÞ0 ci þ li ~xg2i ~xl3 i 1 cib
ð31Þ These formulas, and Eqs. (6)–(27) show that any design sensitivity of a shear property Ps can be later computed if additionally to the previous computations, the veloc~03b Þ, and the warping derivatives cðkÞ0 ity field ð~ x02b ; x are b available. Recall that the warping coefficients cðkÞ are computed from the linear system KcðkÞ ¼ FðkÞ , by differentiating it with respect to design ðkÞ0
K0b cðkÞ þ Kcb
ðkÞ0
¼ Fb
ð32Þ
Moving the first term from the left side to the right side, results in ðkÞ0
Kcb
ðkÞ0
¼ Fb K0b cðkÞ
ð33Þ
This system has the same stiffness matrix as the original, and only backward substitutions are required to comðkÞ0 pute cb . Derivatives of the stiffness matrix, load vectors, and warping constants ak , require integration on the contour, i.e. the velocity field in any Gauss point along the contour. This information is obtained according to Section 4.4.
If surface integration is used instead of contour integration, as proposed by Kosmatka, no velocities would be available in the interior of the contour. This disadvantage is even more severe than the higher runtime used by the surface integration, and ultimately influenced the use of contour integration. The algorithm used for the design sensitivity calculation of the shear properties is as follows: • Compute the stiffness matrix K, and the loads FðkÞ . The stiffness matrix is scaled and decomposed after the Gauss forward elimination. • Compute design sensitivities of the stiffness matrix ðkÞ0 K0b , and loads Fb . They are used to build the pseudoloads associated with the right side of Eq. (33). • Solve 3 np times the system (33), where np is the number of design variables, through backward substitutions on the unchanged, decomposed stiffness ðkÞ0 matrix, and get the coefficients cb . • Integrate the design sensitivities of the warping constants a0kb . • Compute the design sensitivities of the shear proper0 ties I1b , h0b , ð~x0 x0 s2b ; ~ s3b Þ. • Compute the design sensitivities of the shear stresses ðkÞ0 sb , for any desired point along the contour. It is worth noticing that the theoretical results presented in this paper together with those published in Apostol and Santos [2] enable the computation of the axial and shear design sensitivity contributions of any constraint that can be written as a function of the axial and shear properties included in both vectors Pa and Ps .
6. Optimization of a thin-walled cross-section subject to shear stress constraints This section presents the optimization of a thin-walled cross-section, with respect to shear stress constraints, for two different parametrizations, direct vertex parametrization, and indirect vertex parametrization using NURBS. 6.1. Direct vertex parametrization A thin-walled cross-section is parametrized using direct non-linear relationships between vertex coordinates and design variables, see Fig. 2. The initial crosssection shape is described as a regular 20-sided polygon, with inner and outer vertices expressed in terms of 18 geometrical dimensions and the wall thickness. The cross-section is made of steel with a Young’s modulus of 2.1e þ 05 MPa, and the Poisson’s ratio is 0.3. The loading consists of a torque of Mx ¼ 1e þ 05 N mm, and two shear forces Py ¼ Pz ¼ 1e þ 04 N.
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
397
Table 1 Optimization results of a thin-walled cross-section
Fig. 2. Definition and loading of a thin-walled cross-section.
The optimization problem to be solved is the minimization of the cross-section area subject to shear stress constraints along the cross-section boundary. The shear stresses are defined at nine equidistant points within each boundary segment. A total of 234 stress performances are included in the optimization model. The shear stress constraints are defined as s ¼ ðs2xy þ s2xz Þ1=2 , where smax ¼ 100 MPa. The design parametrization requires the wall thickness to be constant along the wall cross-section contour. The horizontal vertex coordinates must be kept constant, and the vertical vertex coordinates are allowed to vary linearly with respect to the design variables, for the inner vertices, and non-linearly for the outer vertices. The analytical expressions for the vertices may be easily obtained, and therefore they are not presented here. Whenever a new design is investigated, the vertex coordinates must be recomputed as well as the nodal velocity field. The direct differentiation method described in Section 5 is used for nodal velocity computations. The output of these computations is used for the cross-section shear related geometric and performance computations. The velocity field within the contour is interpolated linearly from the nodal velocities. The optimization problem is solved with a sequential quadratic algorithm, see Apostol [18] for more details. Table 1 shows the initial and the final results after 20 iterations. The initial design is feasible since all shear stresses are smaller than 100 MPa. The 15% structural reserve allows the area to decrease from 453.18 to 351.58 mm2 . At the 20th iteration, 16 out of 234 constraints become active, i.e. the corresponding shear stresses are close to 100 MPa. Fig. 3 shows the first and the last shapes, for comparison. The optimizer did not converge yet, but no significant changes occurred after the 20th iteration. Fig. 4 presents the objective function iteration history. However, the optimum shape is not likely to be suitable
DV
Initial (mm)
Final (mm)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
9.270 17.633 24.270 28.531 30.000 28.531 24.270 17.633 9.270
9.270
17.633
24.270
28.531
30.000
28.531
24.270
17.633
9.270
10.057 19.675 30.381 35.018 34.507 29.810 27.302 19.586 15.365
5.303
12.582
32.162
36.486
41.115
39.878
24.541
17.982
5.515
Properties Thickness t (mm) Area A (mm2 ) Tors. mom. I1 (mm4 ) Max. s (MPa)
2.000 453.18 3.842eþ05 85.267
1.395 351.58 3.240eþ05 100.020
for manufacturing since it has a very irregular contour. A new parametrization, typical of many CAD systems, is used next to overcome the above problems and verify the appropriateness of the proposed formulation for indirect vertex parametrization. This example clearly shows the application of the theoretical developments presented in Sections 4 and 5. Despite the highly non-linear relationships between the shear stress constraints and the design variables, still, regular optimization paths and convergence were achieved, demonstrating that the power-series method is a practical tool for structural optimization. 6.2. Indirect vertex parametrization using NURBS The approach of direct vertex definition with shape variables was earlier intensively used. Nowadays, interpolation curves are preferred in modern shape optimization, having at least two advantages: they are controlled by fewer design variables, and they assure the smoothness required for manufacturing. In particular NURBS were used by Schramm and Pilkey [8,9] for design parametrization and integration with structural shape optimization. The thin-walled structure illustrated in Fig. 2 is now parametrized with two NURBS, see Fig. 5. The inner vertices are defined by the two independent sets of control points: p1 , p2 , p3 , p4 , p5 , p5 , and p5 , p6 , p7 , p8 , p1 , p1 , with coordinates ðzj ; yj Þ given under the column ‘Initial’
398
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
Fig. 3. Initial and final shapes of a thin-walled cross-section.
in Table 2. Remark that the last vertex has to be entered twice. This is a characteristic of B-splines. The polynoTable 2 Optimization results using NURBS
Fig. 4. Optimization history of a thin-walled cross-section.
Fig. 5. Definition of the thin-walled cross-section using NURBS.
DV
Initial (mm)
Final (mm)
z1 y1 w1
30.00 0.00 1.00
30.00a 0.00a 1.00a
z2 y2 w2
30.00 19.60 1.00
32.00a 29.68 1.92
z3 y3 w3
0.00 34.00 1.00
4.49 37.46 1.25
z4 y4 w4
30.00 17.60 1.00
31.96 10.09 1.96
z5 y5 w5
30.00 0.00 1.00
30.00a 0.00a 1.00a
z6 y6 w6
30.00
17.60 1.00
32.00a
40.00a 1.19
z7 y7 w7
0.00
34.00 1.00
20.00a
40.00a 1.96
z8 y8 w8
30.00
17.60 1.00
20.00a
5.00a 1.99
2.00 456.40 3.977eþ05
1.21 304.81 2.154eþ05
83.29
100.03
Properties Thickness t (mm) Area A (mm2 ) Tors. mom. I1 (mm4 ) Max. s (MPa) a
Design variable bounds imposed for constructive reasons.
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
399
Fig. 6. Initial and final shapes of inner vertices using NURBS.
mial order is fixed to n ¼ 3, whereas each B-spline is transformed into a polygon with m ¼ 31 vertices. Fig. 6 left shows that, due to the non-uniform knot sequence, the intervals between vertices are not equal. The shape is approximately circular. However, the initial cross-section area defined through NURBS differs by less than 1% relatively to the previous optimization with direct vertex parametrization. The outer vertices are analytically computed such that the thickness t remains constant over the crosssection. Each horizontal flange does not change its dimensions, except for the thickness. Control points p1 and p5 are kept constant during the optimization process. The other six control points define 18 design variables (yj ; zj ; wj ), where wj are the weighting factors. The two curves act independently, and therefore, sharp corners could appear at the connection between the upper and the lower curves. The bounds on the adjacent control points may control this undesirable effect. The thickness t is the 19th design variable. The optimization problem to be solved is the same as that of Section 6.1. It is required to minimize the cross-
section area subject to shear stress constraints, along the cross-section boundary. The shear stresses are defined in this case as three shear stresses at equally spaced points along each segment, on the outer contour. A total of 198 constraints of the type s 6 smax , with smax ¼ 100 MPa, are specified. The convergence of the optimization process was rather slow, requiring 40 iterations. The final results of the model are given under column ‘Final’ in Table 2. Some of the design variables reached their bounds. This fact shows that slightly improved shapes could be found if the bounds were relaxed. However, in practical optimization the objective is to find meaningful shapes rather than a mathematical optimum. For comparison, Figs. 6 and 7 present the initial and the final shapes of the inner contour, and of the entire cross-section, accordingly. Fig. 8 shows the optimization history of the crosssection area. The jump in the objective function value observed at the 21st iteration was caused by the stop of the optimizer after the 20th iteration. After restarting the optimizer, the second order Hessian was lost, and a poorer solution was reached at the 21st iteration.
Fig. 7. Initial and final shapes using NURBS.
400
V. Apostol et al. / Computers and Structures 80 (2002) 391–401
Fig. 8. Optimization history using NURBS.
Although the optimizer did not converge even after 40 iterations, since the variations in design were at least four orders of magnitude smaller than the units in millimeters, the process was stopped. Some comments are due to this point, regarding the optimization of the beam cross-section. As already verified in classical shape optimization, the parametrization using NURBS leads to a final smooth shape, whereas the parametrization using piecewise linear segments leads to an irregular final design, see Fig. 3. Also, the smoother parametrization associated with the NURBS allows convergence to a significantly better design than the piecewise parametrization of the boundary, leading to final cross-section areas of 304.81, and 351.58 mm2 , respectively. The optimal shape associated with the NURBS boundary parametrization illustrates better the optimal material distribution required to reach the highest stiffness under the shear constraints, by aligning the material preferably along two quasi straight-sided segments, almost parallel to each other, and to the 45 oriented shear force. A final comment on the selection of the NURBS weighting factors is appropriate. By defining these weighting factors as design variables it was possible to search for the optimum within a large design space. In Fig. 7a rounded curve at the right-upper side of the cross-section is associated with a weighting factor w3 ¼ 1:25, whereas a sharp corner at the lower-left side corresponds to a weighting factor w7 ¼ 1:96.
7. Conclusions A general approach for detailed analysis and design optimization of cross-sections of arbitrary shape of truss/beam built up structures is proposed. The approach allows for both direct and indirect parametric
definition of an arbitrary cross-section, and may be easily integrated with the parametrization supported by state-of-the-art CAD systems. The truss/beam components are subjected to axial and transversal forces, and to bending and torsion moments. Both shear properties and stresses are considered in both the analysis and the optimization process. The power-series method is used for the solution of the shear problem, and boundary integration is proposed for carrying-out the integrations required for the solution of the warping equations. A method for computing the design sensitivities of the shear properties and stresses with respect to crosssection design parameters is developed. This method uses direct differentiation of the warping equations and only requires backward substitutions on the triangular stiffness matrix, making it a quite attractive alternative for structural shape optimization of truss/beam crosssections of built up structures. A numerical example illustrates the applicability of the method to design optimization of a real-life crosssection.
References [1] Apostol V, Santos JLT. Design sensitivity analysis and optimization of truss/beam structures with arbitrary crosssections. In: Olhoff N, Rozvany G, editors. WCSMO-1 Proceedings of the First World Congress of Structural and Multidisciplinary Optimization, 28 May–2 June, Goslar, Germany. Oxford: Pergamon Press; 1995. p. 353–8. [2] Apostol V, Santos JLT. Sensitivity analysis and optimization of truss/team components of arbitrary cross-section. I––Axial stresses. Comput Struct 1996;58(4):727–37. [3] Timoshenko SP, Goodier JN. Theory of elasticity. 3rd ed. New York: McGraw-Hill; 1982. [4] Banichuk NV. Optimization of elastic bars in torsion. Int J Solids Struct 1976;12:275–86.
V. Apostol et al. / Computers and Structures 80 (2002) 391–401 [5] Pilkey WD, Liu Y. Field theory: A two-dimensional case for not using finite or boundary elements. Finite Elem Anal Des 1993;13(2/3):127–36. [6] Schramm U, Pilkey WD. Structural shape optimization for torsion problem using direct integration and B-splines. Comput Meth Appl Mech Eng 1993;107:251–68. [7] Schramm U, Kitis L, Kang W, Pilkey WD. On the shear deformation coefficient in beam theory. Finite Elem Anal Des 1994;16:141–62. [8] Schramm U, Pilkey WD. The coupling of geometric descriptions and finite elements using NURBs––a study in shape optimization. Finite Elem Anal Des 1993;15: 11–34. [9] Schramm U, Pilkey WD. Optimal shape design for thin walled beam cross-sections. Int J Numer Meth Eng 1994;37(23):4039–58. [10] Mindlin RD. Solution of St. Venant’s torsion problem by power series. In: Solid structures, vol. 11. Oxford: Pergamon Press; 1975. p. 321–8. [11] Kosmatka JB. On the behavior of pretwisted beams with irregular cross-sections. ASME––J Appl Mech 1992;59(1): 146–52.
401
[12] Kosmatka JB. Flexure–torsion behavior of prismatic beams, part I: Section properties via power series. AIAA J 1993;31(1):170–9. [13] Dunavant DA. High degree efficient symmetrical Gaussian quadrature rules for the triangle. Int J Numer Meth Eng 1985;21:1129–48. [14] Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical recipes. Cambridge: Cambridge University Press; 1989. [15] Langelaan JW, Livne E. Analytic sensitivities and design oriented structural analysis for airplane fuselage shape synthesis. Comput Struct 1997;62(3):505–19. [16] Braibant V, Fleury C. Shape optimal design using Bsplines. Comput Meth Appl Mech Eng 1984;44:247–67. [17] Lindby T, Santos JLT. 2-D and 3-D shape optimization using mesh velocities to integrate analytical sensitivities with associative CAD. Struct Optimization 1997;13(4): 213–22. [18] Apostol V. Design sensitivity analysis and optimization of geometrically nonlinear structures in the pre- and postbuckling domains. PhD thesis, Instituto Superior Tecnico Lisbon, Portugal, 1998.