Multidisciplinary Optimization Standardization Approach for Integration and Configurability MOSAIC Project
Task 6 WING–BOX STRUCTURAL DESIGN OPTIMIZATION
Report 6 3D Finite Element Model of DLR-F6 aircraft WingBox Structure, Created in PATRAN and Analyzed in NASTRAN By Mostafa S.A. Elsayed, M.Sc. Ph.D. Candidate Amandeep Sing M.Sc. Student Ramin Sedaghati, Ph.D, P.Eng. Associate Professor Principle Investigator for Task 6 Department of Mechanical and Industrial Engineering Concordia University Sponsor’s Ref. No: CRIAQ 4.1-TASK 6
September 2006
1
3D Finite Element Model of DLR-F6 aircraft WingBox Structure, Created in PATRAN and Analyzed in NASTRAN Abstract The current report is a summary of the work done in stages three and four of TASK 6 of the MOSAIC project. The third stage is completed where a methodology for the design optimization of the spars and spar caps of the wing box is presented. The methodology is based on the incomplete diagonal tension theory where the spars are considered in the diagonal tension state of stress. A set of constraints are applied with an objective function of mass minimization which generated an acceptable results. The data obtained from stage three along with the data of the stiffened panels and the ribs are all employed to generate the 3D finite element model of the wing box structure. MSC.PATRAN is used as a modeler while MSC.NASTRAN is used as analyzer. The generated 3D FEM is validated by testing the performance of the model due to the application of a set of aerodynamic loads representing normal cruising conditions. On the other hand a stick model of the 3D finite element model is generated where the flexibility method is used to evaluate the model stiffness properties. Also, a set of empirical formulas generated by the Bombardier aerospace are used to generate the stiffness properties of the ideal stick model of similar aircraft wing-box. The empirical stick model performance is compared with the performance of the 3D FEM and its stick model which showed a great agreement. The comparison showed that the design methodology followed in this project stage is a conservative design where the model generated is stiff model compared with the ideal one. On the other hand the design methodology succeeded in achieving a weight reduction in the wing-box structure as previously explained in previous reports.
Key Words:
Wing-Box, Spars and spar caps, stick models, Diagonal
Tension, Multi Disciplinary Design Optimization (MOD). 2
Section I
Overview
I-1 Description of the Project: The objective of task 6 in the MOSAIC project is to improve the available structural analysis modules in the Bombardier Aerospace and perform a structural design optimization of the wing box by adding an optimization loop around the analysis code. The objective is to design a wing-box more rapidly and automatically. Task 6 is divided into four stages.
Stage I: Optimization of one skin stringer panel: (finished) Stage I explained in details the procedure to optimize one skin-stringer panel consists of one stringer with one stringer spacing (or pitch) of skin in the chord wise direction and the distance between two ribs in the span wise direction. Skin-stringer panels on the upper and lower wing covers are considered. The load acting on the panels is taken to be constant (i.e. same load acting on all panels) which resulted in identical dimensions for all panels. Stage-I provides a methodology to obtain the optimum dimensions for a skin-stringer compression panel with a minimum mass under six constraints namely crippling stress, column buckling, up-bending at center span (compression in skin), downbending at supports (compression in stringer outstanding flange), inter-rivet buckling and beam column eccentricity. It also provides optimum design variables for panels under tensile loading with fatigue life as a design constraint with same objective function (Minimum mass for panel). A panel on the lower wing cover is designed for Damage Tolerance. (For more details refer to report II and III)
Stage II: Load Redistribution: (Finished) Stage II presented the methodology for calculating the actual load experienced by each skin-stringer panel when arranged on the airfoil profile at any span wise section of the wing. The number of stringers required on the upper and lower wing covers is obtained by dividing the width of the wing-box by their corresponding stringer pitch obtained from stage I. These panels are then re3
arranged on the actual airfoil profile at certain span wise section. Each panel now experiences different magnitude of compressive or tensile load depending on its relative location with respect to the centriodal axes of the section. The optimum dimensions for panels on upper and lower wing covers are thus obtained using stage-I optimization program with the new calculated design load which resulted in a different optimum dimensions for each panel according to its location. (For more details refer to report IV).
Stage III: Optimization of the Spars and Spar Caps: (Finished) This stage is an extension to stage II. In this stage the development of the optimization tools to include the spars thickness and web cap dimensions will be considered.
Stage IV: 3D FE Model of the wing box: (Finished) This stage is the subject of the current report. Please read below for details. The DLR-F6 aircraft has been chosen as a practical example to apply the optimization methodology under investigation.
I-2 DLR-F6 Aircraft Geometry and Wing Details: The geometry and load details are taken from DLR-F6 aircraft [5]. The actual wind tunnel model geometry is shown in Figure (1).
4
Fig. (1) DLR-F6 wind tunnel model [5]
y* y
Fig. (2) DLR-F6 wind tunnel model geometry [5] 5
In Figure (2), Axes x, y and z denote the coordinate system for the aircraft body and axes x*, y* and z* refer to the wing coordinate system. The wing with nacelle is defined in wing coordinate system and is placed in the body system with x and z translations of 13.661 in. and -1.335 in respectively with a dihedral of 4.787 degrees. The nacelle is located at 8.189 in. from the wing origin. The projected wing semispan is 23.0571 in. The wing is defined by a number of airfoil sections at different stations along the wing span as shown in Figure (3). The shape of the airfoil at each station is selected based on the aerodynamics and holds the shape of the wing.
Fig. (3) DLR-F6 wing showing different airfoil sections [5] Figure (3) shows a number of airfoil sections that are defined at different η along the wing span, where η is the normalized coordinate defined as η =
y* s*
.
The front spar is usually positioned at 15% of chord and the rear spar at 65% of chord measured from the leading edge. The enclosed area between the spars as shown is called the wing-box.
6
In order to test the optimization, the wind tunnel geometry is scaled by a factor λ=20 to build an approximately realistic aircraft model. The scaled model dimensions of the wing are given below: The wing reference area for the scaled model is S=90148 in.2 and the semi-span in wing coordinate system is s*= 463.3 in. The average chord length of the wing is C av = 97.746 in. and the mean aerodynamic chord length is C mac = 111.18 in.
7
Initial Sizing of Wing-Box Structure
Section II
For initial sizing of the wing-box structure, it is necessary to define the different loads experienced by an aircraft during its maneuvering.
II-1 Loads Acting on an aircraft wing-box: An aircraft wing structure is mainly subjected to three kinds of loading a) Aerodynamic loads in the form of lift and drag forces and pitching moments. b) Concentrated forces due to landing gear connections, power plant’s nacelle connections, connections to the fuselage, connection with the controlling surfaces structures like ailerons…etc. c) Body forces in the form of gravitational forces and inertia forces due to wing structural mass. The stress analysis of the wing-box requires a complete identification of all the loads acting on its structure.
II-1-a Aerodynamic Loads: Generally, an aircraft flying in air is subjected to aerodynamic loads [2, 3]. The lift produced by the aircraft balances its weight and the drag force balances the thrust produced by the aircraft as shown in Figure (4).
Fig. (4) Lift & Weight and Drag & Thrust balancing the Aircraft [5]
8
Figure (5) shows different rotational motions exhibited by an aircraft. Pitching moment is expressed about the center of gravity of the aircraft.
Fig. (5) Pitch, Yaw and Roll motions of an Aircraft [5] The loads experienced by an aircraft wing are usually expressed in terms of aerodynamic coefficients [2], namely, the lift coefficient ( C L ), the drag coefficient ( C D ), the pitching moment coefficient ( C M ), the normal force coefficient ( C N ) and the tangential force coefficient ( CT ). All these coefficients are usually calculated using CFD solutions, as shown in figure (6), and are verified by wind tunnel tests since testing an actual aircraft is quite cumbersome and expensive.
Fig. (6) DLR-F6 CFD model, used to calculate aerodynamic loads [5]
9
The above mentioned aerodynamic coefficients are all defined as below:
CL = CD = CM =
CN = CT =
L q∞ S D q∞ S M q∞ S c
N q∞ S T q∞ S
(1) (2) (3) (4) (5)
Where S is the wing reference area; for airfoils a reference length is required rather than an area; thus the chord or length of the airfoil section is used for this purpose. q∞ is the free stream dynamic pressure calculated as:
q∞ =
1 ρV 2 2
(6)
Where ρ and V are the density of air and speed of the aircraft (calculated from Mach number, M) respectively. Since the speed of sound varies with the density of air, it is required to determine the density of the air through which the aircraft is flying. To compute this, the chart shown in Table (1), called the International Civil Aviation Organization Table (ICAO) is always used. It can be noticed that as the altitude increases, the density of air decreases and so does the speed of sound.
10
Table 1: Variation of density of air and speed of sound with altitude Altitude (ft)
Density of Air
Speed of Sound
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000
( kg / m 3 ) 1.2249 1.1894 1.1548 1.1208 1.0878 1.0554 1.0239 0.9930 0.9626 0.9332 0.9044 0.7709 0.6524 0.5488 0.4581 0.3798 0.3015 0.2370 0.1865 0.1469
(m/s) 340.4076 339.2758 338.0926 336.9094 335.7262 334.5429 333.3083 332.1250 330.9418 329.7072 328.5239 322.4021 316.1773 309.7982 303.2647 296.6284 295.1880 295.1880 295.1880 295.1880
When a wind tunnel is used to collect aerodynamic data, first the actual lift force L is measured then it is converted to a non-dimensional coefficient C L using equation (1). All the complex aerodynamics has been hidden away in the lift coefficient. It is noticed that C L depends on the angle of attack ( α ), Mach number (M) and Reynolds’s number (Re). To summarize, the lift coefficient it becomes a function of three variables, C L = f ( α , M, Re)
(7)
The CFD solution for wing-body-pylon-engine (wing-mounted engine) case giving the lift coefficient and pitching moment coefficient for DLR-F6 aircraft wing at test conditions of Mach = 0.75; CL=0.5 (CL is the overall lift coefficient);
α = -0.0111o and Re = 0.300E7 is given in Tables (2) and (3). 11
Table 2 Variation of Lift Coefficient (vs) η
η=
y* s*
0.1274 0.1651 0.2029 0.2409 0.2793 0.3180 0.3572 0.3971 0.4377 0.4792 0.5219 0.5657 0.6111 0.6582 0.7074 0.7589 0.8133 0.8711 0.9330 1.0000
CL
c
0.4328 0.4580 0.4784 0.4908 0.4926 0.4864 0.5141 0.5483 0.5698 0.5899 0.6068 0.6212 0.6340 0.6439 0.6502 0.6553 0.6504 0.6353 0.5861 0.4832
158.8756 149.9319 140.9653 131.9525 122.8697 113.6916 104.3917 94.9413 91.3885 88.1641 84.8560 81.4495 77.9280 74.2717 70.4573 66.4566 62.2348 57.7487 52.9427 47.7443
Tables (2) shows the variation of the local lift coefficient at different stations along the wing span where “ C L ” is the local lift coefficient at a specific span coordinate and “c” is the local chord length at that span coordinate. Figure (9) shows the variation of the lift coefficient along the wing span. From table (2) and by using equation (1), the lift force per unit length along the wing span can be calculated, as shown in figures (7) and (8).
12
Fig. (7) lift coefficient C L (vs) normalized wing span coordinate η
Fig. (8) Lift Force L per unit length (vs) η
13
Table (3) Pitching Moment Coefficient about Local Quarter Chord (vs) η
η=
y* s
C Mqc
*
0.1274 0.1651 0.2029 0.2409 0.2793 0.3180 0.3572 0.3971 0.4377 0.4792 0.5219 0.5657 0.6111 0.6582 0.7074 0.7589 0.8133 0.8711 0.9330 1.0000
-0.0958 -0.0929 -0.0939 -0.0987 -0.1071 -0.1201 -0.1374 -0.1461 -0.1381 -0.1325 -0.1280 -0.1249 -0.1231 -0.1221 -0.1228 -0.1222 -0.1203 -0.1165 -0.1128 -0.1093
Table (3) shows the values of pitching moment coefficient about quarter chord length along the wing span. These data are represented graphically in Figure (9). From table (3) and by using equation (3), the pitching moment about quarter chord length can be calculated, as shown in Figure (10).
14
Fig. (9) Pitching moment coefficient about quarter chord C Mqc (vs) η
Fig.(10) Total Pitching Moment (about Quarter Chord) (vs) η Integration of the curve in Figure (8) along the spanwise direction gives the shear force distribution on the wing as shown in Figure (11). The bending moment distribution along the wing span can also be obtained by integrating the shear force distribution as shown in Figure (12).
15
Fig. (11) Shear Force (vs) η
Fig. (12) Bending moment (vs) η The loads calculated in Figures (10), (11) and (12) are still not the actual DESIGN loads. They need to be scaled up by applying suitable scaling factors as these
16
loads are too small to use for sizing the wing box. The conditions of Mach = 0.75; CL=0.5 and Re = 0.3E7 is a cruise condition. Hence, a design condition of 2.5g maneuver is considered here and the obtained loads are multiplied by a factor 2.5 to make them the actual DESIGN loads. Also an additional safety factor of 1.5 is applied over these loads. All these external aerodynamic loads will be resisted by internal reactions in the wing structure. The design of the stiffened panels is based on the assumption that the stringers are the members which are responsible about the bending resistance, while the skin is designed to just carry in plane stresses in the form of in plane shear stresses and tensile stresses, but its resistance to compressive stresses is very limited due to its instability under slightly compressive loads. The variation of the bending stress along the stiffened panels will generate a flexural shear flow in the plane of the airfoil.
II-1-b Concentrated Loads: Concentrated forces acting on the wing-box structure are acted primarily on the wing ribs, which by its turn redistribute these forces in to the wing section in the form of shear flow. The following is a summary for different cases of concentrated loads and the corresponding rib stiffeners arrangements: 1) If the concentrated force is applied in the plane of the rib, then the stiffener should be aligned with the line of action of the force. 2) If placing the stiffener to be aligned with the load is impossible due to some openings in the rib, cutouts…etc, then placing two inclined stiffeners is also acceptable, since each stiffener will carry a component of the load in its direction. 3) If the load is out of plane of the rib, then placing three stiffeners perpendicular to each other is also acceptable since each stiffener will carry a component of the force in its direction, as shown in figure (13)
17
Fig. (13) the stiffeners arrangement in a shear web subjected to out-of-plane concentrated force [6] 4) If the load is normal to the web, then design of stronger flanges to carry the load in bending then transfer it to the web.
II-1-c General definition of wing station external loads: A wing station “j” is subjected to two type of loading 1- Shear forces in the form of: a.
Vertical shear force VZj This vertical shear force includes:
(i)
The total lift summation from the wing tip till the ‘jth’ wing station which can be obtained from figure (11) for the DLF-6 aircraft.
(ii) Wing structural weight (body forces) included in the wing portion extending from the wing tip till the ‘jth’ station. It is important to note that in the conceptual design stage the size of the wing parts is not yet determined. Accordingly, the weight of the wing portions will not be available. Alternatively, an approximate value for the distribution of the wing weight along the wing span can be obtained from previously designed
18
airplanes, or the weight may not be included in the initial sizing process, then it can be included later then an iteration design process can be conducted for a suitable convergence for the wing weight. (iii) Inertia forces (body forces), where the mass of the wing portion structure must be multiplied by the acceleration of flight in the vertical direction. (iv) Non-structural mass forces due to the fuel tank weight…etc. in the form of weight and inertia forces. b.
Horizontal shear force V Xj This horizontal shear force includes:
(i)
The total drag summation from the wing tip till the ‘jth’ wing station.
(ii)
Inertia forces (body forces), where the mass of the wing portion structure must be multiplied by the acceleration of flight in the horizontal direction.
(iii)
Non-structural mass forces due to the fuel tank mass…etc. in the form of inertia forces.
2. Twisting moment a wing station “j” is subjected to twisting moment ‘ M j ’ the sources of this twisting moment are (i)
The pitching moment M qcj . The pitching moment about quarter chord location for DLR-F6 can be obtained from figure (10).
(ii)
Twisting effect of lift forces.
The lift force is always calculated with respect to the aerodynamic center
of
the
wing
cross-section
which
with
an
acceptable
approximation considered as the airfoil quarter chord location. This lift force at the quarter chord has a twisting effect with the value of VZj e X j
19
where e X j is the horizontal distance between jth station quarter chord and shear center. (iii)
Flange forces twisting moment ‘ M fj ’:
Since the aircraft has a tapered wing, which implies that the stiffeners are not perpendicular to the airfoil cross-section but they have inclination angle in the Y-Z plane as well as in the Y-X planes. These inclinations generate a flange force components in the three space directions F fj
X ,i
, F fj
Y ,i
and F fj
Z ,i
.
In the calculation of the shear flow around the airfoil cross-section the in-plane forces are of quite importance to the calculations. F fj
X ,i
F fj
Z ,i
and
are the ‘ith’ stringer in the ‘jth’ wing station flange forces in the X
and Z directions respectively. These forces are generating a flexural shear effect as well as a twisting effect on the airfoil cross-section. (iv) Twisting effect of the drag forces: The general shape of the airfoil is shown in figure (14) and the drag forces are always considered as acting horizontally through the airfoil chord line, as shown in the following figure
Fig. (14) Airfoil main lines
20
If the line of action of the drag forces is not passing with the airfoil shear center, then a twisting effect takes place with magnitude V Xj e Z j where eZ j is the vertical distance between wing station shear center and its chord line. (vi) Twisting effect due to wing portion weight: Once the wing weight included in the design process, a twisting effect of the wing portion weight must be introduced, the magnitude of this twisting moment is ‘ W j e0 ’ where W j is the weight of the wing portion extending from the wing tip till wing station ‘j’ and e0 is the horizontal distance between the airfoil centroid (center of gravity) and its shear center at that wing station.
II-2 Initial Sizing of Stiffened Panels Refer to reports one, two, three and four for details of stiffened panels sizing.
II-3 Initial Sizing of Wing-Box Ribs Refer to report five for details of wing rib stress analysis and design.
21
II-4 Initial Sizing of Wing-Box Spars and Spar’s Caps II-4-1 Literature reviews about spars and wing box Generally, there are two categories of approaches to deal with the design and analysis of spar and wing box. One is the traditional method, which was introduced in detail by Kuhn and al.[1], Bruhn [2], and Niu [3]. By using column and plate buckling and crippling analysis theory combined with empirical equations or curves, the dimensions of spars and wing box can be decided and optimized. The advantages of this method are simple and quick calculation with coarse accuracy. The main obstacle is that the empirical equations and curves can only be used to the specific materials and the configuration of structures. Therefore, the applications are limited. Another category is finite element analysis, which has been used more and more, especially for composite wing box structure. With the developments of integrated CAD/CAE software and reducing cost of computer hardware, a considerable amount of research has been conducted in this field, especially for Multi-Disciplinary Optimization (MDO). The main obstacle of this method is the increased computational cost and the unacceptable solution time because of the nonlinear analysis and the so many iteration procedures. Therefore, the objective of stage 3 is the section sizing optimization of spars and wing box by using traditional method and stage 4 is whole wing box optimization by using FEA method based on the results of stage 1 to stage 3.
22
II-4-2 Design of Spars Web and Caps Each wing box structure has affront spar assembly, a rear spar assembly and the ribs. The air load act directly on the wing covers which transmits the loads to the rib. The rib transmits the load in shear to the spar web and distributes the load between them in proportion to the web stiffness. The spar web and caps are mainly subjected to bending and shear loading. The depth of the web is usually large as compared to depth of the cap, therefore bending stress in the web are neglected. It is assumed that caps develop the entire bending resistance and shear flow is constant over the web.
Fig. (14) Spar Cap Assembly The spars are approximately located early in the design phase during the selection and layout of the wing box size. A natural tendency is to locate the front spar at a constant chord location, between 5% to 20% chord. The front spar location should be selected to ascertain the space provisions in the leading edge device and to maximize the box volume for fuel containment and structural rigidity. The rear spar is usually located between 60% to 80% chord. The rear spar location is subject to as many or more influences as the front spar. Spar caps
23
are used to connect spar web to the skin of the wing box as shown in Figure 1. Generally, T-type and L-type caps are used for the aircraft structures. The cross sectional area and other design parameters for these sections are: For a T-section as shown in Figure 2
Fig. (15) T-Type Spar Cap
Acap = B1T1 + B2T2 Y=
Ix =
(8)
(B T
+ T2 B2 (2T1 + B2 ) ) (2(T1 B1 + T2 B2 ))
2 1 1
B1 B3 ( B2 + T1 )3 − 2 ( B1 − T2 ) − Acap ( B2 + T1 − Y ) 2 3 3
(9)
(10)
where Acap , Y and I x are cap cross sectional area, second moment of area about centroidal axis and centroidal distance from the top of cap. For a L-section as shown in Figure 3
Acap = B1T1 + B2T2 Y=
(11)
B22 + ( B1 − T2 )T1 2( B1 + B2 − T1 )
(12)
24
(T1 ( B2 − Y ) 3 + B1Y 3 − ( B1 − T2 )(Y − T1 ) 3 ) Ix = 3
(13)
Fig. (16) L-Type Spar Cap The he , hu and he are calculated as:
he = h − 2Y ,
hu = he − ( B2 − Y ) ,
hc = hu −
( B2 − Y ) 2
(14)
The applied shear flow in the web can be written as:
q=
V he
(15)
where V is the applied shear force on the beam. The applied shear stress in the web can be calculated as:
fs =
q t
(16)
where t is the thickness of web.
II-4-2-a Objective Function The objective of the optimization problem is to minimize the mass of the spar web and caps assembly while preventing against any type of failure. The design
25
variables are the thickness of the spar and dimensions of the caps. The objective function can be stated as: Objective Function= (t h + 2 Acap ) L
(17)
where h is the height of the spar web and L is the length of the bay.
II-4-2-b Constraints In order to optimize the wing box structure, the design must satisfy a set of constraints, e.g. material failure and buckling must not occur anywhere within the configuration. The present work is mainly concentrated on following constraints II-4-2-b -1) Spar Web Failure II-4-2-b-2) Spar Cap Failure II-4-2-b-2-1) Crippling failure II-4-2-b-2-2) Bending failure
II-4-2-b -1 Spar Web Failure The two basic types of web design are shear resistant type and diagonal tension field type. A shear resistant web is one that carries its design load without buckling of the web. The design shear stress is not greater than the buckling shear stress for the individual web panels and the web have sufficient stiffness to keep the web from buckling as a whole. It is realized that the buckling web stress is not a failing stress as the web will take more before collapse of the web take place, thus in general web is not loaded to its full capacity for taking load.
26
Therefore, diagonal tension type web are generally used for the design of spars of an aircraft. In the diagonal tension webs, buckling of the web is permitted with shear loads being carried by diagonal tension stresses in the web. At the buckling load panel buckles into the diagonal folds and additional loading is taken by diagonal tension produced in these folds. The equations required for the analysis are presented here and the detail description of theory of incomplete diagonal tension can be referred from [1]. The stresses in the web subjected to incomplete diagonal tension depend on the diagonal tension factor which is measure of degree of loading of structure above its buckling strength. The diagonal tensional factor can be calculated as:
⎧⎪ ⎛ f ⎞⎫⎪ K = tanh ⎨0.5 log10 ⎜ s ⎟⎬ ⎜ F ⎟⎪ ⎪⎩ ⎝ s_cr ⎠⎭
(18)
where Fs_cr is the shear buckling stress of the web. The shear buckling stress of the web can obtain by following steps: (i) Calculate the flat plate buckling coefficient K s for inplane shear loading using the following polynomial:
K s = 0.001 r 6 − 0.0293 r 5 + 0.3401 r 4 − 2.0808 r 3 − 13.668 r + 16.497
(19)
where
r=
hc d
(ii) Calculate the following ratio:
FS_R
⎛t⎞ = K S E⎜ ⎟ ⎝d⎠
2
27
(20)
(iii) The shear buckling stress of the web is obtained from the following polynomial approximations:
Fs_cr = 42.5
Fs_cr = −3 × 10−12 (Fs_R ) + 2 × 10− 9 (Fs_R ) − 8 × 10− 7 (Fs_R ) 6
5
4
+ 0.0001(Fs_R ) − 0.0141(Fs_R ) + 0.8472Fs_R + 12.819 3
FS_R > 197.85 FS_R > 0 & FS_R < 197.85
2
Fs_cr = 0
FS_R = 0 (21)
The maximum shear stress in the web corresponding to the above calculated diagonal tension factor can be calculated as: f s_max = f s (1 + K 2 C1 )(1 + KC 2 )
(22)
where C1 and C 2 are the stress correction factors. The factor C1 is to allow for the fact that the angle α of the diagonal tension differs from 45 degrees and can be obtained as:
C1 =
1 −1 Sin (2α )
(23)
The factor C 2 is the stress concentration factor arising from flexibility of the cap and can be obtained as: C 2 = 0.0013(wd)6 − 0.0211(wd)5 + 0.1156(wd) 4 − 0.2434(wd)3 + 0.2267(wd) 2 − 0.0713(wd)
wd < 4
C2 = 1
wd > 4
(24) where wd is the cap flexibility factor and can be obtained as: 1
⎛ ⎞4 t ⎟⎟ wd = d⎜⎜ + 4 h (I I ) ⎝ e T c ⎠
28
(25)
where I T and I C represent the second moment of area of the tension and the compression flanges respectively (with respect to their centroidal axis). The allowable maximum shear stress in the web can be obtained from the following relation: For α = 40o Fs,all = −52.612K 6 + 186.3K 5 − 243.12K 4 + 124.62K 3 + 4.2935K 2 − 27.272K + 34.85 (26)
For α = 45o Fs,all = −1.1689K 6 + 41.756K 5 − 92.897K 4 + 53.707K 3 + 19.277K 2 − 28.05K + 34.924 (27)
The allowable maximum shear stress at α = 45o is good approximation for the most of the problems. The optimization problem is constrained such that maximum shear stress is less than the allowable maximum shear stress. The optimization constraint can be written as:
f s,max − Fs,all ≤ 0 or f s,max Fs,all
−1 ≤ 0
(28) (29)
II-4-2-b-2 Spar Cap Failure The crippling and bending failures are two main modes of failures in the spar cap. The cap is designed such that its resist both types of failures
29
II-4-2-b-2-1 Crippling Failure The cap resists two types of axial compressive stresses, compressive stress caused by bending moment and compressive stresses caused by diagonal tension. The compressive stress in the cap caused by beam bending moment can be written as: fb =
M h e A cap
(30)
The compressive stress caused by the diagonal tension in the web can be written as: fF =
KV 2A cap tanα
(31)
The crippling failure in the cap is caused by combination of f b and f F . To compute the allowable crippling stresses of the cap, the section is broken down into individual segments and each segment n has width a width b and a thickness t and will have either one or no edge free. The allowable crippling stress for each segment n is found from the applicable material test curve or from the following empirical formulas: If segment n has free edge: Fccn
⎡b = 0.6424Fcyn ⎢ n ⎢⎣ t n
Fcyn ⎤ ⎥ E n ⎥⎦
−0.788
(32)
If segment n has no edge free: Fccn
⎡b = 1.1819Fcyn ⎢ n ⎣⎢ t n
30
Fcyn ⎤ ⎥ E n ⎦⎥
−0.7882
(33)
where Fcyn and E n are the allowable compression yield stress and the modulus of elasticity of segment n. The allowable crippling stress for the entire section is computed by taking a weighted average of the allowable for each segment: Fcc =
∑b t F ∑b t
n n ccn
(34)
n n
II-4-2-b-2-2 Bending Failure In addition to the compressive stress, the cap is also subjected to bending moment. The bending moment is known as secondary bending moment and produce by incomplete diagonal tension in the web. The secondary bending moment can be obtained as:
M max
K C3 fS t d 2 tanα = 12
(35)
where C3 is the stress concentration factor and can be obtained from the following equation:
C3 = 5 × 10−5 (wd)6 − 0.001(wd)5 + 0.0084(wd)4 − 0.0332(wd)3 + 0.0341(wd)2 − 0.0092(wd) + 1 C3 = 0.6
wd < 4 (36) wd ≥ 4
The secondary bending stress in the cap can be obtained as: f sb =
M max ( B2 + T1 − Y ) I x , cap
(37)
The cap is subjected to both compressive and bending stress simultaneously, therefore the margin of safety of the cap is combination of both compressive and bending failure. The margin of safety for the cap can be calculated as:
31
MS =
1 f b + f F f sb + Fcc Ftu
−1
(38)
The constraint on the optimization problem is imposed such that margin of safety of cap is always greater than zero. The nonlinear optimization constraint can be written as:
1−
1 f b + f F f sb + Fcc Ftu
≤0
(39)
II-4-3 Numerical Validation The optimization problem formulated above is validated by comparison with design example solved by Niu [3] and Abdo [4]. The design parameters for the problem are:
h = 14in., V = 30000lb, d = 8.0in and M = 560000lb.in The material properties of the web and caps are: Web-7075-T6 bare sheet
Caps-7075-T6 Extrusion
E = 10.5 ×10 6
Ftu = 80000 psi
Fcy = 71000 psi
Fsu = 48000 psi
E = 10.7 ×10 6
Ftu = 82000 psi
Fcy = 74000 psi
Fsu = 44000 psi
In the present case, the optimization problem is solved by considering three different cases: Case 1: All design variables are independent Case 2: Lengths of both flanges are constrained to be equal.
32
Case 3: Lengths and thickness of both flanges are constrained to be equal The results obtained from the present optimization algorithm along with those obtained by Niu [3] and Abdo [4] are presented in Table 4. The optimum dimensions of the spar cap are given in Table 5. Table 4: Results obtained by Niu and M.Abdo
Niu M.Abdo Case 1 Case 2 Case 3
Acap 0.918 1.01 0.7946 0.8177 0.8181
T 0.085 0.085 0.0855 0.0822 0.0821
Total Cross. Area 3.026 3.21 2.626 2.6869 2.6873
K 0.25 0.26 0.2804 0.3047 0.3047
he 12.08 12.12 12.122 12.784 12.798
hu 10.4 10.51 10.85 11.52 11.48
hc 9.7 9.71 10.21 10.889 10.828
Fscr 8772 8145 7680.1 6702.2 6706.3
Fs,all 29700 29383 29255.8 28993.9 28995
fs,max 29387 29242 29273 29009 29994
MS (web) 0.01 0.0048 0.0006 0.0005 0
Table 5: Dimensions of cap Case 1 Case 2 Case 3
B1 1.0986 1.6315 1.669
B2 2.0491 1.6315 1.669
T1 0.1617 0.2396 0.2451
T2 0.3011 0.2616 0.2451
It can be seen that total cross sectional area of the web-caps assembly is reduced significantly by using present method. It can also be observed that the diagonal tension factor obtained at the optimum design using present method is more than that obtained by both Niu [3] and Abdo[4], and web is subjected to large diagonal tension. The minimum cross sectional area is obtained for case 1 where all design variables are independent. For case 2, an additional constraint is imposed on the optimization algorithm such that lengths of both flanges are equal. The additional constraint is imposed to obtain more symmetrical design and results in decrease in the number of design variables.
The additional
constraint results in a little heavier design than the previous case, but still much
33
MS Cap 0.03 0.067 0.0002 0.0004 0.0003
lighter than that obtained by Niu [3] and Abdo[4]. The thickness of the web is decreased and web is subjected to higher diagonal tension field. To obtain more symmetrical design, case 3 is considered. It is assumed that both flanges have equal length and thickness. It can be that very small increase in the mass of the structure is observed by imposing this additional constraint, and insignificant change has been observed in the diagonal tension factor and thickness of the web.
II-4-3.4 Conclusion The optimization problem formulated above generate very accurate results, and even better than other formulations. The optimization algorithm will be used to size web spar and spar caps at each section of wing box. Furthermore, the comparison between T type and L type section will be also be made, and effect of numbers of uprights on the optimum design will also be investigated. The optimum dimensions of spar web and caps obtained from optimization process will be used to build conceptual wing box model.
34
Section III
3D Finite Element Model of DLR-F6 Aircraft Wing-Box Structure, Created in PATRAN and Analyzed in NASTRAN
Since the thickness of the skin as well as the width of the skin-stringer panels are the two design variables of the optimization process of stages one and two, then the output of these two stages are the dimensions of the skin thickness along each panel pitch at each wing station. These dimensions are determined through an optimization process for mass minimization as an objective functions. Considering station 21 as an example, the output of stages one and two is tu=[0.08 0.08 0.08 0.08 0.08 0.08] bu=[3.64 3.64 3.64 3.64 3.64 3.64] tL=[0.07 0.07 0.07 0.07 0.07] bL=[4.22 4.22 4.22 4.22 4.22]
Where “tu” and “tL” are the wing skin thicknesses along the upper and the lower skin panels, respectively. While “bu” and “bL” are the width of the skin-stringer panels along the upper and the lower skin, respectively. Using the stringer’s pitch along each skin-stringer panel, the number of stringers as well as the coordinates of the stringers along the upper and the lower wing skin are determined in wing coordinate system. Taking station 21 as an example, the coordinates of the stringers along the upper and the lower skin are presented as x=[245.02 246.84 250.48 254.12 257.76 261.40 265.04 268.98 268.98 264.01 259.79 255.57 251.35 247.13 245.02] z=[43.96 44.16 44.48 44.68 44.79 44.80 44.71 44.52 40.41 39.70 39.29 39.06 39.05 39.25 39.44]
These two vectors represent the x and z coordinates of each stringer at station 21. It is important to mention that, the fist and last component in the x and z vectors represent the location of the front spar upper and lower caps respectively. While the 8th and the 9th component represent the rear spar upper and lower cap respectively. The rest represent the coordinates of the stringers locations along the upper and the lower skin. To insure moment of inertia maximization, a set of relations are adopted to relate skin thickness and panel width with the rest dimensions of the skin-stringer panel. Recalling from the first report the details concerning the ‘Z’ stringer
35
Fig. (17) Panel geometry definition using ‘Z’ stringer [7] Where
ba = 2.08t s + 068
if t s < 0.3
ba = 1.312
if t s > 0.3
t a = 0.7t s for equal flange stringers b f = ba and t F = ta
(40)
⎛b ⎞ bw = ⎜⎜ e ⎟⎟( Ast − 1.4ba t a ) ⎝ ts ⎠ ‘ Ast ’ is the stringer cross-section area, and it can be represented as
t ⎞ ⎛ Ast = (bw − t f )t w + 2⎜ b f + w ⎟t f 2⎠ ⎝
(41)
And ‘ be ’ is the effective width of the skin [10], where, be = t s K
ηE sk σ sk
(42)
where ‘K’ is the skin buckling coefficient and it takes the values
36
K = 3.62
for
bs < 40 ts
K = 6.32
or
for
bs > 110 ts
Between the above values there is a gradual transition, as plotted in this figure
bs ts
Fig. (18) Variation of compression panel skin buckling constant with skin crosssection aspect ratio [6] ‘η ’ in equation (18) is the plasticity reduction factor which is determined using the following equation
η=
Et sk
(43)
E sk
Where ‘ E sk ’ and ‘ Et sk ’ are the elastic and tangent modulii of the skin, respectively While ‘ σ sk ’ is the skin axial stress. For practical use, the design curves for the skin stringer panels can be used, where
37
bs =1 be
for high value of load index
bs = 1.1 : 1.3 be
N L
N for low values of load index L
(44)
Where ‘N’ is the axial load intensity, and it can be calculated using the equation N=
Mb cs hs
(45)
And ‘L’ is the effective column length, or the distance between two successive ribs. The output of station three is the dimensions of the spars. A wing spar is composed of a spar web, upper and lower spar cap and a group of uprights to reinforce the spar against collapse when it is subjected to incomplete diagonal tension. Also, the output of the optimization process of the wing rib is the thickness of the rib web, number and center location of lightening holes and the diameters of the lightening holes. Refer for the report five for more details. All there dimensions are employed to build the 3D finite element model of the wing.
III-1 Organization of the Finite Element Model The wing-box is divided into 20 bays extending between 21 stations as shown in Figure (20). The stations are named from station 1 at the wing root to station 21 at the wing tip. While the bays are named from bay 2001, extending from station one to station two, to bay 2120, extending from station 20 to station 21. All the elements in a bay are numbered so that the first four digits in its index represent the name of the bay (i.e. an element in bay 2120 has the index 2120xxx).
38
y
Fig. (19) DLR-F6 wing lay-out MSC.PATRAN is used as the modeler to build the finite element model. Each bay is grouped into five groups namely the skin, the stringers, two groups representing each rib bounding the bay in the span wise direction and a group for the load card and its rigid elements which are used to load distribution.
III-2 Building the Finite element model MSC.PATRAN is used to build the finite element model. Following are the steps used to create bay 2120, extending between stations “21” and “20” in the current wing-box 1- in MSC.PATRAN a new data base is created and named DLR-F6-wing-box3D-FEM 2- The model tolerance is set to the default. Analysis type is set to “structure” while the analysis code is chosen to be MSC/NASTRAN. 3- The set of coordinates representing the locations of the stringers and the spar caps at stations “20” and “21” are used to generate these points in MSC.PATRAN as shown in Figure (20).
39
Fig. (20) Points representing the locations of the spar caps and the stringers at stations “20” and “21” 4- A new group is created and named “2120_stringers”, then post this group as the current group. Using the points created in the previous step a group of lines is generated between pairs of points extending from station “21” to station “20”. The order of creating these lines must start by the lines representing the spar caps, then the lines representing the stringers are created in the order starting from the points near the front spar then proceed towards the rear spar. The reason of this order is that the stringers run out always takes place at the rear spar, i.e. a difference in the number of stringers between two stations means that there is a stringer run out equal to the difference between the number of stringers between the two successive stations and these run outs take place at the rear spars, as shown in the next figure.
40
Fig. (21) Group of lines representing the spar caps and stringers in bay 2120 From figure (21) it can be noticed that there is a point on the top skin and another one on the lower skin that are not employed in generating the stringer lines, this indicates that these two points are a run out of two stringers in bay 2019. 5- A group is created and named as “2120_skin” and posted as the current group. Use the lines generated in the previous step to generate surfaces extending between adjacent lines in the chord wise direction as shown in the following figure.
41
Fig. (22) Group of surfaces representing the bay skin and the spars webs 6- Once the geometry of the bay is created, finite elements can be generated. It is well known that increasing the number of elements in the model enhances the accuracy of the results but it increases the model cost. Accordingly, it is required to keep the minimum number of elements necessary to obtain acceptable accurate results. To do so, the number of elements along the bay is selected to be two elements in the span wise direction, and after finishing the whole bay model, the bay is tested for an arbitrary load and the results are obtained. Then the number of elements is increased to three in the span wise direction and the model is resubmitted to NASTRAN for analysis. The results obtained are compared with the results obtained from the pervious step. If a significant change is obtained in the results then, it is required to re-increase the number of elements and test again. A change in the result with in 0.05% doesn’t
42
require additional refining of the model. It has been found that three elements in the span wise direction results in acceptable results. 9- Elements Properties: After creating the finite elements, the elements properties should be applied. The stringers are modeled by beam elements with a Z shape cross-section. The details of the Z-shape cross-section are shown in the next figure.
Fig. (23) The Z-shape cross-section of the beam element used in the PBEAML Card for stringer modeling [8] A comparison between the dimensions of this Z-shape cross-section and the dimensions obtained from the optimization process in stages one and two, shows that these DIM1, DIM2, DIM3 and DIM4 dimensions can be calculated by simple transformations as follows DIM 1 = ba −
tw 2
(46)
DIM 2 = t w
(47)
DIM 3 = bw − t a
(48)
DIM 4 = bw + t a
(49)
Since the dimensions of the stringers vary from one station to the other, an interpolation process is used to obtain the dimensions of all elements between stations.
43
This is done by defining the dimensions in PATRAN as fields, where a local coordinate is created at each station with its Z-coordinate directed in the span wise direction. Set of fields are created in PATRAN defined in the station local coordinate, representing the variation of the dimension in the span wise direction. As an example, consider the dimension DIM1 of the Z-shape stringer extending between stations 21 and 20, this dimension is defined in PATRAN as a field on the form
DIM 1 = DIM 1 _ 20 + ( DIM 1 _ 21 − DIM 1 _ 20) Z
(50)
Similarly for all the other dimensions. The spar caps are also modeled by beam elements but with L-shape cross-section as shown in the following figure
Fig. (24) The L-shape cross-section of the beam element used in the PBEAML Card for spar caps modeling [8]
44
Fig. (25) the model stringers after applying the properties in PATRAN The skin is modeled by SHELL elements, where the thickness of the shells is also defined by fields representing the variation of the skin thickness in the span wise direction. 7- modeling of ribs: a group of points is generated to represent the perimeter of the rib, these points have the same y-coordinate of the station, while its x-coordinate has the same xcoordinate of the corresponding stringer, while its z-c00rdinate can be defined by the following equation
z r = z s − DIM 4 s
(51)
Where z r is the z-coordinate of the rib perimeter point, z s is the z-coordinate of the corresponding stringer while DIM 4 s is a dimension belongs to the stinger corresponding to this rib point.
45
The rib web is modeled by QUAD4 elements with PSHELL card for its properties. While the perimeter of the rib and the lightening holes are reinforced by beam elements. The following figure shows a complete bay modeled in PATRAN.
Fig. (26) Complete bay modeled in PATRAN
III-3 Model Verification Early model verification is very important before proceeding for the whole finite element model. Early detection of errors is very important, since detection of errors in the advanced stages is very costly and time consuming. The model can be verified by either of the following two methods a) Model verification. b) Modeling methodology verification.
46
a) Model verification: The finite element model of the DLR-F6 wing-box bay, created in PATRAN in the previous section, can be tested in NASTRAN for an arbitrary value of loading. Then the result obtained from NASTRAN is compared with the analytical solution of such model with the same loading. The complementary internal virtual work theory of idealized beams is used to calculate the deflection of this wing bay under the effect of a flexural bending force. NASTRAN Analysis The model created in the previous section is loaded by an arbitrary load of 1000 lbs acting at the bay centeroid. To apply a loading to the bay at its centroid, a grid point is created at the section centroid, then this grid point is connected to the skin-stringers connectivity grid points by a group of RBE2 elements, with its independent degrees of freedom are at the centroidal grid point, and its dependent degrees of freedom are at the skin-stringers connectivity grids. A load of 1000 lbs is applied in the negative z-direction at the centroidal grid point. The model is fixed at all the skin-stringers connectivity grid points of the opposite bay station. Then, the model is submitted to NASTRAN for linear static analysis, which resulted in a deflection of 0.0126 in., as shown in the following figure.
47
Fig. (27) Deflection of a wing-box bay due to an arbitrary force loading it is clear from figure (31) that the deflection is the result of superposition of a combined loading. Since the bay section centroid does not coincide with the its shear center, then the force applied at the centroid has a bending, shear as well as a torsion effect. The deflection obtained from the NASTRAN analysis is verified by the analytical results. b) Modeling methodology verification. Another method to verify the finite element model is to verify the modeling methodology it self. By solving the deflection of simple box beam structure subjected to an arbitrary loading, then modeling the same structure and analyzes it in NASTRAN. If the results coincide, then the modeling procedure is correct.
48
III-4 Complete 3D Finite Element Model of the DLRF6 Wing-Box Once the model is verified, the work can proceed towards creating the full finite element model of the wing-box, as shown in the following figure
Fig. (28) DLR-F6 wing-box finite element model The model is submitted to NASTRAN for linear static analysis. The wing-box is subjected to static loading represents normal cruising conditions, applied along the wing-box elastic axis which produced a maximum deflection at the wing tip of magnitude 17.7 in. as shown in the following figure.
49
Fig. (29) DLR-F6 wing-box deflection due to cruising conditions loads Completing the entry of the material card in the bulk data file, to include the ultimate stresses of the material in tension and in compression, generates the margins of safety of the finite elements. The analysis showed that the margins of safety are in the zero one interval which indicates proper sizing of the model.
50
III-5
Post-processing
of
the
Wing-Box
Finite
Element Model: Deformation of an aircraft wing during flight has significant consequences on the aerodynamic performance. Predicting an accurate value of the bending and twisting of the wing in flight depends on the fidelity of the finite element model of the wing-box. Validation of the finite element model means making sure that the structural response of the model reproduces the structural response of the real wing within an acceptable accuracy. To find the deflection and the twisting experienced by the current wing-box due to the applied aerodynamic loads, the deflection and the twisting experienced by the wing-box elastic axis are plotted against the normalized span wise coordinate
Deflection in Z-direction
“η ” as shown in the following figure
η Fig. (30) Deflections in z-direction (vertical) experienced by the DLR-F6 wingbox elastic axis under the effect of cruising conditions aerodynamic loads
51
Deflection in X-direction
η Fig. (31) Deflections in x-direction (in plane bending) experienced by the DLR-F6 wing-box elastic axis under the effect of cruising conditions aerodynamic loads
52
Twisting around y-axis
η Fig. (32) Twisting angle around the y-direction (torsional) experienced by the DLR-F6 wing-box elastic axis under the effect of cruising conditions aerodynamic loads Since these deformations experienced by the wing-box are just an interpretation of the structural stiffness properties, then it is more convenient to calculate the model stiffness properties while the deformations can vary based on the loading conditions. To evaluate the equivalent moment of inertia and torsional rigidity of the model, two shear center nodes are created at the extremities of each wing bay, those two nodes are attached to the structure, as previously explained, by rigid bodies whose its dependent degrees of freedom are specified at an arbitrary number of grid points of the skin-stringers connectivity points, while its independent degrees of freedom are specified at a single grid point of the shear center. The next step is to load the node, which is towards the wing tip, by three sets of unit load moments. The first set moment is along the x-axis to calculate the vertical
53
bending moment of inertia. The second set is the moment along the y-axis to calculate the torsional stiffness rigidity and the third set is along the z-axis to predict the horizontal bending stiffness. The wing-box rotations in the x, y and z directions due to the corresponding applied moments are computed using NASTRAN and then the corresponding values of the stiffness are calculated using the equations
( EI x ) i→ j =
(GJ y ) i→ j =
( EI z ) i→ j =
(η j − ηi )S *
(θ
xj
− θ xi
)
(52)
(η j − η i )S *
(θ
yj
− θ yi
)
(53)
(η j − η i )S *
(θ
zj
− θ zi
)
(54)
Where S * = 463.2" is the semi-span of the DLR-F6 wing. η i and η j are the normalized coordinates of the two stations i and j respectively. The stiffness properties of the 20 wing bays are calculated and plotted against the wing normalized coordinateη , as shown in the following figures.
54
EI x
η Fig. (33) Distribution of the vertical bending stiffness of the DLR-F6 wing-box along its span
GJ y
η Fig. (34) Distribution of the torsional stiffness of the DLR-F6 wing-box along its span
55
EI z
η Fig. (35) Distribution of the horizontal bending stiffness of the DLR-F6 wing-box along its span
III-6 Model Stiffness Validation A methodology to estimate the stiffness distribution of a new wing using the stiffness distributions of Bombardier’s existing wings was developed by M. Abdo et.al. [9]. This methodology is based on a set of empirical relations that are generated for predicting the ideal stiffness of an arbitrary aircraft wing-box. To obtain those empirical relations which are applicable to all these wings, the data of the stiffness properties of a group of existing Bombardier’s aircrafts were normalized. One of the normalization techniques used is that, the stiffness of the existing wing structure is divided by the stiffness of a solid block material bounded by the leading and trailing edge of the wing, which referred to as ( EI ) CATIA , this is because CATIA was used for the calculation of these solid wing stiffness, as shown in the following figure.
56
Fig. (36) DLR-F6 wing airfoil sections Figure (39) shows the airfoils sections of the DLR-F6 wing at 21 wing stations. At each wing station, the airfoil section is padded, and then the measure tool bar is used to calculate the stiffness of each wing section. The data obtained from CATIA are used along with the empirical relations to predict the ideal stiffness properties of such aircraft wing-box. For EI x the behavior of the normalized stiffness appeared to be different outboard and inboard of the break in the plan form, consequently different relations were used to fit the data of the empirical relations, as follows
(EI x )FEM (EI x )CATIA (EI x )FEM (EI x )CATIA η Break =
=
R Break − RRoot (η − η Root ) + RRoot η Break − η Root
for η Root ≤ η ≤ η Break
(55) for η Break ≤ η ≤ 1
= RBreak
y Break
(56)
S* 57
η Root =
y Root
(57)
S*
Where R Root is the (EI x )FEM And RBreak is the (EI x )FEM
(EI x )CATIA
(EI x )CATIA
ratio at η = η Root
ratio at η = η Break
It was determined that RRoot = 0.03 and R Break = 0.1 provides an acceptable fit for the airplanes. For GJ y the following empirical relations was developed
(GJ y )FEM (GJ y )CATIA = 0.002
(58)
For EI z the following empirical relations was developed
(EI z )FEM (EI z )CATIA
= 0.0103η + 0.007
(59)
The previous empirical relations are used to predict the stiffness properties of the ideal wing-box of the DLR-F6 aircraft, and it is compared with the current FEM stiffness properties, as shown in the following figures.
58
+
+
+
Ideal FEM Model
o o o o Current FEM Model
EI x
η Fig. (37) Comparison between distributions of the vertical bending stiffness of the DLR-F6 wing-box along its span in the ideal FEM and current FEM
+
+
+
Ideal FEM Model
o o o o Current FEM Model
GJ y
η Fig. (38) Comparison between distributions of the torsional rigidity of the DLRF6 wing-box along its span in the ideal FEM and current FEM
59
+
+
+
Ideal FEM Model
o o o o Current FEM Model
EI z
η Fig. (39) Comparison between distributions of the horizontal bending stiffness of the DLR-F6 wing-box along its span in the ideal FEM and current FEM It is clear that there is a great agreement between the current NASTRAN FEM and the ideal FEM of such a wing, predicted by the empirical relations.
III-7 Wing-Stick Model of the DLR-F6 Wing-Box Full Finite Element Model Different levels of wing structural models have been used for wing static aeroelastic analysis and optimization, ranging from simple models based on analytical or empirical expressions to complex finite element structural models. The difficulty is to find or develop aeroelastic models that are sufficiently simple to be called thousands of times during optimization, but are sophisticated enough to accurately predict wing deformations in bending as well as twisting. Simplified beam finite element model of aircraft wing-box structure, also known as stick model, are often used for aerodynamics-structure interaction, such models can be used for either static or dynamic aeroelsatic analysis [9].
60
To build the stick model of the DLR-F6 wing-box, each wing bay is modeled by only one bar element representing the bay elastic axis, where the entry of that bar element are extracted as equivalent parameters from the full finite element model of that bay. Following is the entry necessary for the CBAR card used to construct the stick model in NASTRAN along with its PBAR card.
Where EID…element unique ID PID…element property unique ID GA, GB…grid point identification number, the two shear centers identification number of the wing bay two stations. X1, x2, x3…components of the element orientation vector MID…material unique identification number A…area of the element, equivalent area of the wing-box bay I1, I2, J…equivalent bending stiffness in two planes, and torsional stiffness of the wing-box bay K1, K2…equivalent area factors for shear of the wing-box bay All these equivalent parameters must be extracted from the wing-box full FEM. the equivalent stiffnesses of the wing-box are already obtained in the previous section where their values can be calculated from the following relations, based on a unit load deformations, where the loads are applied in the directions of the principle inertia of the wing-box cross-section
61
(
)
(
)
(
)
( I x ) i→ j
* 1 η j − ηi S = E θ x j − θ xi
( J y ) i→ j
* 1 η j − ηi S = G θ y j − θ yi
(
( I z ) i→ j =
(
)
(58)
)
(59)
* 1 η j − ηi S E θ z j − θ zi
(
)
(60)
While the equivalent area and equivalent area factors for shear needs further processing in NASTRAN and PATRAN to be calculated.
III-7-1 Evaluation of the Equivalent EA’s and GK’s of the Wing-Box The process of evaluating the equivalent area and shear factors span wise distribution is different from that of the bending and torsional stiffness evaluation. In this process NASTRAN is executed for three load subcases for each wing-box bay. In this process, the three degrees of freedom related to rotation in all skin-stringers connectivity grid points at all wing station are frozen. While the translation degrees of freedom are kept free. The elastic axis nodes are connected to the skin-stringer connectivity grid points by RBE2 elements with its dependent degrees of freedom specified at the skin and its independent degrees of freedom specified at the shear center grid points. The first subcase, a unit force in the y-direction, the axial direction, is applied at the bay shear center grid point, in order to calculate the equivalent area, as follows
( A)i→ j
(
)
* 1 η j − ηi S = E D y j − D yi
(
)
(61)
In the second subcase, the structure is loaded by a unit load at the shear center in the z-direction, to calculate the shear factor in the z-direction as follows
62
(K1 )i→ j
(
)
* 1 η j − ηi S = GA D z j − D zi
(
)
(62)
Similarly, the structure is loaded by a unit force in the x-direction to calculate the shear factor in the x-direction, as follows
(K 2 )i→ j
=
(
)
* 1 η j − ηi S GA D x j − D xi
(
)
(63)
Where “D” denotes the displacement deformation in certain direction. This process is applied to the wing box 20 bays, where all the elastic properties of the 3D model are extracted and applied to construct the stick model, the results of the calculations are shown in the following table Table (6) the PBAR card entries necessary to generate wing stick model
I1 1.0e+003 * 4.2220 1.0e+003 * 3.3290 1.0e+003 * 2.5298 1.0e+003 * 1.9032 1.0e+003 *1.3495 1.0e+003 *1.0322 1.0e+003 *0.7285 1.0e+003 *0.5558 1.0e+003 * 0.4592 1.0e+003 *0.3814 1.0e+003 *0.3071 1.0e+003 *0.2367 1.0e+003 *0.1850 1.0e+003 *0.1440 1.0e+003 *0.1156 1.0e+003 *0.0812 1.0e+003 *0.0610 1.0e+003 *0.0497 1.0e+003 *0.0406 1.0e+003 *0.0329
I2 1.0e+004 *2.5894 1.0e+004 *2.1679 1.0e+004 *1.6881 1.0e+004 *1.3030 1.0e+004 *0.8151 1.0e+004 *0.7390 1.0e+004 *0.5305 1.0e+004 *0.4432 1.0e+004 *0.3676 1.0e+004 *0.3013 1.0e+004 *0.2379 1.0e+004 *0.1741 1.0e+004 *0.1356 1.0e+004 *0.1059 1.0e+004 *0.0841 1.0e+004 *0.0576 1.0e+004 *0.0434 1.0e+004 *0.0347 1.0e+004 *0.0280 1.0e+004 *0.0224
J 1.0e+004 *1.0549 1.0e+004 *0.8692 1.0e+004 *0.6786 1.0e+004 *0.5353 1.0e+004 *0.3961 1.0e+004 *0.3099 1.0e+004 *0.2335 1.0e+004 *0.1774 1.0e+004 *0.1477 1.0e+004 *0.1213 1.0e+004 *0.0958 1.0e+004 *0.0716 1.0e+004 *0.0553 1.0e+004 *0.0431 1.0e+004 *0.0348 1.0e+004 *0.0233 1.0e+004 *0.0177 1.0e+004 *0.0145 1.0e+004 *0.0119 1.0e+004 *0.0098
A 31.0505 31.1960 29.5022 28.3718 26.8366 26.9245 26.2704 23.3872 20.9122 18.8147 16.3768 13.7504 11.6117 9.9506 8.9070 6.8692 5.7850 5.3652 5.0415 4.9649
K1 0.2303 0.2423 0.2374 0.2356 0.2311 0.221 0.2190 0.2048 0.2075 0.2056 0.2014 0.1935 0.1879 0.1890 0.1827 0.1792 0.1788 0.1790 0.1814 0.1728
K2 1.0244 1.0150 0.9807 0.8337 0.8602 0.8203 0.8243 0.8379 0.8386 0.8344 0.8366 0.8411 0.8454 0.8389 0.8357 0.8333 0.8280 0.7729 0.7679 0.7155
These data are used to generate the wing stick model in NASTRAN, as shown in the following figure
63
Fig. (40) Stick model of the DLR-F6 wing-box structure Another stick model is created to represent the ideal stick model of this aircraft, the two stick models are loaded with the same loading that was previously applied to the 3D FEM to test the behavior of the models, as shown in the following figure.
64
Fig. (41) the stick model deformed under the effect of the cruising conditions aerodynamic loads A comparison diagrams are generated to compare between the performances of the three models, namely the 3D FEM, the stick model representing the 3D FEM and the ideal stick model generated using the Bombardier aerospace empirical formulas, as shown in the following figure.
65
Vertical Deformations (y direction) [inch]
+
+
Ideal Stick Model
+
o o o o 3D FE Stick Model 3D FE Model
> > > >
η
Axial Deformations (y direction) [inch]
Fig. (42) Comparison between the 3D FEM, 3D FE stick Model and the Ideal stick model vertical deformation [inch]
+
+
+
Ideal Stick Model
o o o o 3D FE Stick Model > > > >
3D FE Model
η Fig. (43) Comparison between the 3D FEM, 3D FE stick Model and the Ideal stick model axial deformations (y direction) [inch] 66
Twisting deformations (around y-axis) [degree]
+
+
+
Ideal Stick Model
o o o o 3D FE Stick Model > > > >
3D FE Model
η Fig. (44) Comparison between the 3D FEM, 3D FE stick Model and the Ideal stick model twisting deformation (y direction) [degree] The comparison of the three models reveals that the design of the 3D FEM is a conservative design where its deformations are always less than the deformations experienced by the ideal stick model, generated using the Bombardier Aerospace empirical formulas, when the two models are loading with the same loading. Also the comparison curves validate the stick model generated using the stiffness properties extracted from the 3D FEM which is very useful in optimization and aeroelastic analyses processes.
67
References 1. Kuhn P., Peterson J. and Levin L. “A summary of diagonal tension, Part1methods of analysis,” NACA Technical Note 2661 2. Bruhn, E.F. “Analysis and Design of Flight Vehicle Structures”, Jacobs & Associates Inc., June 1973 3. Niu M. “Airframe Stress Analysis and Sizing,” Hong Kong, Conmilit Press Ltd., 1997. 4. Abdo M., Piperni P.,Isikveren A.T., Kafyeke, F. “Optimization of a Business Jet,” Canadian Aeronautics and Space Institute Annual General Meeting, 2005. 5. http://aaac.larc.nasa.gov/tsab/cfdlarc/aiaa-dpw, 3rd AIAA CFD Drag Prediction Workshop, San Francisco, 2006. 6. Bruhn E.F, Analysis and Design of Flight Vehicle Structures, Jacobs & Associates Inc., June 1973 7. M. Abdo, P. Piperni, F. Kafyeke “Conceptual Design of Stinger Stiffened Compression Panels”. 8. http://www.mscsoftware.com/support/online_ex/Library.cfm 9. M. Abdo, R. L’Heureux, F. Pepin and F. Kafyeke “Equivalent Finite Element Wing Structural Models Used for Aerodynamics-structures Intraction”, Canadian Aeronautic and Space Institute 50th AGM and Conference, 16th Aerospace Structures and Materials Symposium 28-30 April 2003. 10. M. Abdo, P. Piperni, F. Kafyeke “Conceptual Design of Stinger Stiffened Compression Panels”.
68
