Design Optimization for Concrete Beam Subjectedto Blast Loads

Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 10th International Conference on Civil and Architecture Engineering ICCAE-10-2014 Military Technical College Kobry El-Kobbah, Cairo, Egypt DESIGN OPTIMIZATION FOR CONCRETE BEAM SUBJECTEDTO BLAST LOADS Amr A. Shaheen1, Sameh Y.Mahfouz1, Mostafa S. Amin1, Metwally Abu-Hamed2 P P P P P P P Abstract This paper presents a design optimization of a concrete beam subjected to blast loading due to nuclear explosion. It focuses on sizing optimization problem of a concrete beam having elastic foundation and fixed ends. In the formulation of the optimization problem, the objective function is to minimize the maximum existing displacement. The design variables of the optimization problem are the depths of the concrete beam. The constraints are the concrete beam mass, the derivative of the depth and the plastic strain. A MATLAB code has been developed to obtain the pressure time history. This code is linked to the finite element software ANSYS to determine the displacement time history , strains, stresses and the best cross section shape of the concrete beam. Keywords: Shape Optimization, Blast loading, Nuclear Explosion and Concrete Beam 1. Introduction Few papers, available in literature, have focused on blast mitigation for concrete and steel sections. Dharaneepathy and Sudhesh [1] investigated the stiffener patterns on a square plate subjected to blast loads modeled using Friedlander‟s exponential function. Xue and Hutchinson [2] and Fleck and Deshpande [3] compared blast resistance of solid versus sandwich panels which were assumed to be ductile to withstand deformation caused by impulsive blast loads. Yen, et al [4] presented an experimentally validated dynamic analysis procedure utilizing LSDYNA and the ConWep The numerical results indicate that significant reduction in the maximum stress amplitude propagating within the protected components can be achieved by suitable selection of honeycomb material. Liang et al [5] investigated the optimum design of steel panels under blast loading by using a combined algorithm of the feasible direction method (FDM) and the Backtrack Program Method (BPM). This work ignores the effect of nonlinearity and deals with the blast loads as static loads. Main and Gazonas [6] studied the effect of an air blast on uniaxial crushing of cellular sandwich plates. Initial numerical results show that the capacity of the sandwich plates of shock mitigation can be improved by varying mass fractions of the front and back face sheets. Further, the authors investigate variation in geometry and 1 Egyptian Armed Force 2 Professor of steel structures, Dept. of Civil Engineering, Cairo University Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 shock mitigation capacity. This leads to an optimization formulation to maximize mass distribution of front sheet while minimizing back-face acceleration. From the literature, it is clear that no application in optimizing shape of concrete beams is studied to resist blast loads. Therefore, an approach and methodology to optimize the shape of concrete beam to withstand blast loading due to nuclear explosion is presented in the current paper. Here, difficulties are also treated in the optimization problem. These difficulties are: i) the nature of blast loading ii) transient dynamic response iii) non-differentiable, nonconvex and computationally expensive functions iv) the nonlinearities are. pressure time load, deformations and plasticity. The transient analysis generally requires more computer resources. Some preliminary work should be done to understand how nonlinearities affect the structure's response by carrying a static analysis first.Here, a MATLAB computer program has been developed to implement the nuclear explosion.In the shape design optimization process, two different optimization methods (i.e. zero order method and first order method) are used to find the optimum shape of a concrete beam. 2.Pressure loads of nuclear explosion The effects of an explosion can be distinguished in three ranges: i) contact detonation ii) near zone and iii) far zone of the detonation. The size of all these zones depends on many parameters, the most important one is the quantity of the explosive charge. Additional effective parameters are given in curves for the description of the different air blast using a rich body of experimental data see, Kingery [7] in 1984. The parameters are presented in double logarithmic diagrams with the scaled distance (Z), and are also available as polynomial equations. The pressure at a known point can be described by the modified Friedlander equation from Baker [8] and depends on the time of the arrival of the pressure wave (t = t 0 -ta).The other effects of nuclear explosions as thermal, electromagnetic pulse and radiation are disregarded. The peak over-pressure can be calculated using equation (1). This equation is valid for weapon yield (1kt or 1 Mt). Using scaling laws,the pressure can be determined for other yields. All parameters of the pressure time curve are normally written in terms of scaled laws utilizing equation (2). (1) , , (2) Where (W) is the weapon yield and (R)is the distance between the center of the charge and the object. The pressure attains its maximum very fast (extremely in short time called arrival time (t a ) of the shock wave to the point under consideration). Then, the pressure starts decreasing until it reaches the reference pressure (P 0 ) as shown in figure 1. The time of arrival (t a ) depends on the weapon yield and the distance from the point of burst. The time of arrival can be calculated using: (3) The positive phase duration (t d ) is the time for reaching the reference pressure. After that the pressure drops below the reference pressure until the maximum negative pressure (P min ) is reached. Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 Figure.1 Field pressure- time history (TM 5-855-1, 1986[9]) The negative phase duration is denoted as (t n ). The over-pressure impulse is the integral of the overpressure curve over the positive phase (t d ) and this can be computed using: (4) The peak value of the dynamic pressure (P ds ) as a function of peak pressure is: (5) In the current study, A MATLAB code has been developed to simulate the nuclear explosion by calculating the load time history as shown in figure 2.This program contains yield of the weapon in kt of TNT, type of blast (surface or air)and charge location. MATLAB calculations Explosive: TNT Quantity: 80Kt Figure 2: Pressure time history for 80 kt TNT at 800 m from MATLAB program Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 3 Theoretical Aspects of Design Optimization Let us mark O a family of admissible domainsΩ ∈ Rn. supposing that O is contained in some larger family Ō. With any Ω ∈O a Hilbert space V (Ω) of functions is associated supposing that a convergence of domains in Ō and of functions in V (Ω) is defined. If {Ωn⊂O}, Ω ∈Ō, then the symbol Ω n → Ω, as n → ∞ is used for the convergence of domains. By the symbol u n → u, as n → ∞ the convergence of function u n ∈ V (Ω n ) is detonated to the function u∈ Ω, where Ω n, Ω ∈Ō.. (6) (P) u : Ω ∈ O → u(Ω) ∈ V (Ω) Let be a mapping which associates with any domainΩ ∈ O a solution u(Ω) of a state problem (given by equations or inequalities in the domain ), and let G := {( Ω, u(Ω)) : Ω ∈ O }. Finally, let J : G → R be an objective (cost) function. The abstract optimal shape design problem is then stated as follows (P) (7) The exact solution of the problem is only in exceptional cases, so it is needed to approximate the problem (P), see [1], for example: Let h > 0 be a parameter of discretization tending to zero. In practice, any domain O h is determined by finite number of parameters, defining its boundary. With any O h a finite dimensional space V h (Ω h ) ⊂V (Ω h ) of functions defined on h is associated. The convergences are defined in the same way as in the continuous case. (8) Let be a mapping which associates with the domain (h), the solution U h (h) of the approximated state problem, and let G h := {( Ω h , u h (Ω h )) : Ω h ∈O h } be a graph of this mapping. The discrete optimal shape design problem can be stated as follows: (P h ) (9) Such problems need more computational time and requirements because of the evaluation of the cost functional or its gradient. It is needed to solve the state problem first, and another complication may bring the fact that the cost function is not necessarily convex and continuously differentiable. Moreover, the optimal solution may not be unique or exists. This is the reason why the detailed mathematical analysis of the shape design optimization problem is important. 4. Design Optimization in ANSYS The optimization routine termed OPT is an integral part of the ANSYS program that can be employed to determine the best solution. ANSYS also provides the user with a Parametric Design Language (APDL) that can be used to define design variables, objective function and constraint conditions. The ANSYS program offers two optimization methods to accommodate a wide range of optimization problems. These methods are the advanced zero–order method (Subproblem approximation method) and the first order method. For both of these methods, the program performs a series of analysis–evaluation–modification cycles. That is, an analysis of the initial design is performed, the results are evaluated against specified design criteria, and the Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 design is modified as necessary. This process is repeated until all specified criteria are met. ANSYS also provides helpful tools that can be used to obtain a good initial guess. 4.1. Advanced zero order method ANSYS uses optimization methods. These techniques are described in this section. In the advanced zero order method, the program uses only the values of objective function and constraints during the optimization (minimization) process. There are two concepts which and , and the conversion play a key role in this method: the use of approximations for of the constrained to unconstrained optimization problem using the penalty approach. The approximation is achieved by least square fitting between the data points. This results in curves (or surfaces) that are created in each optimization loop. The program allows the user to control curve fitting for the approximation. The user can request a linear fit, quadratic fit or quadratic plus cross terms fit. Using the advanced zero order method, the optimization procedure within ANSYS can be described in the following steps: Step 1. Define the design variables (design set) and initiate their values x o . These initial values can be either user defined or can be obtained using one of the program tools as described later. Step 2. Model the structural system (preprocessing stage). Step 3. Solve the structural system (solution stage) and obtain the values of the constraints (post– processing stage) and objective function F o . Step 4. Evaluate the feasibility of the design set between the current Gscu and the previous Gscu −1 using the feasibility tolerance ε sTol, G according to (10) Where ε sTol, G is either user–defined or the default is (11) Where GsU GsL are the upper and lower bounds of the constraint Gs . and Step 5. Check the convergence criteria at the end of each loop. The problem is said to be converged if the current , previous , or the best and any of the following conditions are satisfied. The change in the value F from F best to the current F cu is less than a certain tolerance ε Tol, F This can be formulated as (12) Where the value of ε Tol, F is either user–defined or by default: (13) The change in the value from F cu to F cu −1 is less than ε Tol, F , therefore (14) cu The changes in all design variables from x to x cu −1 , are less than ε iTol, x as (15) Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 Where ε iTol, x is user defined or the default: (16) U Where xi and xiL are user–defined upper and lower bounds of xi . The changes in all design variables from x cu to x best are less than their respective tolerances ε iTol, x . (17) Step 6. Check of termination conditions. The termination occurs when any of the following criteria is true: i. Any of the convergence criteria is reached. ii. The prescribed maximum number of loops is performed. iii. The number of consecutive infeasible solutions has reached the specified limit. Step 7. If the termination does not occur, then the software formulates an approximate sub– problem to determine best design vector ( x̂ ) by minimization of the objective function, based on current information. Using a series of response surfaces as a starting point for the next loop. This algorithm is called Sequential Unconstrained Minimization Technique SUMT (see Fiacoo and McCormick, 1968 [12]). The minimum of each response surface is found by a series of searches in design space, starting at the previous best design set. The search consists of minimization in the following directions: i. Directions along each x i , ii. Directions tangent to each G s ( x ) and iii. Direction of steepest descent of F ( x ) . This sequence continues until the change in the minimum of the response surface is less than a small tolerance that is described in Step 5. based on a combination of best design to date Step 8. Calculate a new design vector with the best design ( x̂ ) determined from last minimization. The new design variables are evaluated using the formula: (18) Where η is constant and is defined as: η = 1.0 − C o − C r ricu (19) o Where C is a fraction of the design variable in the next loop related to the current loop and it takes its value form: (20) 0 < C o < 0.9 , The random contribution fraction C r is determined by (21) 0 < C r < 1.0 − C o And the randomly generated number ricu applied to each xicu is − 0.5 < ricu < 0.5 Step 9. Repeat Steps 2–8 until the termination condition is achieved. (22) Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 4.2. First order method Using the same procedure illustrated for the advanced zero order method, the optimization procedure using the first order method can be described except: 1. In the first order method, the program uses the gradients ∂ F ( x ) ∂ x i and ∂ G s ( x ) ∂ x i . For each iteration, the gradient calculations are performed in order to determine a search direction, and a line search strategy is adopted to minimize the unconstrained problem. Thus, each iteration is composed of a number of sub–iterations that include search direction and gradient computations. 2. For the convergence checking, the program continues until either convergence is achieved or termination occurs. The problem is said to be converged if one of the following conditions is satisfied: • The change in the value F from F best to the current F cu is less than a certain tolerance ε Tol, F as described in the advanced zero order method. • The change in the value from F cu to F cu −1 is less than ε Tol, F as illustrated in the advanced zero order method. 5. Sizing Optimization of a Beam 5.1 Problems Statement In this section the mathematical model describing bending of an elastic beam with varying thickness is described. The concrete beam is made of homogeneous isotropic material with constant Young’s modulus of elasticity E. Body of the beam occupies the domain (D) (23) The blast load acts only in the y-axis direction. The concrete beam meets all requirements of Euler-Bernoulli hypothesis. The beam is supported along its entire length by a continuous elastic foundation with a stiffness function Q(x). Using Kantoroviˇc’s method the original displacement vector (u x , u y , u z ) is replace by its approximation (−ὠ (x)y, w(x), 0). The classical formulation of the beam bending problem in continuous case can be stated as follows: (4) (1) Find a deflection w C (Ω)∩ C ( )such that: (24) where f(x) is a continuous function in (0, L) and represents the intensity of the given external load dependent only on the coordinate x. Young’s modulus of elasticity of material E(x) and the moment of inertia of the cross section of the beam J(x) are continuous functions whose derivatives of J(x) up to the second order are continuous in the interval (0, L). Stiffness function of the elastic foundation (subsoil) of the beam Q(x) is a continuous function in (0, L). The bending moment M and the shearing force T corresponding to the deflection w = w(x) are defined as follows: (25) Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 (26) Some possible variants of given boundary conditions are introduced. Four types of classical boundary conditions may be described in points x = (0, L) (27) (28) The cross section of the beam is rectangular having breadth b(x) and thickness t(x) the beam is divided into ten elements. The state problem for fixed t(x) and Q(x) is obtained using: (29) Where J(x) = 1/12 t3(x)b(x)is the moment of inertia and K(x) = 1/12b(x)E(x) is a continuous function in (0, L). For simplicity K(x) equals constant value K; i.e. the beam is clamped at both ends. The definition of the space functions is: (30) Generally (Ω) ⊆V ⊆H2(Ω) holds. The following bilinear form and linear function can be obtained from the classical formulation using suitable Green’s formula. (31) 5.2. Classical continuous formulation For the optimization purposes, it is needed to formulate the set of admissible design variables. The design variables appear in coefficients of the corresponding differential operator, while the domain of integration remains fixed. Let us define the (32) The thickness is represented by functions which are continuous and uniformly bounded. The set of admissible thickness also preserves the beam volume. The condition |t\(x)| ≤ causes that t(x) is uniformly Lipschitz continuous. The uniform Lipschitz constraint appearing in the prevents thickness oscillations and is sufficient for definition of ad to be compact. The second set ad contains conditions for the second design variable Q where Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 (33) , , , Parameters t , , us define the design variable p(x) , and are chosen in such a way that U ad ≠0 ;. Now, let (34) It means that the optimal thickness distribution t(x) and the optimal distribution of the foundation stiffness Q(x) shall be found. By the cost functional a “quality” of the design variable is measured using the following function (35) The function has a quite practical meaning; J represents “stiffness” of the beam. For another cost function see [7], [8]. The corresponding sizing optimization problem can be then stated as follows: (P 1 ) (36) Where w(p(x), x) is a solution of the given state problem (P 1 ) for corresponding p(x)= (t(x),Q(x)). It can be proved that the optimization problem (P 1 ) has at least one solution, see [9]. 5.3. Approximation of the optimization problem and discrete formulation The exact solution of problem (P 1 ) is usually hardly available. Discrete formulation of the problem is the objective of this section. Let D n(h) is an equidistant division of [0, L] with step h>0 According to ANSYS possibilities the approximation of U ad is chosen by a piecewise linear respective piecewise constant functions on the same division D n(h) . First, sets and are defined as: (37) (38) The sets , and U ad can be approximated as follows: (39) Creation of a parametric model which will be directly implemented in finite element software ANSYS will be the next step in our analysis. From the approximation of the sizing optimization problem listed above now, to represent it in terms of real scalar parameters. The design variables t h and Q h with (n(h) + 1) are associated respective (n(h))- dimensional vectors. Components of these vectors can be obtained as function values at nodal points Dn(h) (40) (41) Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 Now the sets and SA-1 can be stated as follows: (42) where the relations hold for i = 0, . . . , n(h) respective i = 1, . . . , n(h) in the case of for i = 0, . . . , n(h) − 1 respective i = 1, . . . , n(h) − 1 in the case of and . Let us define Using the finite element method, the state problem takes form of the system of linear algebraic equations with positive definite stiffness matrix. It means that the state problem reduces for every fixed p h = (t h ,Q h ) to finding coefficients w i from the system K(P h )w(p h )= F. Now, the algebraic form of the cost function and the discrete form of the problem (P 1 ) are presented. (43) (P 1,h ) (44) Wherew(P h ) is the solution of the state problem K(p h ) w(p h ) = F for corresponding p h = (t h , Q h ). 6. Solution process In the finite element problem the concrete material properties and geometry properties are given in Table .1. Key points, element type and Young’s modulus of elasticity are defined in the analysis file. One set of real constants as AREA1, IZ1, AREA2, IZ2 for each element is defined. Table 1: Input data of the state problem Concrete material properties E= 3E10Pa ʋ = 0.2 Geometric properties I = 0.009365 m4 H = 0.762m B = 0.254m L = 6.096m Loading F=200kPa T d =650msec For the concrete beam of length (L) having fixed ends with blast load function P(t) . ANSYS element BEAM54 is chosen to model the beam element. Element (BEAM54) allows tapered cross section and constant stiffness of the foundation along the element. The thickness may vary between nodes.BEAM54 is conforming element for the current problem and corresponds to the choice of the space V h . In this case the supports, load and material properties are symmetrical ones. The number of variables is 6 due to symmetry. Each element (i) has an initial thickness (Ti) Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 ParameterT1 represents the thickness t 1 in the first node and along the first element. In the same way the remaining parameters T2, T3, . . . , T6 are defined. For each element (i) parameters express the constraints of the derivative and parameters of the type (T2/T1)3 are defined to prevent the computational error caused by the element BEAM54 assumption and restriction. The solution phase begins by the definition of displacement boundary conditions and blast load. Half symmetry boundary conditions have been prescribed at node number 6.The optimization file contains commands which specify the parameters as the design variables, state variables , the objective function and the maximum and minimum values of these variables are defined as shown in Table 2. Table.2. Typical values of input parameters used in the input file Parameter Value Design set Parameter Min Max Tolerance Design variable Thickness of beam 0.1905m 0.381m 0.01 Volume of beam 1.78 m3 1.9 m3 0.1 Derivative of the thickness -1 1 0.01 State variable Plastic strain limit Objective 0.15 Displacement Minimum 0.00001 1st order method, Sub-problem approximation Optimization method Remaining parameters are set in the same way. The convergence tolerance was set to 10−5. The value of the cost function J for the initial constant thickness t i = 0.25, ฀i. In Tab.3 and Tab. 4. The first order method was run with variable input parameters. The percent step size was set to 100 each time. The convergence tolerance is set to 10−5 and with the percent forward difference ∆D = 0.1 .The best minimum displacement value 11.019 mm is reached. In the second run ∆D = 0.2 is set. The run was terminated after 3 loops because the convergence conditions were satisfied. The convergence tolerance was changed to 10−6 but better solution than in the first run was not reached. The same situation occurred when the forward difference is changed to 0.3 or 0.4, see Table. 3. Table 3: First Order Method – Results Trial number Input data Tolerance ∆D Number of iter Minimum displacement(mm) 1 10−5 0.2 2 16.023 2 10−6 0.2 125 11.959 3 10−5 0.1 87 11.019 −6 4 10 0.3 71 14.202 5 10−6 0.4 127 13.203 Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 The sub problem approximation method was run each time with the tolerance parameter set to 10−5. After several runs of this method, it was clear that the values of input parameters are very important according to the final result. Especially the maximum number of iterations (MNI) and maximum number of infeasible sets (MNIS). When running the method with MNI=140 and MNIS=60. The best result with minimum displacement value 15.816 mm is reached after every possible combination of other input parameters. The process was always terminated because of the maximum number of iterations was reached. Value 15.816 is much worse than the best result obtained by the first order method. After some experiments MNI=300 and MNIS=300 are set. The method reached its best sets somewhere between 200 and 320 iterations, see Table 4. Table 4: Sub problem Approximation Method – Results Trial number 1 2 3 4 5 Tol. 10−5 10−5 10−5 10−5 10−5 OBJ approx. Linear Quadratic Quadratic Quadratic Quadratic Input data CON app. Linear Quadratic Quadratic Quadratic Quadratic Weights Compound Compound Feasibility Objective Distance Iter. 300 300 265 254 300 Minimum displacement(mm) 10.563 11.312 10.462 11.389 10.494 By setting of these values better results with this method are obtained than with the first order method. From Tab. 2 and Tab. 3, it can be seen that the best results obtained by the first order method and the sub problem approximation method can be compared in figure 3 and in figure 4. The best result for the first order method is presented on the left figure. Figure 3: The optimal thickness Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 Figure 4: Graphs of functional values for both methods 7. Conclusions The main aim of this paper is the presentation of the procedure suitable for shape design optimization of concrete beam subjected to blast loads. This problem has wide usage in technical and engineering applications, including fortifications and bridges. Results obtained by two different optimization methods are summarized into tables and analyzed. The first order method, compared to the sub problem approximation method, usually demands more computations and is more accurate. The first order method proceeds systematically to the optimum. On the other hand, the sub problem approximation method is generally much faster. It is clear that both methods are dependent on the choice of their parameters and tolerances. 8. References [1] Dharaneepathy, M.V. and Sudhesh, K.G. Optimal Stiffening of square plates subjected to airblast loading. Compuers& Structures.Vol 36 No 5, 891-899 (1990) [2]. Xue, Z. and Hutchinson, J.W., A Comparative study of impulse-resistant metal sandwich plates, International Journal of Impact Engineering, 30, 1283–1305 (2004) [3]. Fleck, N.A. and Deshpande, V.S., The Resistance of Clamped Sandwich Beams to Shock Loading, vol. 71, J. of Applied Mechanics, Transactions of the ASME, 386-401, (2004) [4]Yen, C.F., Skaggs, R., Cheeseman, B.A., Modeling of shock mitigation sandwich structures for blast protection. The 3rd First International Conference on structural stability and dynamics, Kissimmee, Florida, June 19-22 (2005) [5]. Liang, C.C.,Yang, M.F., Wu, P.W ,Optimum design of metallic corrugated core sandwich panels subjected to blast loads. Ocean Engineering, 28, 825–861 (2001) [6] Main, J.A. and Gazonas, G.A, Uniaxial crushing of sandwich plates under air blast: Influence of mass distribution, International Journal of Solids and Structures. 45, 2297–2321 (2008) [7] Kingery, Charles N.: Bulmash, Gerald: Airblast Parameters from TNT Spherical Air Burst and Hemispherical Surface Burst, Defense Technical Information Center, Ballistic Research Laboratory, Aberdeen Proving Ground, Maryland, 1984. Proceedings of the 10th ICCAE-10 Conference, 27-29 May, 2014 SA-1 [8] Baker, Wilfrid E.: Explosions in the Air, University of Texas Pr., Ausint, 1973. [9]TM5,”Structures to Resist The Effects of Accidental Explosion”, USA (1990).. [10] Release 11.0 Documentation for ANSYS [11] Moaveni, S. : Finite Element Analysis: Theory and Applications with ANSYS, Prentice Hall, New Jersey, 1999. [12] Fiacoo and McCormick, 1968 Nonlinear Programming: Sequential Unconstrained Optimization Techniques ,Willy , New York.