Primary Response Assessment Method for Concept Design of Monotonous Thin-Walled Structures

The concept design methodology for monotonous, tapered, thin-walled structures (wing/fuselage/ship/bridge) is presented. The problem solution is based on the OCTOPUS program [1, 2]. It contains: (A) response and feasibility analysis modules (FIN-CREST), (B) decision making-synthesis modules (DeMak) and (C) interaction/visualization programs (MAESTRO MM/MG and DeVIEW) that irerate in the design cycle. The modules are summarized in Table 1 as modules 1a–8c.


Introduction
The concept design methodology for monotonous, tapered, thin-walled structures (wing/fuselage/ship/bridge) is presented.The problem solution is based on the OCTOPUS program [1,2].It contains: (A) response and feasibility analysis modules (FIN-CREST), (B) decision making-synthesis modules (DeMak) and (C) interaction/visualization programs (MAESTRO MM/MG and DeVIEW) that irerate in the design cycle.The modules are summarized in Table 1 as modules 1a-8c.
(A) The analytical (CREST) modules and methods are fully described in [2].Module-1a INDAT is used for data generation, combined with the MAESTRO FEM MODELER [8].Module-2 LOAD is used for design load generation.Module-1b MIND is used for determining the minimal scantlings based on prescribed rules (LR, DnV, ABS, etc.) Module-3a LTOR is used for direct calculation of the primary strength (shear flow and corrected stresses in bending and warping torsion).They are calculated using an original extended beam theory [5,6].Module-3b TOKV is used for the transverse strength calculation (newly developed 8-node stiffened panel macro-elements are used for modeling the transverse structural response).Module-4 is the PANEL library of structural serviceability and ultimate strength criteria for the structural adequacy calculation [4], using the response fields generated in modules 3a and 3b.Module-8a is VB-SHELL for the designer-model interaction.Module-8b MAESTRO GRAPHIC [8] is used for presenting the model (loading, response, adequacy, etc.) (B) Synthesis (DeMak) modules and methods are documented in [1,3].Local variables for substructures (s = 1, …, NS), are denoted x s = {x i } s = {t plating , n stiffeners , h web …} s .Substructure areas X = {X s }are intermediate (global) variables, where X s = X s (x s ).Project k is defined as P k = {x 1 , …, x NS , x fixed } k .Design criteria (attributes, objectives, constraints) are formulated as a library of mathematical func-tions/procedures for driving the optimization process or feasibility check.OCTOPUS metamodeling of failure surfaces is based on the most unsatisfied constraint from each local problem.They are added to the set of global constraints.The value function for a global level is a multicriterion combination of normalized attribute functions.The solution strategy involves generation of designs, using (a) a Random Number Generator in the first cycles of design space exploration, and (b) Fractional Factorial Experiments for subsequent cycles.Coordination is performed by modifying v(x s ) respective to its divergence from globally optimal substructure area X s .Special provisions: -generation of promising designs using 27 designs obtained from Orthogonal Array L27.
-extensive usage of tables of optimized profiles to speed up the generation process.
Modules 6 and 7c (GAZ) are used for calculating the sensitivity of the structural response with respect to the design variables, based on the global strength module FIN/LTOR [7].Module-7a GLO is used for global level MODM optimization (level 1) of the cross-section.Module-7b LOC is used for local coordinated MADM decision making via sequential application of stochastic search methods and Theory of experiments.
(C) Module-8c DeVIEW is used for the designer-model interaction with graphic presentation of designs in the design and attribute spaces.The stratified distances from the target (or the ideal) design, calculated by L p metric are used as a means of visualizing the multidimensional space of the design attributes and/or free variables.Visualization is the most powerful tool for the designer's understanding of the Decision Support Problem.It generates expert knowledge about the problem for all participants involved, and helps the designer to identify advantageous combinations of variables, feasible options and clusters of non-dominated designs (Pareto frontier, thus enabling realistic decision support to the head and structural designer).

Modeling philosophy for primary response in concept design
Classical FE modeling, giving good insight into stresses and deformations, is not capable of giving the efficient and fast answers regarding feasibility criteria (buckling, fatigue, yield) required by the Rules.However, structural feasibility and compliance with the Rule requirements are of primary interest, not stresses or deformations.
Most of the local failure criteria, e.g.various buckling failure modes of stiffened panels, require specified force and displacement boundary conditions.They are available only if logical structural parts, such as complete stiffened panels between girders and frames, are modeled (macro-elements).
For the concept design structural evaluation of the primary response (longitudinal strength, torsional strength), the beam idealization of a wing/ship/bridge is often used.A primary strength calculation provides the dominant response field (Demand) for design feasibility assessment.The evaluation is based on extended beam theory, which needs cross-sectional characteristics.These are obtained using analytical methods, which can be very complicated for real combinations of open and closed cross-sections.Application of energy based numerical methods gives an opportunity for an alternative approach to the given problems.The method is based on decomposing a cross-section into the line finite elements between nodes i and j with coordinates (y i , z i ), (y j , z j ); element thickness t e ; material characteristics (Young's modulus  Using the FEM approach, a procedure is developed for calculating the set of cross-sectional geometric and stiffness characteristics at position x denoted G x with the following elements: l Cross-section area A S Center of gravity Y CG , Z CG , S Shear/torsion center Y CT , Z CT .l Moments of inertia with respect to the center of gravity: I Y , I Z , I yz , I p ; principal: I 1 , I 2 , f 0 -angle of axis-1 w.r.t.Z-axis, l Horizontal and vertical bending: S Flexural stiffness EI Z , EI Y S Shear stiffness GA V , GA H , l Cross-section axial stiffness EA S Torsional stiffness GI T S Warping stiffness EI W . Standard stiffness matrices (alternatively with geometrical nonlinearity [4]), for axial (k u ) and flexural (k v , k w ) and torsional response modes are given as functions of the geometric set G x : k w is obtained similarily.Stiffness matrices for free torsion (k T ) and restrained warping (k w ) read:   ters the element axial, bending and torsion parameter distributions, based on the applied shape functions, can be derived (e.g.q(x), q ,x (x), q ,xx (x), q ,xxx (x) for torsion).The key element for calculation of the response of the complex thinwalled structure is therefore determination of elements of G x .A simple and elegant FEM procedure for such a calculation is presented in the sequel.
3 Calculation of response for a transverse strip with a complex cross section The shear flow and geometrical characteristics of the cross section in bending and torsion is usually calculated using analytical methods.Such calculations become rather complicated for multiple-connected cross section graphs with a combination of open and closed (cell) contours.Application of numerical methods based on the energy approach offers an elegant alternative.The procedure is based on section decomposition into finite elements, as first introduced by Herman and Kawai.In the sequel, the method of calculation as in [5,6] is presented.It has been successfully used in practical calculations since its development for [10].The simplest decomposition of thin-walled cross-section (symmetric or not) into line finite elements (segments) is shown in Figs. 1 and 2. These elements form stiffened panel macro-elements for the feasibility evaluation.
The methodology is based on applying the principle of minimum total potential energy (P) with respect to parameters which define the displacement fields of the structure.The primary displacement field (following classical beam theory) is defined via displacements and rotations of the cross section as a whole.Secondary displacement field u u 2 ( , , ) ( , , ) x y z x y z º represents warping (deplanation) of the cross section.For piecewise-linear FEM idealization of the cross-section, divided into n elements, with shape functions N in the element coordinate system (x, s), the warping field reads: Element strain and stress fields ò and s are obtained from the strain-displacement and stress-strain relations: where RS is the prescribed shear efficiency.

Cross-sectional shear stress distribution due to bending
In the case of bending, the net external load (due to bending moments M(x+Dx) and M(x)) is the normal stresses: ( ) is distance from the point to N. A. The load vector for a nonsymmetrical cross-section in, e.g., bending about the z axis reads: For bending around the Y and Z axes, the matrix relations The approximate value of normal stress for simultaneous bending about axes y and z for node I reads: ) .é ë ê ù û ú

Calculation of warping and primary shear stresses due to pure torsion
A transverse strip of a thin-walled beam of length Dx is subjected to torsional loading.The displacement field of the middle line of thin walled elements can be expressed using the warping function u s t ( ) =0 , rotation v x s t ( ) =0 around the centre of twist, twist rate (q x,x ) and angle (q x ) of the twist reads: where d T is the normal distance from the element to the center of torsion.The strain (with e s » 0) and stress fields read: , The total potential energy of a section is given by the standard expression: After summation of all elements and transformation of local element displacements u N u e e

= ×
T and loads F e into global displacements u and loads F we get: Minimization of total potential energy leads to two sets of equations: (1) d q P ¢ = 0 (1D beam torsion) and (2) dP u = 0 (2D cross-section warping).
A second set of equations, d d P u u U = = ® = 0 Ku F , enables determination of the unit warping field.
The primary shear stresses on the elements which are parts of closed contours (cc) and open sections (os) can now be calculated as functions of 1D twist rates q ,x (x) (to be obtained from the first relation for 1D beam torsion):

Calculation of torsional and warping stiffness of thin-walled structures
To solve the equation for 1D beam free torsion, the torsional stiffness of elements which are parts of the open e o and closed cells e c can now be calculated using the known unit warping field u: Warping stiffness is calculated using the expression: Using GI T and EI W , the matrix K 1D for the 1D beam problem can be formed and relevant parameter distributions q(x), q ,x (x), q ,xx (x), q ,xxx (x) can be determined for use in shear stress calculations.

Normal and secondary shear stresses due to restrained warping
Restrained warping of a thin-walled beam will induce (a) normal stresses in a cross-section and (b) secondary shear stresses which will balance the longitudinally non-uniform distribution of normal stresses.This additional mechanism will influence the strain energy and the work of expression, so an iterative solution may be needed for greater accuracy.
Let u x s u s x x x ( , ) ( ) ( ) , = ×q be the warping field in the cross--section calculated from the case of free torsion.Normal stresses are caused by restraining the warping, and vary along the x axis.They are given by: ( ) Let u 2 (x, s) be the secondary displacement field containing a displacement correction due to restrained warping.The total potential energy of a transverse strip consists of the internal energy generated from the fields ò The shear stress distribution can be calculated more accurately along the element (similar to the bending case)

Examples
The first example is based on reports from the Advanced Subsonic Technology (AST) program.In the course of this program an experimental model of a composite wing box was made and tested (Fig. 3).[11] gives the loads carried by the hydraulic actuators which simulate the in-flight loading conditions.
The example shows a way of rapidly modeling and calculating the overall response of a similar metal wing box with linear behavior during the early stages of wing structural design.The analyzed wing box with reference to the AST box was shortened to 9.8 meters and modeled with high strength aluminum alloy 7075.The loads are decreased for the shortened wing box.The wing box is modeled by OCTOPUS 1D/2D combination.The accuracy of the method is demonstrated in Figures 5a  and 6.Strake 5 of Fig. 5a is located in the middle of the upper skin.Fig. 6 presents the application to the U-channel and two standard ship structures.It can be seen that the accuracy of the shear stress distribution based on FEM (constant per element) and analytical formulae (continuous line) in examples 6a and 6c is very good, even without parabolic correction.
The verification examples are taken from [5].

Conclusions
A simple and practical method for calculating the primary response of monotonous structures (wings, ships, bridges) has been presented.All cross-section parameters are easily determined for complex stiffened thin-walled structures using a special FEM procedure.It could successfully replace classical, often cumbersome, analytical calculations.The method has been in constant use since 1980, applied to many real structures for concept design (in the OCTOPUS system) or as the generator of the force boundary conditions for partial 3D FEM models.

.
Global stifness matrix K 1D is obtained by the combination of modal stiffnesses corrected for centroid and shear centre position relative to the position of the origin of global C.S.
, the global stiffness matrix K 1D is obtained as the sum of the element stiffness matrices k T k T G L e T e = with ap-propriate node numbering.The system K 1D a = F can now be solved for unknown displacements a which in turn enable determination of element parameters a L .From these parame-

Fig. 2 :
Fig. 2: OCTOPUS/MAESTRO GRAPHIC shear stress fields in wing and bulk-carrier transverse strip energy of the Dx -long transverse strip of the beam, with the cross-section divided into n elements, reads: is the external loading on two cross sections (S 1 and S 2 ) of the strip.Minimization of P leads to the classical FEM matrix relation K 2D u 2D = F 2D (shortened to K u = F).The element stiffness matrix for the proposed linear displacement distribution along the line element (the same for bending and torsion) reads:

2 2 1 . 5
In this case, the sectional characteristics and shear center are easily obtained.The shear/torsion center position reads: C is the normal distance from the centroid to e.The shear stiffness for bending about the Y and Z axes, GA V , GA H reads: Corrected normal stresses due to the influence of shear (shear lag)The normal stress must be corrected for (a) stress arising from a longitudinal change of the warping field and (b) normal stress due to correcting bending moment (M c ), compensating for the loss of cross section equilibrium: 2 and s 2 (based on u 2 ) and the additional work done by the strip axial load p x on the secondary displacements u 2 .If the change of u 2 along strip length Dx is neglected, the total potential energy reads: load Dp x due to restrained warping reads: potential energy of the element, using the same shape functions as before, reads: total potential energy with respect to the unknown displacement field u 2 leads to: is the global stiffness matrix as before, u 2 is the global vector of unknown displacements u u load and the secondary shear stresses (constant on element) read:

Fig. 4 Fig. 3 :
Fig. 3: MD 90 airplane and the test model of the wing box

Table 1 :
Summary of OCTOPUS modules x Dxp x+Dx x qy q(x)