C O M P U T A T I O N A L M E C H A N I C S IN STRUCTURAL ENGINEERING
Recent Developments
Elsevier Science Internet Homepage http://www.elsevier.nl (Europe) http://www.elsevier.com (America) http://www.elsevier.co.jp (Asia) Consult the Elsevier homepage for full catalogue information on all books, journals and electronic products and services.
Elsevier Titles of Related Interest CHENG & SHEU Urban Disaster Mitigation The Role of Engineering and Technology ISBN. 0080419208
GODOY ThinWalled Structures with Structural Imperfections: Analysis and Behavior ISBN. 008042266 7
CHENG & WANG PostEarthquake Rehabilitation and Reconstruction ISBN: 0080428258
FUKUMOTO Structural Stability Design ISBN. 0080422632
SRIVASTAVA Structural Engineering World Wide 1998 (CDROM Proceedings with Printed Abstracts Volume, 702 papers) ISBN. 0080428452 OWENS Steel in Construction (CDROM Proceedings with Printed Abstracts Volume, 268 papers) ISBN." 0080429971
USAMI & ITOH Stability and Ductility of Steel Structures ISBN. 0080433200 GUEDESSOARES Advances in Safety and Reliability (3 Volume Set) ISBN." 0080428355 FRANGOPOL, COROTIS & RACKWITZ Reliability and Optimization of Structural Systems ISBN." 0080428266
Related Journals Free specimen copy gladly sent on request." Elsevier Science Ltd, The Boulevard, Langford Lane, Kidlington, Oxford, 0)(5 1GB, U.K.
• • • • •
Advances in Engineering Software Composite Structures Computers and Structures Construction and Building Materials Computer Methods in Applied Mechanics and Engineering • Engineering Analysis with Boundary Elements
• • • • • • •
Engineering Failure Analysis Engineering Structures Journal of Constructional Steel Research Probabilistic Engineering Mechanics Reliability Engineering and System Safety Structural Safety ThinWalled Structures
COMPUTATIONAL MECHANICS IN STRUCTURAL ENGINEERING Recent Developments
Edited by
Franklin Y. Cheng
UniversiO' of MissouriRolla, USA and
Yuanxian Gu
Dalian University of Technology, China
1999
ELSEVIER
An Imprint of Elsevier Science Amsterdam • Lausanne • New York • Oxford • Shannon • Singapore • Tokyo
ELSEVIER SCIENCE Ltd The Boulevard, Langford Lane Kidlington, Oxford OX5 1GB, UK
© 1999 Elsevier Science Ltd. All rights reserved.
This work is protected under copyright by Elsevier Science, and the following terms and conditions apply to its use: Photocopying Single photocopies of single chapters may be made for personal use as allowed by national copyright laws. Permission of the publisher and payment of a fee is required for all other photocopying, including multiple or systematic copying, copying for advertising or promotional purposes, resale, and all forms of document delivery. Special rates are available for educational institutions that wish to make photocopies for nonprofit educational classroom use. Permissions may be sought directly from Elsevier Science Rights & Permissions Department, PO Box 800, Oxford OX5 1DX, UK; phone: (+44) 1865 843830, fax: (+44) 1865 853333, email:
[email protected] You may also contact Rights & Permissions directly through Elsevier's home page (http://www.elsevier.nl), selecting first 'Customer Support', then 'General Information', then 'Permissions Query Form'. In the USA, users may clear permissions and make payments through the Copyright Clearance Center, Inc., 222 Rosewood Drive, Danvers, MA 01923, USA; phone: (978) 7508400, fax: (978) 7504744, and in the UK through the Copyright Licensing Agency Rapid Clearance Service (CLARCS), 90 Tottenham Court Road, London W1P 0LP, UK; phonei (+44) 171 436 5931; fax: (+44) 171 436 3986. Other countries may have a local reprographic rights agency for payments. Derivative Works Tables of contents may be reproduced for internal circulation, but permission of Elsevier Science is required for external resale or distribution of such material. Permission of the publisher is required for all other derivative works, including compilations and translations. Electronic Storage or Usage Permission of the publisher is required to store or use electronically any material contained in this work, including any chapter or part of a chapter. Contact the publisher at the address indicated. Except as outlined above, no part of this work may be reproduced, stored in a retrieval system or transmitted in any form or by any means, electronic, mechanical, photocopying, recording or otherwise, without prior written permission of the publisher. Address permissions requests to: Elsevier Science Rights & Permissions Department, at the mail, fax and email addresses noted above. Notice No responsibility is assumed by the Publisher for any injury and/or damage to persons or property as a matter of products liability, negligence or otherwise, or from any use or operation of any methods, products, instructions or ideas contained in the material herein. Because of rapid advances in the medical sciences, in particular, independent verification of diagnoses and drug dosages should be made.
First edition 1999 Library of Congress Cataloging in Publication Data Computational mechanics developments / edited ed.
In s t r u c t u r a l engineering : recent by F r a n k l i n Y. Cheng and Y u a n x l a n Gu.

1st
p. cm. ISBN 0080430082 (hc.) I. Structural ana]ysls (Engineering)Data processlng. 2. Computeralded englneerlng. I. Cheng, Franklin Y. II. Gu, Yuanxian. TA647.C64 1999 9851822 624.1'7dc21 CIP
ISBN: 0080430082
~ ) The paper used in this publication meets the requirements of ANSI/NISO Z39.481992 (Permanence of Paper). Printed in The Netherlands.
CONTENTS
Preface
ixx
Acknowledgements Symposium Activities Resolutions
xiiixviii xix=xx
KEYNOTE PRESENTATIONS Multiobjective Optimum Design of Structures with Genetic Algorithm and Game Theory: Application to LifeCycle Cost Design Franklin Y. Cheng New Solution System for Plate Bending W.X. Zhong and W.A. Yao
17
Finite Element Algorithm Based on TwoScale Analysis Method J.Z. Cui
31
Control of LateralTorsional Motion of Naming TV Transmission Tower J.C. Wu and J.N. Yang
43
PAPER PRESENTATIONS Fast Computation of Stationary/NonStationary Random Responses of Complex Structures J.H. Lin and W.X. Zhong
57
Response of Dynamical Systems Driven by Additive Gaussian and Poisson White Noises L.A. Bergman, S.F. Wojtkiewicz and M. Grigoriu
71
Solving Large Systems of Equations on IntelParagon P. Chen, P. Tong, T.Y.P. Chang and S.L. Sun
85
Dynamic Behavior of Railway Bridges under Random Loading and Assessment of VehicleRunning Safety YingJun Chen, He Xia and Daqing Wang
97
A Review on the Numerical Solution Schemes for Localization Problems Z. Chen, X. )(in and D. Qian
111
EContinuation Approach for Truss Topology Optimization J(u Guo and Gengdong Cheng
125
Recent Developments in Basic Finite Element Technologies Carlos A. Felippa
141
vi
Contents
Structural Optimization for Practical Engineering: Software Development and Applications Yuanxian Gu, Hongwu Zhang, Zhan Kang and Zhenqun Guan Parallel Integration Algorithms for Dynamic Analysis of Structures in the Clustered Network System J.G. Cheng, Z.H. Yao, X.P. Zheng, Y. Gao, Z.J. Kou and W.B. Huang
157
169
Maximum Entropy Principle and Topological Optimization of Truss Structures B. Y. Duan, Y.Z. Zhao and H. Liu
179
Practical Issues in the Application of Structural Identification Burcu Gunes, Behnam Arya, Sara WadiaFascetti and Masoud Sanayei
193
ER Devices for Control of Seismically Excited Structures Henri P. Gavin
207
Parallel PCG Algorithm on Distributed Network by PVM Z.C. Hou, J.G. Cheng, Z.H. Yao and Z.C. Zheng
221
Computer Aided Design for Vibration Isolation Systems with Damped Elastic Stops H.Y. Hu and F.,V. Wang
231
Calculation of Thin Plate on Statistical NonUniform Foundations Yi Huang, Yuming Men and Guansheng Yin
245
Computer Simulation of Structural Analysis in Civil Engineering Jiang JianJing, Guo WenJun and Hua Bin
259
A Mixed Finite Element for Local and Nonlocal Plasticity Xikui Li
267
Behavior of InFilled Steel Plate Panels Subjected to Cyclic Shear Ming Xue and LeWu Lu
281
A Finite Element Model for Geometrically Nonlinear Analysis of MultiLayered Composite Shells Ari[" Masud and Choon L. Tham
295
Nonlinear and Buckling Analysis of Complex Branched Shells of Revolution J.G. Teng
309
Structural Optimization for Seismic Loads: PseudoStatic, Response Spectra and Time History Kevin Z. Truman
325
Monitoring of Cable Forces Using MagnetoElastic Sensors M.L. Wang, D. Satpathi, S. Koontz, A. Jarosevic and M. Chandoga
337
Fracture Estimation: Bound Theorem and Numerical Strategy ChangChun Wu, QiZhi Xiao and ZiRan Li
349
Finite ElementBased Buffeting Analysis of Long Span Bridges Y.L. Xu, D.K. Sun, J.M. Ko and J.H. Lin
361
Contents
ODEOriented SemiAnalytic Methods
vii
375
Si Yuan
Index of Contributors
389
Keyword Index
391
This Page Intentionally Left Blank
PREFACE
Since the First SinoUS Joint Symposium/Workshop on Recent Developments and Future Trends of Computational Mechanics in Structural Engineering was held September 2428, 1991, in Beijing, both countries have been actively engaged in cooperative research benefitting both countries and the research communities at large. Researchers from the US and China have a strong desire and commitment to exchange stateoftheart information and to advance scientific progress. The Second SinoUS Symposium Workshop on Recent Advancement of Computational Mechanics in Structural Engineering was held May 2528, 1998, in Dalian, China. The objectives were: to share the insights and experiences gained from recent developments in theory and practice to assess the current state of knowledge in various topic areas of mechanics and computational methods and to identify joint research opportunities •
to stimulate future cooperative research and to develop joint efforts in subjects of common needs and interests to build and to strengthen the longterm bilateral scientific relationship between academic and professional practicing communities
Topic areas discussed cover the entire field of computational structural mechanics. These topics have advanced broad applications in the engineering practice of modern structural analysis, design and construction of buildings and other structures, and in natural hazard mitigation. Speakers were invited through steering committees in both countries. The Chinese steering committee consisted of: Professor Wanxie Zhong (chair)
Department of Engineering Mechanics, Dalian University of Technology Professor Junzhi Cui
Institute of Computational Mathematics and ScienceEngineering Computation, Academia Sinica Professor Yuanxian Gu
Department of Engineering Mechanics, Dalian University of Technology Professor Jianjing Jiang
Department of Engineering Mechanics, Tsinghua University
x
Preface
T h e U S steering c o m m i t t e e c o n s i s t e d of: Dr. J a n n N. Y a n g (chair)
Professor of Civil Engineering, University of California at Irvine Dr. F r a n k l i n Y. C h e n g
Curators' Professor of Civil Engineering University of MissouriRolla Dr. Z h e n C h e n
Assistant Professor of Civil Engineering, University of MissouriColumbia C h i n e s e d e l e g a t e s included: Professor Professor Professor Professor Professor Professor Professor Professor Professor
P u Chen, . . . . . . . . . . . . . Y i n g j u n Chen, . . . . . . . . . Gengdong Cheng, . . . . . . J u n z h i Cui, . . . . . . . . . . . B a o y a n Duan, . . . . . . . . . Y u a n x i a n Gu, . . . . . . . . . Z h i c h a o Hou, . ........ H a i y a n Hu, . . . . . . . . . . . Yi H u a n g ,
Professor Professor Professor Professor Professor Professor Professor Professor Professor
J i a n j i n g Jiang, . . . . . . . . . X i k u i Li, . . . . . . . . . . . . . J i a h o a Lin, . . . . . . . . . . . . J.G. Teng, . . . . . . . . . . . . C h a n g c h u n Wu, . ...... Y o u l i n Xu, . . . . . . . . . . . Z h e n h a n Yao, . . . . . . . . . Si Y u a n , . ............ Wanxie Zhong, . .......
Beijing University Northern Jiaotong University Dalian University of Technology Chinese Academy of Sciences Xian Universityof Electronic Science & Technology Dalian University of Technology Tsinghua University Nanjing University of Aeronautics & Astronautics Xian University of Architectural Science & Technology Tsinghua University Dalian University of Technology Dalian University of Technology Hong Kong Polytechnic University University of Science & Technology of China Hong Kong Polytechnic University Tsinghua University Tsinghua University Dalian University of Technology
US delegates included: Professor Professor Professor Professor Professor Professor Professor Professor Professor Professor Professor
L a w r e n c e A. B e r g m a n , .. Z h e n Chen, . . . . . . . . . . . F r a n k l i n Y. C h e n g , . . . . . C a r l o s A. Felippa, . . . . . . H e n r i P. Gavin, . . . . . . . . L e  W u Lu, . . . . . . . . . . . . Arif Masud, . .......... K e v i n Z. T r u m a n , . . . . . . Sara W a d i a  F a s c e t t i , . . . . M i n g L. W a n g , . . . . . . . . J a n n N. Y a n g . . . . . . . . . .
University of Illinois at UrbanaChampaign University of MissouriColumbia University of MissouriRolla University of Colorado at Boulder Duke University Lehigh University University of Illinois at Chicago Washington University Northeastern University University of Illinois at Chicago University of California at Irvine
All technical papers were carefully r e v i e w e d by the steering c o m m i t t e e s . P r o f e s s o r F r a n k l i n Y. C h e n g c o m p l e t e d the final editing o f the papers in the p r o c e e d i n g ' s v o l u m e for p u b l i c a t i o n . T h e i r c o n t r i b u t i o n s are s p e c i a l l y a c k n o w l e d g e d . S.C. L I U , P r o g r a m D i r e c t o r K.P. C H O N G , P r o g r a m D i r e c t o r Civil & M e c h a n i c a l S y s t e m s National Science Foundation
ACKNOWLEDGMENTS
The Second SinoUS Symposium/Workshop on Recent Advancement of Computational Mechanics in Structural Engineering was jointly sponsored by the US National Science Foundation and the Chinese National Natural Science Foundation and Ministry of Education. As Editors of the proceedings, we gratefully acknowledge Drs. S.C. Liu and K.P. Chong of the National Science Foundation and Professors L.X. Qian and W.X. Zhong of Dalian University of Technology for their encouragement and advice. Papers were reviewed by both steering committees as well as Professors Lawrence Bergman, Zhen Chen, Carlos Felippa, Henri Gavin, Arif Masud, Kevin Truman, Sara WadiaFascetti, and Ming Wang. Workshop resolutions were integrated by Professor Bergman. We thank them all for their cooperation and assistance. Our gratitude goes to Dr. James Milne, Senior Publishing Editor, Engineering and Technology, Elsevier Science Ltd, Oxford, England for his guidance and to Ms. Elizabeth Farrell at the University of MissouriRolla for her technical support. Papers assembled in the proceedings volume are intended to serve as a reference tool and as supplementary course material for educators, researchers and practitioners in the international engineering community.
FRANKLIN Y. CHENG
Curators' Professor of Civil Engineering, University of MissouriRolla, USA Honorary Professor, Harbin University of Architecture & Engineering, China Honorary Professor, Xian University of Architecture & Technology, China Honorary Professor, Yunnan Polytechnic University, China YUANXIAN GU
Professor and Head, Department of Engineering Mechanics Research Institute of Engineering Mechanics Dalian University of Technology, China
This Page Intentionally Left Blank
Second SinoUS Joint Symposium on Recent Advancement of Computational Mechanics in Structural Engineering 2528 May 1998, Dalian, China
Front row: Back row:
K. Fu, J.H. Lin, X.S. Xu, M.L. Wang, P. Chen, Z.C. Hou, Y.X. Gu, J.G. Teng, Y.L. Xu R.F. Wu, Q.G. Meng, C.C. Wu, Y.J. Chen, H.Y. Hu, G.D. Cheng, W.M. Dong, A. Masud, Z.H. Yao, X.K. Li, W.X. Zhong, K.Z. Truman, C.A. Felippa, S.C. Liu, L.W. Lu, F.Y. Cheng, L.A. Bergman, J.Z. Cui, L.X. Qian, S. WadiaFascetti, J.J. Jiang, H. Gavin, J.N. Yang, Z. Chen, B.Y. Duan, S. Yuan
X. .
xiv
'Symposium Activities'
Opening remarks by G.D. Cheng, President, Dalian University of Technology, with speakers (left to right) L.X. Qian, W.X. Zhong, S.C. Liu, F.Y. Cheng, J.N. Yang, C.A. Felippa, Q.G. Meng
Conference Room
'Symposium Activities'
xv
Welcoming remarks by Junwen Liu, Deputy Mayor of Dalian (left to right) L.X. Qian (partially hidden), W.X. Zhong, S.C. Liu, Mayor Liu, F.Y. Cheng, J.N. Yang, C.A. Felippa, Q.G. Meng
........ .............. ::..~.:: ..: ~::.;~:::.::.::.:~:,;,:.~.:;~::~,:,:,::,:,~.:.~:
=========================================== ~,:: ~,:.: :::,~::::~ :, ?:: :~:::::::~:::~....
i!!ii! ;i::i!
'
.s:ii::i ~
iI
G.D. Cheng and F.Y. Cheng outside conference building
xvi
'Symposium Activities'
US delegates (left to right) S. WadiaFascetti, M.L. Wang, L.W. Lu, J.N. Yang, L.A. Bergman, H.P. Gavin, K.Z. Truman, C.A. Felippa, F.Y. Cheng, Z. Chen, A. Masud
Delegates at reception
'Symposium Activities'
Luncheon for delegates (clockwise) K.Z. Truman, L.A. Bergman, S. WadiaFascetti, J.Z. Cui, W.H. Zhong, M.L. Wang, F.Y. Cheng, A. Masud, C.A. Felippa
Longtime friends (seated) X.S. Li, F.Y. Cheng, W.X. Zhong, (standing) J.H. Lin
xvii
xviii
'Symposium Activities'
W.M. Dong and S.C. Liu en route to site visit
W.A. Anderson, NSF Senior Advisor, and F.Y. Cheng in front of hotel
R E S O L U T I O N S OF S E C O N D SINOUS JOINT S Y M P O S I U M / W O R K S H O P ON R E C E N T A D V A N C E M E N T OF C O M P U T A T I O N A L M E C H A N I C S IN S T R U R C T U R A L E N G I N E E R I N G
The second SinoUS Joint Symposium/Workshop on Recent Advancement of Computational Mechanics in Structural Engineering, jointly chaired by Professor Wanxie Zhong of the Dalian University of Technology and Professor Jann N. Yang of the University of California at Irvine, was held May 2528, 1998, in Dalian, China. Presentations and discussions during the symposium mainly focused on the areas of structural optimization, stochastic structural dynamics, structural control and system identification, linear and nonlinear finite element analysis, largescale and parallel computing, and structural design. At the conclusion of the workshop, which followed the symposium, the participants unanimously passed the following resolutions: Encourage collaborative research projects between Chinese and US individual investigators and groups of investigators, through shortterm exchange of science and engineers; Seek NSF (US) and NNSF (China) joint research support to promote and facilitate the evaluation of joint research proposals; Establish future SinoUS joint symposia/workshops on a biennial basis, with venues alternating between the US and China; Propose that the Third SinoUS Joint Symposium/Workshop be scheduled for the year 2000 and that the theme of the symposium be "Methods and Applications of Optimization in Design" with St. Louis, Missouri as a possible location. Professor Kevin Truman of Washington University in St. Louis will investigate such possibilities and will coordinate his efforts with Professor Yuanxian Gu of the Dalian University of Technology. Resolutions thanking the cochairs for their efforts and the local organizing committee for their outstanding arrangements and hospitality also passed unanimously. Topics of mutual interest that were recommended for further investigation as potential joint research activities by the three working groups are as follows.
Working Group 1: Optimization of Civil Infrastructure Systems Structural Design Using Optimization Theory PerformanceBased Design  application of optimization theory to performancebased design including objectives and constraints that reflect reliability, serviceability, and strength. LifeCycle Costs  development of optimizationbased system (integrated) design including
xix
xx
Resolutions
elements such as multiple objectives, damage costs, strength degradation, and environmental loads. Upgrade, Repair and Retrofit  development of advanced optimization theory to assess structural enhancements (designs) for improved safety, risk mitigation, and potential loss reduction. Nonlinear Systems  development of advanced optimization theory for application to nonlinear systems due to conditions such as material nonlinearity, large deformations, environmental loads, dynamic loads, optimal control, deterioration or fracture. SocioEconomic Impact integration of socioeconomic factors into the optimal structural design theory based on a given structure's impact on elements such as commerce, transportation networks, and community living. US/China Joint Meeting  exchange of practical applications of structural optimization through a joint meeting between the ASCE committee on Optimal Design, other qualified US researchers, and their Chinese counterparts. HighPerformance Computing  application of highperformance computing for the optimal design of largescale systems.
Working Group II: New Technologies in Structural Engineering HighPerformance Structures  objectivebased issues include durability, damage control, and survival. Random VibrationBased Design  development of computational methods for random vibration analysis, statistical identification of loading, and development of probabilitybased design criteria. Advanced Monitoring Systems for Large Structures  new sensor technologies, algorithm development, damage identification, and prediction models for lifecycle design. Advanced Materials  homogenization for composite materials, engineeringbased design of advanced materials, and design of structural components (i.e., connections). Passive/Hybrid/Active/SemiActive Technologies for Structural implementation, theoretical developments, and device development.
Control

practical
Working Group III: Mathematical and Computational Aspects of Design and Optimization LongSpan Bridges  advanced methods for analysis and design under multiple loadings (earthquake, wind, traffic). Integrated Structural and Material Design development of unified approaches. HighPerformance Computing  application of highperformance computing and communication, and virtual reality modeling and visualization, to the analysis and engineering design of largescale structures. Multiphysics  development of a computational testbed for multiphysics, multiscale, material modeling and testing. Symbolic Computation integration of symbolic computation into largescale multidisciplinary problems in infrastructural engineering.
M U L T I O B J E C T I V E O P T I M U M DESIGN OF S T R U C T U R E S WITH GENETIC A L G O R I T H M AND G A M E THEORY: APPLICATION TO L I F E  C Y C L E COST DESIGN Franklin Y. Cheng Curators' Professor, Department of Civil Engineering Senior Investigator, Intelligent Systems Center University of MissouriRolla Rolla, MO 654090030 USA
ABSTRACT Loss of life and property from possible future earthquakes as well as the expense and difficulty of postearthquake rehabilitation and reconstruction strongly suggest the need for proper structural design with damage control. Design criteria should balance initial cost of the structure with expected losses from potential earthquakeinduced structural damage. Lifecycle cost design addresses these issues. Such a design methodology can be developed using multiobjective and multilevel optimization techniques. Presentation here focuses on genetic algorithm and game theory as well as a lifecycle cost model for this innovative design methodology. Genetic algorithms (GAs) have the characteristic of maintaining a population of solutions, and can search in a parallel manner for many nondominated solutions. These features coincide with the requirement of seeking a Pareto optimal set in a multiobjective optimization problem. The rationale for multiobjective optimization via GAs is that at each generation, the fitness of each individual is defined according to its nondominated property. Since nondominated individuals are assigned the highest fitness values, the convergence of a population will go to the nondominated zone  the Pareto optimal set. Based on this concept, a Pareto GA whose goal is to locate the Pareto optimal set of a multiobjective optimization problem is developed. In this GA, to avoid missing Pareto optimal points during evolutionary processes, a new concept called Paretoset filter is adopted. At each
2
Cheng, F.Y.
generation, the points of rank 1 are put into the filter and undergo a nondominated check. In addition, a niche technique is provided to prevent genetic drift in population evolution. This technique sets a replacement rule for reproduction procedures. For a constrained optimization problem, a revised penalty function method is introduced to transfer a constrained problem into a nonconstrained one. The transferred function of a point contains information on a point's status (feasible or infeasible), position in a search region, and distance to the Pareto optimal set. Three numerical examples are provided: 1) optimum design of a seismicresistant structure with/without control, 2) optimum design for a final structural system selected from steel frame, reinforced concrete, or composite system, and 3) sensitivity analysis of the effect of cost function on structural probability failure. It is concluded that multiobjective and multilevel optimization is essential to determine target reliability and seismic building code performance.
KEYWORDS Genetic algorithm, game theory, lifecycle cost, multiobjective and multilevel optimization, earthquake, probability failure, fuzzy logic, Pareto set filter, niche technique, control.
INTRODUCTION In the current engineering design community, major design efforts are based on a conventional trial and error approach for which the relative stiffness of a structure's constituent members must be assumed. If preliminary stiffness is misjudged, then repeat analysis, even with a sophisticated computer program, will usually not yield an improved design. The optimum design concept is recognized as being more rational and reliable than the conventional design approach. Considerable literature has been published on the subject of optimal structural design for singleobjective function (Cheng and Truman, 1985; Cheng [ed.], 1986; Cheng and Juang, 1988; Frangopol and Cheng [eds.], 1996) Most realworld design optimization problems in structures are multimodal. There often exist several objectives to be considered by the designer. Usually these objectives are conflicting rather than complementary. A singleobjective optimization formulation does well with respect to an optimal objective, but the design may not always be a "good" design. Consider a hypothetical example. If a structure is optimized for minimum weight subject to constraints such as stress, displacement, buckling and vibration period, a structure is then obtained with minimum constructed materials. However, the structure may have a poor performance of dynamic response under the action of seismic loadings. If the minimum earthquake input energy is also included as objective, a more rational, compromise design will be produced (Cheng and Li, 1996; Cheng and Li, X.S., 1998)). This combined formulation is a multiobjective optimization problem (MOP). Multiobjective optimization offers the possibility to consider effectively all the different, mutually conflicting requirements inherent in a design problem.
Multiobjective Optimum Design with Genetic Algorithm and Game Theory
3
In game theory, if players agree to cooperate, a Pareto optimum will be an ideal solution because it has the property that if any other solution is used, at least one player's performance index is worse, or all the players do the same. This study demonstrates how game theory as a design tool applies to an MOP, and describes the relationship between cooperative game theory and Pareto optimal solution. Three genetic algorithms for multiobjective optimization are proposed based on game theory. In the Pareto GA, whose goal is to find a representative sampling of solutions along with the Pareto optimal set, two new techniques are investigated: a new operator called Paretoset filter is introduced to prevent the loss of Pareto optimum points in evolutionary progress; and niche technique is created by putting limitations on reproduction operators. Pareto GA for a constrained MOP is further studied te include fuzzylogic scheme. Lifecycle cost model is introduced along with multilevel optimiziation concept. The proposed multiobjective optimization techniques are applied to the optimum design of a seismic structure with/without control and applied to evaluate a structural system as wether it should be steel frame, reinforced concrete frame, or composite steelandreinforcedconcrete flame. Numerical results show that multiobjective optimization is essential to produce a good seismic structural design.
MULTIOBJECTIVE OPTIMIZATION AND PARETO OPTIMUM
Multiobjective optimization can be defined as determining a vector of design variables that are within the feasible region to minimize (maximize) a vector of objective functions and can be mathematically expressed as follows Minimize
F(x)  {fl(x), fz(X)...... fm(X)}
S u b j e c t to
g(x) oo end boundary conditions. Now the end boundary condition of x = 0 is applied to determine the constants f , , (n = 1,2,...). Only the first k terms of Eqn. 70 are selected to solve, the variational equation of x  0 end boundary reads (71) Due to some of eigenvalues are complex numbers, the real valued canonical equation of Eqn.7071 (Zhong, 1995) is applied. As an example, let v = 0.3, k  11 and k = 21 respectively. The distributions of bending moment at the clamped end are shown in Figure 2. It is clear that the comer has stress singularity ( M x >oo ). The wavy form of bending moment appears due to choosing only finite number of terms in expansion solution. 1.4 1.2 1 0.8 0.6 0.4 0.2 0
o. 2 0
o'2
o'4
  k : 2 1
Figure 2:
o;6
o:8
1
....... k:11
Distribution of bending moment at clamped end (x = 0 )
The method proposed in this paper can also be applied to other boundary conditions. A series of analytic solution can be derived.
30
Zhong, W.X. and Yao, W.A.
CONCLUDING REMARKS In this paper, a new formulation of fundamental equations and the respective solution methodology for classical theory of plate bending is presented, which contrasts with the traditional methodology (Timoshenko and WoinowskyKrieger,1959). The traditional solution uses displacement method and biharmonic equation; the new methodology applies bending moment function vector and differential Eqn.15. The traditional method is semiinverse solution; thus only a few of solutions available. For example, only solution for both lateral simply support plate have been given for rectangle plate, and it is difficult to solve for other lateral boundary conditions. However, the new methodology presents a direct solution via introducing Hamiltonian system. Therefore the efficient mathematical methods, for example separation of variables and eigenfunction expansion etc., can be applied. The new methodology breaks the limitation of traditional semiinverse solution. The solution for both sides free plate given in this paper can't be solved by the traditional method. The analogy between plate bending and plane elasticity can be applied to not only analytical solutions but also to plate bending finite element; therefore the plate bending finite element can be improved to the same level as plane elasticity (Zhong and Yao, 1998). Generally speaking, the new advances in plate bending mathematical theory present many opportunities. More research is anticipated. ACKNOWLEDGMENTS
This work was supported by the National Natural Science Foundation of China. REFERENCES Hu H.C. (1981). Variational Principle of Elasticity and Its Application, Science Press, Beijing, P.R. China (in Chinese) Muschelishvili N.I. (1953). Some Basic Problems of the Mathematical Theory of Elasticity. P. Noordhoff, Groningen, Netherlands Timoshenko S.P. and WoinowskyKrieger S. (1959). Theory of Plates and Shells, McGrawHill, New York, USA Timoshenko S.P. and Goodier J.N. (1970). Theory of Elasticity, 3rd Ed. McGrawHill, New York, USA Zhong W.X. (1995). A New Systematic Methodology for Theory of Elasticity, Dalian University of Technology Press. Dalian, P.R. China (in Chinese) Zhong W.X. and Yang Z.S.(1992). Partial Differential Equations and Hamiltonian System, Computational Mechanics in Structural Engineering, Elsevier. SINO/US Joint Symposium on Computational Mechanics in Structural Engineering, Sept. 1991, Beijing, Zhong W.X. and Yao W.A. (1998). Similarity Between Finite Element in Plane Elasticity and Plate Bending, Chinese J. of Computational Mechanics, 15(1): 113(in Chinese)
FINITE E L E M E N T A L G O R I T H M B A S E D ON T W O  S C A L E ANALYSIS METHOD 1
J.Z. CUI Institute of Computational Mathematical & Science/Engineering Computing, Academia Sinica, Beijing, 100080, ER.China
ABSTRACT The FE algorithm based on twoscale analysis method for the structural problems of composite materials and the structures with small period in 2dimension case is briefly presented in this paper, and some numerical results are shown. It is an effective method in computational mechanics to be developing.
KEYWORDS TwoScale Analysis Method, Composite Material, Structure with Small Period, Finite Element Algorithm.
INTRODUCTION The structures shown in Fig. 1 are often encountered in structural engineering and in the design of new industrial products. They are made from woven composite materials or composed of numbers of same basic configurations. The analysis problems for this kind of structure have some common properties as follows: The material parameters vary sharply and periodically, and the period ~ of material change is very
Project supported by National Natural Science Foundation of China
32
Cui, J.Z.
small; it means that, in mathematics and in mechanics,
a~hk(X+o~):a,jhk(X) ,
the coefficient of materials satisfy
1
=bjb>bj) are, respectively, the first and last blocks that couple block bj. They are determined by the profile of matrix A and the processor labeling.
PARALLELIZATION WITH ONESTEP OVERSHOOT We shall analyze the algorithm given above. For simplification, we assume that there is enough core memory to hold the matrix A across all the nodes of the computer. The computation can be classified into the sequential part BCTask(bj,bj) and the parallel part BCTask(bi,bj), bi>bj, where bi's are distributed across all nodes. The key to get higher speedup is to reduce the sequential part (including computing and waiting). We found that the bottleneck occurred in communication. At each loop of bj, all processors except the master processor iam [=mod(bjl,Np)+l] of block bj waited for the factored master block bj. In other words, all processors are synchronizedby the sending and receiving activities. In practice, the idle time due to communication and working load could be much longer than the computation time for most of the processors. This limits the extent of possible speedup improvement. We must reduce the idle tilne in order to improve the parallel performance. As an example, we consider an execution with 4 nodes, that are labeled from 1 to 4. Assume that the extraneous loop for bj=5, processor 1 is clearly the owner of block bj=5. While doing selfreduction BCTask(5,5), processors labeled 2,3 and 4 wait for receiving block 5 sent by processor 1. Once processor 2 receives block bj=5, it starts its own works BCTask(6,5), BCTask(10,5), BCTask(14,5) and so on. Assume the height of each column takes a constant value; processor 2 will have the heaviest working load in the parallel part of tasks. So the sequential part of step bj=6 will start ,,hen all other processors have finished their own tasks and have been waiting for the next master block 6. Clearly, the sending and receiving syuchronization of the master block is a bottleneck of parallel computation on distributed memory machines. For this reason, the parallel computation efficiency of LDL T in the JIK form on distributed memory machines is much lower than expected. The prerequisite to start BCTask (bj+l,bj+l) (A in Figure 2) depends only upon the complement of BCTask (bj+l, bj) (or in Figure 2) and not upon the complement of BCTask
92
Chen, P. et al.
(bj+l,bj+Np), BCTask (bj+l,bj+2Np).. (13,X in Figure 2). In conclusion, the blocked JIK form is not the only way to arrange the LDL r factorization.
m
xaiiiill .........
iiiiiTii iii £1111!
3
2 4 1
N~C~ 2
\
3
\4
Figure 2: Onestep Overshoot Concept (4 nodes) Using the asynchronous communication model for receiving, we adopt an algorithm with overshoot to minimize the idle time as shown in Fig. 1. In the design, the task order by the processor 2 follows the route AA'ot13X... rather than the route Aot13X...A' in the algorithm of Table 1. Clearly, the A' corresponding task BCTask (bj+l=6,bj+l=6) has been overshot into the parallel part of the step bj=5. Once the block bj+l =6 is factored, the master processor for block bj+l will immediately send the block to the others that may still work on the parallel part for step bj=5. This leads to the onestep overshoot scheme. Please note that in this algorithm synchronous communication model cannot be used.
NUMERICAL RESULTS AND REMARKS In order to evaluate the computational time and parallelism of the onestep overshot algorithm, we consider the flex joint of a tensionleg of an offshore platform. We discretize the flex joint using finite elements of different mesh sizes. Because of symmetry in geometry, material and loading, only half of the flex joint is considered in the analysis. Table 3 summarizes the parameters of the analysis, where mmmp denotes the number of nodes, neq denotes number of equations, and m is the average halfbandwidth.
Solving Large Systems of Equations on IntelParagon
93
TABLE 3 FINITE ELEMENT MODELS OF A FLEX JOINT model
nr
no
n~
numnp
1 2 3
8 1 10
20 35 35
23 24 32
4536 8100 13068
neq
m
12764 23155 37557
554 667 1071
In Table 4, the average wallclock time, Tp, required for onestep overshoot LDL T factorization with p processors running at HKUST is listed in Table 4. The corresponding results of the parallel speedup, Sp=TI/Tp, and the parallel efficiency, Ep=Tp/p, calculated on Tp are shown in the same table. The wallclock time reflects the parallel computation performance. In order to test the speedup and the efficiency equitably, all cases were executed four times for each given number of processors. In the table, BC and RS represenl broadcasting and ringsending models of communication[6,8,9]. All communications ar~ asynchronous. We have to mention that the recorded wallclock times are very discretized. Ir some cases of small number of processors, the standard deviation were as high as 30%, so th~ performance measurement was a difficult job. TABLE 4 PERFORMANCE OF ONESTEP OVERSHOOT ALGORITHM performance model sending term BC 1 RS
BC 2 RS
BC 3 RS
Tp Sp Ep Tp Sp Ep Tp Sp Ep Tp Sp Ep Tp Sp Ep Tp Sp Ep
number of processors 4 8 i 16
1
2
32
64
278.3 1.000 1.000
153.5 1.812 0.906 138 2.014 1.007 344.8 1.956 0.978 344.8 1.956 0.978
89.1 3.124 0.781 85.6 3.251 0.813 224.5 3.004 0.751 210.1 3.21 0.803
70.1 3.969 0.496 65.2 4.268 0.534 151.9 4.439 0.555 125.5 5.375 0.672
59.7 5.086 0.424 64.6 4.308 0.269 132 5.109 0.319 215.1 3.136 0.196
41.7 6.675 0.209 51.6 5.393 0.169 103 6.547 0.205 137 4.923 0.154
33.1 8.400 0.131 53.6 5.192 0.081 85.7 7.871 0.123 123.1 5.479 0.086
1188.1 1.937 0.968 1191.7 1.930 0.965
632.9 3.636 0.909 680.1 3.383 0.846
479.4 4.800 0.600 388.6 5.921 0.740
326.3 7.051 0.441 700.1 3.286 0.205
187.4 12.28 0.380 524.4 4.387 0.274
112.3 20.48 0.320 289.9 7.936 0.124
674.5 1.000 1.000
2300.8 1.000 1.000
94
Chen, P. et al.
As a reference, the wallclock time, can observe great difference.
Tp, required by the basic parallel algorithm is listed. We
TABLE 4 WALLCLOCK TIME FOR BASIC PARALLEL ALGORITHM wallclock time model sending BC 1 RS BC 2 RS BC 3 RS
1
2
278.3
171.5 161.0 370.1 369.8 1308.4 1295.1
674.5 2300.8
number of processors 4 8 16 115.6 265.5 231.1 681.9 681.1
506.3 497.4
467.7 450.4
32
433.1 358.4
64 78.5 65.6 199.0 144.4 389.2 308.8
Because of the size of problems, unavoidable I/O has influence on the performance, although the asynchronous read and write control are used in the coding. Research on quantitative influence of I/O has not been reported in the literature for MIMD machine. We obtained the maximum computation rate of 382.65 MFLOPS (multiplication and addition) for FE model 3 using 64 processors, which corresponds to 7.99% of the theoretical peak performance of 64x75 MFLOPS. This is 24.66% of the peak performance of dotproduct of Intel Paragon processor measured in the same machine. The best single node computation rate takes 18.7 MFLOPS for FE model 3, that is, 12.5% of the theoretical peak performance and 77.7% of the measured dotproduct performance. From this view, the influence of I/O cannot be over 22.3% of the wallclock time. We have the following observations: •
• •
Unlike the basic parallel algorithm reported in [6,8,9], broadcasting communication is in general better than the ringnode sending for the proposed algorithm. The exception is for the smallscale problem using only a small number of processors. The upper limit on the number of processors is 64, beyond which there is no additional speedup for the LDL T decomposition. The wallclock time is dominated by interprocessor communication when larger number of processors are used.
In summary, we have implemented the LDL v parallel factorization on a distributed memory machine such as the IntelParagon. Compared with conventional approaches, two improvements have been made to increase sequential and parallel calculation, respectively. Since interprocessor communication across network dominates the efficiency of parallel calculation, using an asynchronous communication algorithm with onestep overshoot can greatly reduce the processor idling time and enhance the performance of parallel computation.
Solving Large Systems of Equations on IntelParagon
95
ACKNOWLEDGMENTS The authors wish to acknowledge the support of this research effort by UPGC research funding. The first author wishes to thank Professor Duc T. Nguyen at Old Dominion University, Norfolk, Virginia, for his kind discussion.
REFERENCES [1 ] E.L. Wilson and H. Dovey. (1978). Solution or Reduction of Equilibrium Equations for Large Complex Structural Systems, Adv. Engrg. Software 1, 1925 [2] D. Zheng and T.Y.P Chang. (1995). Parallel Cholesky Method on MIMD with Shared Memory, Computers & Structures, 56:1, 2538 [3] J. Ortega (1978). Introduction to Parallel and Vector Solution of Linear Systems, Plenum Press, US A [4] O.C. Zienkiwicz and R.L. Taylor. (1990) The Finite Element Method, McGrawHill [5] P. Chen, T.Y.P. Chang and P. Tong, Implementation of a General Purpose FE code NFAP on Intel Paragon Platform, to be published [6] C. Farhat and F.X. Roux (1994). Implicit Parallel Processing in Structural Mechanics, Comp. Mech. Adv. 2, 1124 [7] C. Farhat and E.L. Wilson (1988) A Parallel Active Column Equation Solver, Computers & Structures, 28, 289304 [8] J. Qin and D.T. Nguyen. (1994). A ParallelVector Equation Solver for Distributed Memory Computers, Computing System in Engineering Journal, 5:1 [9] Qin and D.T. Nguyen. (1993). A New ParallelVector Finite Element Analysis Sottware on Distributed Memory Computers Proceedings of the AIAMASME/ASCE/AHS 34th Conference, La Jolla, CA [ 10]Intel Paragon TM User's Guide, Intel Corporation (1994) [ 11]Intel Paragon TM FORTRAN Compiler User's Guide, Intel Corporation (1994)
This Page Intentionally Left Blank
D Y N A M I C BEHAVIOR OF RAILWAY BRIDGES U N D E R RANDOM LOADING AND ASSESSMENT OF V E H I C L E  R U N N I N G SAFETY YingJun Chen, He Xia and Daqing Wang
Civil Engineering Department, Northern Jiaotong University, Beijing 100044, China
ABSTRACT
The dynamic model of a trainbridge system consists of the train model and the bridge model. Selfexcitations of the system are the track irregularity and the vehicle hunting movement. As the main external excitations, there are wind load and earthquake action. This paper studies this problem of different cases for the purpose of assessment of vehiclerunning safety and serviceabiliy; three methods are used and the theoretical characteristics of them are studied from the viewpoint of engineering application.
KEYWORDS
railway bridge, dynamic behavior, wind loading, earthquake, computer simulation, dynamic interaction, vehiclerunning safety, stochastic process, dynamic serviceability
A S S E S S M E N T BY E X P E R I M E N T AND S I M U L A T I O N
The
" s n a k e  h u n t i n g " movement by vehicle tracking is a principal selfexcitation
source. It can be expressed as follows
97
98
Chen, YingJun et al.
Y, =
2zrVt L, + ~';)
A,sin{
(1)
where A, and L, are the amplitude and wave length respectively, and ~ is the phase angle of ith wheelset; it is a random variable. A , = 3 . 0 m m ,
L,=20.0m. I
\
( LI
t/ x
~X~

64
=
~_%2m
_,
Figure 1.. Trainbridge system model In order to increase train speed, a truss bridge of 64m span is studied; the dynamic interaction model is shown in Fig. 1. Degreesoffreedom for movement are 23 and 17 for 6axle diesel locomotive and 4axle passengerfreight cars respectively. i
i
.
!
i
I
I
~
.
jw. ,
I
.
.
.
1
Z~.~
,
~~. ~ t~m ll.z.. J
,
it
.
 
,~,~;H t,!
~
It
'~t " it
!©©@ •
,,e
2SI,
' i
I
,...f
1
Figure 2: Model of diesel locomotive The vehicle model is shown in Fig. 2 E13. Horizontal torsional vibration is coupled with vertical vibration for it is a twotrack bridge, so a spatial model should be used. In order to simplify the calculation, the modal analysis method is used. The train on the bridge is directly concerned with the stringer. Suppose that there is no relative displacement between the stringers and the track.
When the vibration
Dynamic Behavior of Railway Bridges under Random Loading
99
modes of the bridge are normalized as ~rnrp,   1 , the nth order modal equation is expressed as follows
(2) i
j
in which F,,~,,F,,~ and Fo~j are horizontal, vertical and torsional force components when the j th wheel of the i th car acts on the bridge position, co and ~ are natural frequency and damping ratio respectively, ~;:,fL~and ~.~are horizontal, vertical and torsional displacement components at the stringer of nth mode at the position of ith car and jth wheel, respectively. When a train is moving on the bridge, lateral displacement Yw;j , torsional angle 0~;1, vertical displacement Z~;j of the wheelset, displacements of the girder Yb(3Cij),Ob (:C;j) and Zb(:r~j), hunting movements and irregularities of rails Y,(:c;j) ,O,(x~j) and Z,(scij) should satisfy the following relations
[Y,.~j~ [Z~J
[Yb(x~j) +
HiOb(Xij)
"Jr Y,(x;j)}
Zb(Z~j) + BOb(xij) + Z,(x u)
q,(~,j + Hi~i~) + Y,(xij) i Iq.:o,j + O,(x,j) I
l
,=1 q,(~ij + B ~ j ) + Z,(xij)
(3)
,
in which distances H and B are shown in Figs. 1 and 2. Dynamic spatial equilibrium equations can be obtained by combining the vehicle models and the bridge modes with the wheelrail relations (3).
MA 4:CA + K A   F
(4)
Matrices of mass M, damping C and stiffness K,generalized displacement A and generalized force F in Eqns. 4 and 2 can be written in detailed .form ~23. A similar problem is also studied in this paper for Fig. 10. The equations are solved by Newmark/9 method, ¢ / = 1/4. In this example (Fig. 1), the train consists of two locomotives and eight passengerfreight cars. In o r d e r t o investigate the effect of phase angle ~"in Eqn. 1 which is a random number of uniform distribution in 0 ~ 2n, the Monte Carlo simulation method is used; 30 sets are calculated for each train speed. Dynamic deflexions of the bridge are analyzed statistically.
Figs. 3 and 4 show the dynamic factor and horizontal amplitude(maximum, mean and standard deviation ~) at the midspan lower chord panel point respectively. But values higher than 200km/h are only reference values. They are also compared with the field
100
Chen, Y i n g  J u n et al.
I~5.
t~°
120 ~ (l+~t) i t.l, ~ Z
maximum _ mean 
1.05

/
"~'
/
6.0
~
t°°V .... ~o
maximum
9.0
7

amplitude(mm)
~
3.0 ~ .... ioo
i'
'~
~ ~
~
15o :oo :~o. Train speed (km/h)
:~
0.Q
30o
[ .... 1oo
~... r , ,,~.r 130 20O
.r .r F.r r ~ ~ ' 77.5O 30(1
Train speed (kin/h)
Figure 3 : Dynamic factor ( 1 ÷ # )
experiments.
....
Figure 4 .Horizontal amplitude of bridge lower chord
From the Chinese code, it is known that vehiclerunning safety is as
sured when maximum
train speed is 160km/h.
ASSESSMENT OF VEHICLERUNNING SAFETY UNDER EARTHQUAKE
BY USING THE WHEELRAIL INTERACTION MECHANISM Due to lack of data, the method by Wakui H. [3]is used. But in our case, maximum train speed is 300kin/h, degreesoffreedom of movement of the vehicle are 17 in order to study lateral responses, and earthquake acceleration record during the Tangshan earthquake is used. Vehiclerunning safety is controlled by horizontal bent angle 0y of the track~ track deformation yR can be expressed as follows(Fig. 5)
yR  0, O, (
YR =
X¢.
x
ZyX ,
(5)
Xc 2
L1
/ straight line ~ransitEon curve
.
X c
T
0
.
.
X
.
straight line
.
.
.
:
T
Figure 5: Track shape of bent angle portion
Figure 6: Tread surface slope of wheel.
When the train passes through this portion, a relative displacement between the wheel
Dynamic Behavior of Railway Bridges under Random Loading
101
creep force ~ K p flange n compressio
~
~tiiting
action
Figure 7: Cause of horizontal compressions[3] axle and the track may happen from the inertial force of the train; if this displacement is greater than clearance u, the flange may accept a shock ;the wheel radius also varies owing to the tread surface slope "/(Fig. 6); then the wheel axle may accept an action in the direction of the bent angle. This is the mechanism of the "snake hunting" movement. Fig. 7 shows the action o{ horizontal forces; creep forces also play an essential role in determining the lateral dynamic performance of ra{1 vehicles [4~. Fig. 8 shows a mechanical model of vehicle; it consists o{ one car body, two bogies and •
four wheelsets connected with various springs; all data of the Japanese Shinkansen car [3] are used in this study. The total degreesoffreedom o{ movement are 17 for one car, for car body: lateral moving (yB) , rolling (asB) , yawing ( grB ) ; f o r bogies (i = 1 , 2 ) : lateral moving (YT;), rolling (aST;) , yawing (grT;) ; for wheelsets (i = 1 , 2 , j 1,2) :lateral moving (y~;j), yawing ( g r )
, w h e r e / = 1 , 2 r e p r e s e n t s front and rear
bogie respectively, j = 1,2 represents former and latter axle respectively ; the effects of rolling as~,~ and vertical motion zwij can be neglected; it is different from [3] for earthquake acceleration ag is considered in this case. The springs are divided into two kinds (Fig. 8 ) . one is the nonlinear spring with damping; it consists of vertical support spring K2 , horizontal spring K3 , axle box spring K1 ; the other is linear spring without damping. For the bolster anchor spring K0 , suppose the relation between the relative yawing angle and the moment has bilinear hysteresis characteristic. The resultant F from spring stiffness K and damping C can be obtained from spring displacement q and its velocity; q has a limited value q, by a stopper, for example q~;~ = ~ y~ ± ( 1);L~B ::Ph2¢B 4 YT~ =~ h,CPT~
102
Chen, YingJun et al.
b,
K,. C, \ j }
X
Z.T~
X l
b,
K,.
.i
Side A
SKIa B
mwl side
A
sKJe B
F i g u r e 8: Mechanical model of v e h i c l e [ 3 ] [qai~ I > qa,,
Fa,~
q~i~ [qa,; I { K , q 3 , + K s , ( Iqai~ ]  q3,) } + C,c).;~
(6)
T h e o t h e r n o t a t i o n s are 2m~: mass of car b o d y , 2IB=and 2IB, : m o m e n t of inertia of the car a b o u t x axis and t h a t a b o u t z axis, m r : mass o{ bogie f r a m e , IT= and Iv, : m o m e n t of inertia of bogie a b o u t x axis and t h a t about z axis, m w : mass of w h e e l a x l e , 1~ : mom e n t o{ inertia of w h e e l a x l e a b o u t
its center o{ g r a v i t y , 2L: distance of t w o b o g i e s ,
c e n t e r to c e n t e r . T h u s the d y n a m i c equilibrium equations of the car are as follows: 2mByB = F3.1a  F3.1~ + F3.ZA  F3.zB + 2tuBa,
I
car body 2 I n z # B = (F3.1A  F3.1B  Fa.ZA "Jr F3.zB)L 
(M1 + M2)
(7)
2I~@B = (Fz ,A  Fz ,B + Fz zA  Fz zs)bz + (F3 ,A  F3.,B + F3,ZA   F3.z~)h2
rnTYTi =  
] bogie~ [
[
l,
F3.,A + F3.,B + F w y . i , t "Jc" Fwy.izA   Fwy.iiB   Fwy,izB + m T a g
ITz'~Ti = Mi + (Fw..i,A  F v ¢ = , B  Fw,..IzA + Fw..izB)bl +
(8)
(Fwy.i,,1  F w y . i , ~  Fwy.i2A + Fwy.i~s)a I T ~ T i  (   Fz iA + Fz in)bz + (F3 iA  F3.is)h, + (FI.qA  Fl.i,s + Fi.iza   F l , i ~ ) b 1 + ( F w y . i , A   Fwy.; B "Jr Fwy.iza   Fwy.izS) ( H T
 r)
Dynamic Behavior of Railway Bridges under Random Loading
w h e e l s e t j m w y . . j =  Fwy.ijA { Fwy.us  QijA Jr" QijB Jr mwa. I~~,,.j = (   1)~(Fw..;jA F w . . u s ) b l  T~ u
103
(9)
in which horizontal compression Qu; = QJ~ + Qs;;; ,Q,,j~ .... is horizontal creep force , Qfu,~ is the flange compression (Fig. 7), Qf,j~  K p • Y,o~ , where K p is the spring constant f r o m tilting action of the rail, y,;j; is the railtilting displacement when the flange collides against the rail. Restoring moment M of the bolster anchor spring is M =
2KobZoCE
, relative yawing angle aF = aYE  g"r , a n d aFE is a value of ~ within the elastic limit, when g" > We. . . . M has a bilinear relation [3]. T~,,.i is yawing m o m e n t from the creep action. Equilibrium equations (7) ~ (9) are a timedependent s y s t e m ; on the basis of the above idea, simulation analysis is used. The initial condition is that the train passes in uniform speed at the position of the beginning point of the bent angle portion when t  0 , and the earthquake acts at the same time.
lii 1.6
.,_.1 (.
(1)
.  7...'.:"; '''''''; r'''''" r''':'''''"
B1
(1) o 0.8 i
.':¢::/:
, ::.::':"
(.)
:::iiii lt
~0.6 u 0.4
::;.:'
A 1 .... .:::?t.:::""

0.2
00 F
o.ooo
J
I
0.002
~
I
,,
a
I
,
I
i
o.29
I
0.004 0.006 0.008 0.010 (v=3OOkm/h, ag=O. 1 gO. 3g)
J
bent angle
Figure 9: Horizontal bent angle versus derail coefficient
Fig. 9 shows the relation b e t w e e n the horizontal bent angle (in radian) and the derail coefficient when train speed v = 3 0 0 k m / h and earthquake acceleration a, = 0 . l g " ~ 0 . 39; in the figure, A indicates inner rail, B indicates outer rail. In the same m a n n e r , this angle can be calculated for various parameters (derail coefficient, side compression of wheel flange, horizontal vibration acceleration of c a r ) a n d for other vehicle models ( T G V , ICE) and other design earthquake actions. Then the allowable bent angle value can be decided in order to.satisfy vehiclerunning safety. In our case, the calculated
104
Chen, YingJun et al.
results are shown in Table 1; earthquake acceleration is a, = 0 . 2g, it approximately corresponding to Richter Magnitude M  6 , Modified Mercalli Intensity I = 7 ~ 8.
TABLE
1
A L L O W A B L E BENT ANGLE
Horizontal Bent Angle 0 (1/1000)
Train
Horizontal Moving (1)
Speed km/h
L x is
P(X > x , V < v) = f
f (x,v)dvdx
:r 
crz
x
oo
=f where
f(x,v)
f ( x tv)
f(xlv)f(v)dvdx= ~
[1  F ( x l v ) ] f ( v ) d v
(15)
0
is the conditional probability density function of the extreme responses,
is the joint probability density function.
Structural serviceability is defined as the probability of which structural response extreme X does not surpass its.threshold within service period T .
This probability is an
evaluation index for the bridge'train system reliability and can be called dynamic serviceability. Since
FxT(X
%
X It " ( T )
=
] 
PX,T(X > x,t < T) ,
we have
v(T)
FxT(X < x It < T) = ]  f
exp[
e",~"~o>]eO,("~o)a,,
0
×
{ 1   e x p [   e x p {   a o + a l v 41_ a 2 v 2ix__ (bo_4_blv+b,vZ)])]}dv
(16)
where the upper limit of the integration is the maximum expected mean wind speed de
108
Chen, YingJun et
al.
termined by the r e t u r n period. Taking a bridge in south China for e x a m p l e , the local annual e x t r e m e parameters for wind speed a r e / 1 , ,   1 4 . 9 ~ , m / s and a,, = 0. 417.
For a
given return period T  100year, the maximum mean wind speed is 2 6 . 0 m / s . T h e average passage time
T,(x) is the inverse of the probability of which an event hap
pens in one unit period 1 T.(x)
1
 Px,T
(17)
1  FX,T(X)
The threshold for response X can be calculated from a given service period To .
TABLE
2
D Y N A M I C S E R V I C E A B I L I T Y I N D E X E S AND T H E I R REFERENCE THRESHOI.I)S
Horizontal
Bridge
Vehicle
Accleration (g)
Acceleration(g)
24.8
0.2
0.50
Amplitude (mm)
Deflexion ( m m ) 56
28
Japan
UIC
China
Japan
Japan
7" = ] o o
0. 9291
0. 4262
0. 2657
O. 9024
O. 9997
jI
O. 9999996
V =30.0
O. 9208
0.4178
0. 2574
0. 8947
O. 9993
tI
o. 999993
7.8 ¸
Threshold 1
China
The calculated dynamic serviceability indexes ( D S I ) are shown in Table 2 for 7" =
l OOyear, V = 26 m / s and V = 3 0 m / s ; the latter is the m a x i m u m wind speed for an empty car. For such a complicated system as bridgetrain interaction, since wind and train loads do not correlate with each o t h e r , all the DSIs are neither fully correlative nor independent. The system DSI can therefore be estin~ated by the following equation
l[ Fxi ~ F, ~ min{Fx, }
(18)
In this e x a m p l e , for T = 1 0 0 years and the threshold of Japanese code, the system DSI is b e t w e e n 0.8381 and 0. 9024. Figure 11 shows the distribution of the calculated DSI versus m a x i m u m mean wind
Dynamic Behavior of Railway Bridges under Random Loading
109
speed. T h r e s h o l d s for lateral deflexion, amplitude and acceleration of this bridge for service p e r i o d T =
1 0 0 y , inversely calculated by Eqn. 17, can be 7 0 m m , 45mm and
160gal respectively. 1.0
. . . . . .
. .
,
0.2g(acceleration) 0 1g(acceleration) 56mm(deflexion) 24.8(amplitude)
. . . . .
"... ', •. ',
0.8
";
0.6
• %
•.
%
• .
0.4
.
28mm(deflexion)
'".,... °.
0.2 0.0
0 !
~
5 I
,
1'o
,
15 I
,
..................
2XO
,
. .....
2z5
,
7.8mm(amplitude)
3XO
,
v(m/s)
Figure 1 1 : System DSI versus wind t h r e s h o l d
For this kind of b r i d g e , m e t h o d s for design and a s s e s s m e n t aiming at a s s u r a n c e of bridge serviceability for train and wind interaction should be emphasized. In this example, the e s t i m a t e d DSIs for lateral bridge deflexion and amplitude are only 0. 9 or so upon the local wind speed distribution and design period 7 " =
lOOyear, not suitable to
the Chinese code. CONCLUSION A l t h o u g h the c o m p u t e r simulation can be used for the research of t r a i n  b r i d g e s y s t e m , some difficulties remain to be further studied, for i n s t a n c e , the power s p e c t r u m of rail i r r e g u l a r i t y , the vehicle hunting m o v e m e n t function and its a m p l i t u d e , wave length and phase a n g l e , unless it is a case study. It is more rt:~tsonable to use the m e c h a n i s m of wheelrail i n t e r a c t i o n , but the exact data of p a r a m e t e r s are necessary. For the ass e s s m e n t of s e r v i c e a b i l i t y , vibration thresholds to ensure t r a i n  r u n n i n g stability and h u m a n c o m f o r t , d e t e r i o r a t i o n of s t r e n g t h caused by external loading and e n v i r o n m e n tal c o n d i t i o n s , etc. are now further studied.
ACKNOWLEDGEMENT This research is s u p p o r t e d by National N a t u r a l Science F o u n d a t i o n and Ministry of Railways of China.
Chen, YingJun et al.
110
REFERENCES
1. Chen YJ. , Wang D. (1992),Dynamic Behavior of Bridges Under Random l.oading and Dynamic Reliability Problem, Computational MecJ;anics in Structural engineer
ing, edited by Cheng F. Y . , Fu Z., Elsevier Science Publishers, England, 398"~ 411 2. Xia H . , Chen Y  J . , Zhang D . , Ke Z. (1996), Dynamic Analysis of Steel Truss Bridges Under Increased Train Speeds (in Chinese), Journal o f the China Railway
Society 18 : 5, 75~82 3. Wakui H. (1978), Allowable BentAngle of LongSpanned Suspension Bridges Determined by Running Property of Shinkansen Car (in Japanese), Report o f Railway
Technical Research Institute, No. 1087, 1 ~ 5 9 4. Shen Z. Y . , Hedrick J. K . , Elkins J. A. (1983), A Comparison of Alternative Creep Force
Models for Rail Vehicle Dynamic Analysis, Proc.
8th I A S V E
Symposium, Swets ~ ZeitlingerLisse, 591"~605 5. Chen YJ. (1992), Development in Dynamic Reliability of Bridge Structures (in Chinese), Proc. 6th Symposium on Bridge Technology,Wuhan, China, 76"~81 6. Chen YJ. et al. ( 1 9 9 3 ) , Development of Wind Load Standard for Chinese Bridge Code, Proc. Third AsiaPacific Symposium on Wind Engineering, Hong Kong, 793"~798 7. Xia H. , Chen YJ. (1994), Dynamic Reliability of TrainBridge System Under Wind Action (in Chinese), Journal of China Civil Engineering Association 27 : 2,14"21 8. Chen YJ. , Xia H. , Gao R. (1991), Dynamic Analysis of Truss Bridge Stiffened by Flexible Arch Unter Train Loads, Proc. E A S E C 
3, Shanghai, 1451"1456
A REVIEW ON THE NUMERICAL SOLUTION SCHEMES FOR LOCALIZATION PROBLEMS
Z. Chen, X. Xin and D. Qian Department of Civil and Environmental Engineering, University of Missouri Columbia, Missouri 65211, USA
ABSTRACT An introduction is given on the current status of research on localization problems. Numerical solution schemes to simulate the evolution of localization are then reviewed in detail. First, material integration schemes are discussed for both plasticity and damage models for which a continuum tangent stiffness tensor can be formulated for bifurcation analysis. Second, structural integration schemes are reviewed for the postlimit responses including snapback and snapthrough. Finally, one of the "meshless" methods, the material point method (MPM) is introduced for localization problems. Sample problems are then considered to demonstrate the robustness and potential of a newly proposed procedure based on the moving jump forms of conservation laws. Conclusions and future work are given based on the review paper.
KEYWORDS Softening with Localization, Plasticity and Damage, Postlimit Response, Jump Conditions, Transition between Governing Equations, Failure Wave
INTRODUCTION For modem engineering design, the limit design methodology based on elastic or yield strength might not be suitable in many cases. In dynamic cases such as impact or seismic disturbances,
lll
112
Chen, Z. et al.
the loads are of very short duration so that collapse may not occur even if the limit point is reached. Engineering structures are often designed to be statically indeterminate which also helps to preclude collapse. In addition, the current interest in explosionresistant design requires the understanding of the postlimit structural responses. Since a significant part of energy dissipation in the postlimit regime is associated with the evolution of localization, much research has been conducted to investigate the experimental, analytical and numerical aspects of localization problems, as reviewed by Chen [1996b], Chen and Schreyer [1994], and Xie et al. [ 1994]. Localization is manifested by softening or degradation of material properties, which is accompanied by localized large deformations in a finite zone. To predict the evolution of inhomogeneous interactions among material particles within the localization zone, several kinds of unconventional constitutive models have been proposed with the use of higher order terms in space and/or time, such as nonlocal (gradient or integral), Cosserat continuum and ratedependent approaches. Although the advantages of these unconventional models over conventional (local) ones have been demonstrated for many academic problems, there still exist some pressing limitations that prohibit the routine prediction of localization phenomena in practical applications. In particular, it might not be feasible, with current computational facilities, to perform largescale simulation of structural failure due to the use of higher order terms. Since the shift from a testbased to a simulationbased design environment requires an efficient numerical procedure for localization problems, an altemative approach is to catch the essential feature of localization phenomena without invoking higher order terms for constitutive modeling. As can be found from the literature review [Chen and Schreyer, 1994; Chen and Sulsky, 1995], the key component of various higher order models is an attempt to predict the evolution of inhomogeneous interactions among material particles, with meshindependent results. Before the initiation of localization, there is no need to use higher order terms. In a macromechanical sense, however, the evolution process might be equally well characterized by the formation and propagation of a moving material boundary that is associated with a local change in material properties. With the introduction of a moving material boundary, a partitionedmodeling approach has been proposed for localization problems [Chen, 1993a]. The basic idea of the approach is that different local constitutive models are used inside and outside the localization zone, with a moving boundary being introduced between different material domains, if localization occurs. As a result, the extrapolation of material properties beyond the limitations of current experimental techniques can be avoided in identifying the evolution of localization, and simplified governing differential equations can be formulated in the partitioned domains for given boundary and initial conditions. To establish a sound mathematical foundation for the partitionedmodeling approach, an attempt has been made to investigate the use of the jump forms of conservation laws in defining the moving material boundary, with onedimensional analytical illustrations for rateindependent local models [Chen and Sulsky, 1995]. By taking the initial point of localization
Review on the Numerical Solution Schemes for Localization Problems
113
as that point where the type of the governing differential equations changes, i.e., a hyperbolic to an elliptic type for dynamic problems and an elliptic to another elliptic type for static problems, a moving material boundary can be defined through the jump forms of conservation laws across the boundary. Jumps in density, velocity, strain and stress can be accommodated on this moving surface of discontinuity between two material domains. Interestingly, the problems involving the change in the types of governing differential equations also occur in other areas such as fluid mechanics [Chen and Clark, 1995] and thermal shock wave propagation [Tzou, 1989]. It has been shown that the transition from continuum damage mechanics to fracture mechanics might be linked through the moving jump forms of conservation laws so that a complete failure evolution process might be simulated with the use of simple models [Chen, 1996a and b]. Based on the previous research results, numerical solution schemes for localization problems are reviewed in this paper. First, material integration schemes are discussed for both plasticity and damage models for which a continuum tangent stiffness tensor can be formulated for bifurcation analysis. Second, structural integration schemes are reviewed for the postlimit responses including snapback and snapthrough. Finally, one of the "meshless" methods, the material point method (MPM) [Sulsky et al., 1994; Zhou et al., 1998] is discussed for localization problems. To demonstrate the robustness and potential of a newly proposed procedure that is based on the moving jump forms of conservation laws, sample problems are considered for both quasistatic and dynamic cases. Conclusions and future work are then given based on the review paper. A direct notation is employed to describe the constitutive models, with boldfaced letters denoting tensors of first or higher orders.
EXISTING NUMERICAL PROCEDURES Structural solution schemes consist of constitutive model solvers, spatial and temporal discretization methods. With a focus on the localization problems, the numerical solution schemes are reviewed as follows for nonlinear structural analyses including failure simulation.
Material Integration Schemes To perform largescale computer simulation, a simple stressstrain relation, that can predict the essential feature of nonlinear material responses, must be formulated for structural analyses. Nonlinear material behaviors arise from two distinct modes of microstructural changes: one is plastic flow and the other is the degradation of material properties. Plastic flow, which is reflected through permanent deformation, is the consequence of a dislocation process along preferred slip planes as in metals, or particle motion and rearrangement as in geologic materials. Because the number of bonds between material points is hardly altered during the flow process, the material stiffness remains insensitive to this mode of microstructural motion, and change of strength is reflected through plastic strain hardening
114
Chen, Z. et al.
and apparent softening. On the other hand, the nucleation, crushing and coalescence of microcracks and microvoids result in debonding, which is reflected through the damaging of material stiffness and strength. In general, both modes are present and interacting although some mode might dominate at some stage of the evolution process. Final rupture occurs when macrocracks form and propagate from the cluster of microcracks. A systematic procedure, which satisfies thermodynamic restrictions, has been used to formulate local plasticity and damage models [Chen and Schreyer, 1990b and 1994]. A family of incrementaliterative integration rules has also been given with or without the use of a tangent stiffness tensor. In general, both the generalized trapezoidal and midpoint rules can be employed for the integration of inelastic constitutive equations, depending on the feature of a specific model.
Structural Integration Schemes There exist two major computational difficulties in simulating the postcritical structural response. One is the occurrence of an illconditioned tangent stiffness matrix around critical points, namely, limit and bifurcation points. The other is the selection of a suitable constraint on the solution path such that the postcritical response can be traced. Since a robust and efficient solution scheme is necessary to make failure simulation available in a routine manner, several procedures have been proposed to circumvent the difficulties associated with critical points [de Borst et al., 1993; Chen, 1993b and 1996b; Chen and Schreyer, 1990a and 1991; PijaudierCabot and Bode, 1995]. The standard arclength control is still commonly employed in geometrically nonlinear cases, and with some modifications in materially nonlinear cases. Because the material failure zone is localized into a small region, the arclength constraint formulated in the global deformation field is insensitive to the evolution of the localized deformation mode. As a result, a suitable constraint should be constructed in terms of a localized kinematical field if localized failure needs to be simulated. Due to the fact that the localization zone is evolving and there is a sign change of the load increment at a critical point, the control point (element) should also vary in position with a suitable measure of failure, and the localized constraint parameter that reflects the extent of irreversible energy dissipation should be constrained to increase monotonically. Preliminary results obtained for plasticity and damage problems indicate that the use of an evolving localized control is a reasonable choice for localization problems including snapback or snapthrough [Chen, 1993b; Chen and Schreyer, 1990a and 1991 ]. As a remedy to avoid the use of an illconditioned tangent stiffness matrix around critical points, a secant stiffness matrix based on continuum damage mechanics has been used with a dramatic increase in the rate of convergence with respect to the rate obtained using a tangent stiffness matrix. In order to incorporate both damage and plasticity models into one computer code, it has been proposed that an incrementaliterative solution scheme be constructed through the use of an initial elasticity stiffness matrix together with an evolvinglocalization constraint [Chen and Schreyer, 1994]. Thus, only one inverse calculation is required, and the amount of
Review on the Numerical Solution Schemes for Localization Problems
115
computation involved in the iterative loop is dramatically reduced. It can be found that the numerical procedure for tracing the postbifurcation path depends mainly on the choice of a suitable constraint instead of a stiffness matrix. Different constraints generally yield different solution paths. In other words, the constraint imposed on the solution path plays a crucial role whether or not the stiffness matrix is illconditioned. Thus, the dependence of the evolvinglocalization constraint on the location of initial imperfections and on the evolution history of localization detects the solution path following the critical point, while the use of a wellconditioned stiffness matrix guarantees that a numerical solution can be obtained. The Material Point Method
As one of the innovative spatial discretization methods, the Material Point Method is an extension to solid mechanics problems of a hydrodynamics code called FLIP which, in turn, evolved from the ParticleinCell Method. The motivation of the development was to simulate those problems, such as penetration, perforation, metal forming and cutting, which involve large deformations, the transition from continuous to discontinuous failure modes and the creation of new material surfaces, with historydependent internal state variables. The essential idea is to take advantages of both Eulefian and Lagrangian methods. And also, this MPM and other unconventional spatial discretization methods employ the concept of local moving interpolation so that local remeshing can be achieved without the cost of global remeshing. Although the MPM is still under development, sample calculations have demonstrated the robustness and potential of this method [Sulsky et al., 1994]. It is believed that the MPM can be developed into a robust spatial discretization method, combined with solid modeling and postprocessors, for largescale computer simulation of structural failure responses.
A DEVELOPING NUMERICAL PROCEDURE Recently, efforts have been made to simulate the evolution of localization without invoking higher order terms, and to fill the gap between continuum damage mechanics and fracture mechanics. A developing numerical procedure is presented as follows. Based on the previous study [Chen 1996a], a set of moving jump forms of the conservation laws is used here to define a material failure criterion that can predict the initiation and orientation of localized failure, and the transition between continuous and discontinuous failure modes. In a threedimensional framework, the use of the jump forms of conservation of mass and linear momentum would result in Vb'n PlVl.nP2V2 .n
PlP2 and
(1)
116
Chen, Z. et al.
'01'02n'(v1v2)l(vlv2)
(0.1  0"2)" n
(2)
,01  192
if there is a jump in mass density p. In Eqs. (1) and (2), the subscripts 1 and 2 denote field variables on the two sides of a moving material surface that has a velocity vb in the threedimensional space, and n is the unit normal to the material surface. The other variables are v (particle velocity vector) and 0. (stress tensor), with body forces being omitted. For purely mechanical problems where there exists no energy sources or sinks, the conservation of mass and linear momentum implies the conservation of energy. Since the evolution of localization involves jumps in certain field variables, a material failure criterion is defined based on the jump types of the kinematic field variables as follows" Localized Failure :=
v l = v2 a n d i~1 :/: i~2
Discrete Failure :=
vI # v2
and ~1 ~ ~2
with ~ being the strain rate tensor. In other words, the change in the jump type of the kinematic field variables identifies the initiation of different failure modes. As can be seen, the transition from continuous to discontinuous failure modes is characterized by the condition of localized failure. Because there is a jump in the strain rate for localized failure, it makes sense to claim that a corresponding jump must exist in the mass density due to a jump in the volumetric strain rate. The jump in the volumetric strain rate is manifested by microcracking in the localization zone. The use of Eqs. (1) and (2) then yields vl.n
= ve.n
= v b. n
(3)
and the continuity of the traction across the moving material surface, namely
((71
(~2)" n = 0
(4)
for localized failure. Thus, Eqs. (3) and (4) together with a jump in the mass density represent the essential feature of localization. To examine how the jump in the strain rate is derivable from Eq. (4), assume that side 1 initiates localization from a weak material point. The stress tensor on side 1 is then related to that on side 2 by cr1 = cr~ + T l: Aek
(5)
in which T1 denotes a fourthorder tangent stiffness tensor with minor symmetries, and A e k is the incremental strain within side 1 due to the evolution of localization. According to Maxwell's compatibility conditions, A e k must be a rankone tensor of the form
Review on the Numerical Solution Schemes for Localization Problems
1
Aek = d7. ( A m O n + n ® A m ) 2l
117
(6)
with l being the halfband width of the localization zone. The use of Eqs. (4)  (6) then yields the classical necessary condition for a discontinuous bifurcation or loss of ellipticity" Q • Am = 0
(7)
with Q = n . T l . n being the acoustic tensor. An eigen analysis can be performed to find out the orientation of localization, A m , corresponding to a zero eigenvalue of the acoustic tensor. However, the magnitude of A m depends on the evolution of localization. If A m is determined from a constitutive model, then the jump in the strain rate can be found. To determine the jump in mass density, an evolution equation for damage must be defined. It is assumed that the shearinduced cracking in a representative volume is governed by a strainbased damage surface [Chen and Xin, 1997], f
E  ~o(1 + moD )
(8)
where ~' is the second invariant of deviatoric strain tensor, ~o is the critical state parameter, D is damage, and mo is a model parameter. With the use of a standard procedure (Chen and Schreyer 1994), it can be shown that the damage surface satisfies the thermodynamic restrictions and the rate of damage is determined by e
b  ~._
(9)
eomo
To represent the overall effect of shearinduced local dilatation on the change in mass density, an integral average of damage over the representative volume, D, is used to find 01, namely,
P' P2 +(ProP2)(1em'~) with
(10)
m 1 being a model parameter. As can be found from Eq. (10), 191 =P2 if there is no
damage, and P~ will approach the maximum value Pm with an increase in D. If Pm is reached, phase transformation might occur, a further discussion on which is beyond the scope of this paper. Because of P~ > P2 with the evolution of damage, the localization zone will expand based on Eq. (1). A partitionedmodeling approach is employed in the proposed numerical procedure with a totalstraincontrol scheme, as discussed next. Since the condition of v 1 • n = v 2 • n = v b" n
118
Chen, Z. et al.
holds across the moving material surface between the damaging and unloading zone, the increase in the total strain inside the localization zone would expand the zone. A typical example is the formation and propagation of a shear band under uniaxial compression as shown in Fig. 1. As can be seen, local microcracking (dilation) would push the materials inside the band to the boundary of the band so that pl > 192 across the moving boundary.
I()
()
()
~
()I
.°,..e. "'r
o.,e.,o. °°'
Figure 1: Evolution of shear band after failure occurs To implement a moving material surface into the numerical procedure, a local remeshing process is proposed here [Chen et al., 1997]. If the eigen analysis on the acoustic tensor indicates that the material is still in the prelimit stage, no change needs to be made to the original mesh. If material failure occurs, the moving speed of the material surface between the elastic unloading zone and the localization zone is determined based on the jump conditions. For the orientation of the localization zone given by the eigen analysis, a simple remeshing process can be demonstrated through a 2D example as shown in Fig. 2. If the weak element is located in the center, a shear band will occur between the boldfaced lines. The unit normal obtained from the eigen analysis is given as n = (costx, sintx) in the 2D space. Due to the evolution of the failure process, a new material surface is shown by the dashed lines. Assuming the coordinate for a node along the top surface of the shear band is (xl, Yl), and the magnitude of the moving speed for the current incremental step is vb,, then the new coordinate of this node (Xl ,Yl ) is given by X 1  X 1 + V b i A t " COS a
y( = y~
+ VbiAt .
sin a
(11
a)
(11 b)
in which At is the time interval. The same operations will be performed for the other nodes along this moving material surface. It should be noted that when a shear band forms, there are two moving material surfaces, both departing from the original shear band in the contrary direction. Hence, an upper and a lower surface with corresponding normals should be defined based on the location of the node.
Review on the Numerical Solution Schemes for Localization Problems
Ly
0
XX XX XX XX XX
119
Y New Material Surface n
Old Material Surface
!
/
Old Material Surface
New Material Surface
I//// / ~. X
O
Y ~.X
Figure 2 : A 2D demonstration for the moving material surface and the remesh process at a specific node In summary, the proposed solution scheme for simulating the evolution of localization consist of the following steps: Step 1: At the beginning of a new incremental step, an eigen analysis is performed on the acoustic tensor for each element. If failure is detected, find the unit normal to the moving material surface and record the element number, then go to next step; If no failure is indicated, go to step 4. Step 2: Based on the jump conditions, the moving speed of the material surface can be obtained. Step 3" Starting a loop for the elements in which material failure occurs, and update the configuration for the nodes inside the localization zone by using Eq. (11). Attention should be paid to those nodes shared by two or more elements: only one update in the configuration can be made for each of these nodes in a single shear band case. For the nodes inside the elastic regime, no change needs to be made. Step4: Return to the main program which employs the totalstraincontrol algorithm. In the proposed procedure, several assumptions have been made as a result of numerical approximation. First of all, in order to get the displacement of the node at the new location, certain interpolation is performed in the displacement field, depending on the assumption used for a specific type of elements. For example, linear interpolation is applied in the case of classical triangular element. In addition, since there is a change in the nodal position along the material surface, we can expect a change in the shapes of the corresponding elements both inside and outside the localization zone. In this case, we assume that the basic variables such as stress, strain and internal variables are the same as the values from the last step for these elements. For a small incremental step, this is reasonable and therefore the algorithm can be
120
Chen, Z. et al.
greatly simplified. The remeshing process is designed so that the size of the localization zone is updated based on the moving jump conditions.
DEMONSTRATIONS To demonstrate the robustness and potential of the proposed procedure that is based on the moving jump forms of conservation laws, sample results s are presented for both quasistatic and dynamic cases.
A Ly
~X
I~
Lx
Figure 3" The geometry and notation for a plane problem For the quasistatic loading case, the problem geometry is given in Fig. 3. The sizes of the specimen are L x  3m and L y  6m. Rectangular elements composed of four triangles are employed for the simulation. The load is applied such that all the points at the top surface will have the same vertical displacements. In addition, all the points at the bottom are constrained in the vertical direction and no horizontal displacement is allowed for the left most point. An associated von Mises model with bilinear hardening and softening function is used, with Young's Modulus E  50GPa and Poisson's ratio v  0.2. An initial imperfection is introduced by assigning a weaker limit strength in the weak element. With the use of a free mesh, the deformation pattern with the local von Mises model is shown in Fig. 4, and the corresponding deformation pattern with the proposed procedure is given in Fig. 5. As can be seen, the use of a moving material surface expands the localization zone. Figure 6 demonstrates the meshindependence of numerical solutions with the proposed procedure.
Review on the Numerical Solution Schemes for Localization Problems
121
V(m)
0
X(m) 0
1
2
3
Figure 4: Deformation pattern with local von Mises model for mesh III
0
X(m) 0
1
2
3
Figure 5" Deformation pattern with proposed procedure for mesh III k LoadP (N) 2.0e+08
Mesh I Mesh II
1.5e+08.
Mesh III %
1.0e+08.
5.0e+07
O.Oe+O0 0.00
'
I
0.02
'
I
0.04
'
I
0.06
'
I
0.08
'
I
v
Displacement A (m)
0.10
Figure 6" Meshindependent solutions with the proposed procedure For the dynamic case, a bar of length L  l m is considered here [Chen et al., 1997]. The bar is fixed at x=O and loaded at x = L . A Newmark integrator is used in the time domain with the
122
Chen, Z. et al.
time step satisfying the standard stability condition. After the limit state is reached at the fixed end, the evolution of localization is shown in Fig. 7, and the corresponding decrease of stress is given in Fig. 8. Figure 9 demonstrates the convergence of numerical solutions with different meshes at t = 3 0 0 p s .
& 0.010 ~r E
/
0.008 
t = 300 Its
0.006 
0.004 0.002
t = 250 l~s
~ u
0.000
I
0.2
0.0
t = 275 Its
I
0.4
I
0.6
I
I
0.8
1.0
x (m)
Figure 7" The evolution of localization after the limit state is reached O" ( M P a )
100
t = 300 ~s
80
t = 275 l~S
60 40
......
=
20 0
X (m) 0.0
0.2
0.4
0.6
0.8
1.0
Figure 8: The decrease of stress corresponding to Fig.7
0.010 t
E,,
0.008
n e = 10
0.006 
n e = 20
x
,,
0.004 0.002  ~ 0.000 0.0
n e = 40 ~ ,
i
0.2
0.4
I
0.6
I
0.8
I " x (m)
1.0
Figure 9: Convergence study with different meshes at t = 300 ~ts
Review on the Numerical Solution Schemesfor Localization Problems
123
CONCLUDING REMARKS Existing numerical procedures for simulating the evolution of localization have been reviewed in this paper. To simulate the evolution of localization without invoking higher order terms, a developing numerical scheme has also been introduced. As a result, a simple solution procedure with the use of local models can be designed for localization problems, if localization is considered as a phase transition (diffusion) process. In this procedure, the initiation of localization is identified via monitoring the transition between different types of governing differential equations, and the evolution of localization is traced by using a moving material surface of discontinuity. The proposed p~ocedure has been demonstrated through sample problems under dynamic (an elliptic equation inside the localization zone and a hyperbolic one outside the zone) and quasistatic (two different elliptic equations inside and outside the zone, respectively) loading conditions. The numerical results are consistent with the essential feature of localization phenomena. Further study is required to apply the proposed procedure to a general case, which involve multidimensional effects, large deformations, discontinuity, anisotropy and phase transition. Although several promising approaches have been proposed for spatial discretization in localization problems, it is still a challenging task to tackle different types of goveming equations in a single computational domain, because of different temporal and spatial scales. Since the design of advanced engineering structures needs the solutions for this kind of problems, it is hoped that more innovative ideas will come from the international research community in the near future.
ACKNOWLEDGMENTS This work was sponsored (in part) by NSF and AFOSR.
REFERENCES
R. de Borst, L.J. Sluys, H.B. Muhlhaus and J. Pamin (1993). Fundamental Issues in Finite Element Analyses of Localization of Deformation. Engineering Computation 10, 99121. Z. Chen (1993a). A PartitionedSolution Method with Moving Boundaries for Nonlocal Plasticity. Modern Approaches to Plasticity, edited by D. Kolymbas, Elsevier, New York, 449468. Z. Chen (1993b). A SemiAnalytical Solution Procedure for Predicting Damage Evolution at Interfaces. International Journal for Numerical and Analytical Methods in Geomechanics 17, 807819.
124
Chen, Z. et al.
Z. Chen (1996a). Continuous and Discominuous Failure Modes. Journal of Engineering Mechanics 122, 8082. Z. Chen (1996b). A Simple Procedure to Simulate the Failure Evolution. International Journal of Structural Engineering and Mechanics 4, 601612. Z. Chen and H.L. Schreyer (1990a). A Numerical Solution Scheme for Softening Problems Involving Total Strain Control. Computers and Structures 37, 10431050. Z. Chen and Schreyer, H.L. (1990b). Formulation and Computational Aspects of Plasticity and Damage Models with Application to QuasiBrittle Materials, SAND950329, Sandia National Laboratories, Albuquerque, NM. Z. Chen and H.L. Schreyer (1991). Secant Structural Solution Strategies under Element Constraint for Incremental Damage. Computer Methods in Applied Mechanics and Engineering 90, 869884. Z. Chen and H.L. Schreyer (1994). On Nonlocal Damage Models for Interface Problems. International Journal of Solids and Structures 31, 12411261. Z. Chen and D. Sulsky (1995). A PartitionedModeling Approach with Moving Jump Conditions for Localization. International Journal of Solids and Structures 32, 18931905. Z. Chen and T. Clark (1995). Some Remarks on DomainTransition Problems. Archives of Mechanics 47, 499512. Z. Chen and X. Xin (1997). An Analytical and Numerical Study of Failure Waves. To appear in International Journal of Solids and Structures. Chen, Z., X. Xin and D. Qian (1997). A Study of Localization Problems Based on the Transition between Governing Equations. Submitted for publication in International Journal of Mechanics of CohesiveFrictional Materials. G. PijaudierCabot and L. Bode (1995). Arbitrary LagrangianEulerian Finite Element Analysis of Strain Localization in Transient Problems. International Journal for Numerical Methods in Engineering 38, 41714191. D. Sulsky, Z. Chen and H.L. Schreyer (1994). A Particle Method for HistoryDependent Materials. Computer Methods in Applied Mechanics and Engineering 118, 179196. D.Y. Tzou (1989). On the Thermal Shock Wave Induced by a Moving Heat Source. ASME Journal of Heat Transfer 111,232238. M. Xie, W. Gerstle and Z. Chen (1994). Finite Element Analysis of Combined Smeared and Discrete Mechanisms in Rock Salt. Computer Methods and Advances in Geomechanics, edited by H.J. Siriwardane and M.M. Zaman, A.A. Balkema Publishers, Brookfield, VT, pp. 16591664.
CCONTINUATION APPROACH FOR TRUSS TOPOLOGY OPTIMIZATION
Xu Guo and Gengdong Cheng State Key Laboratory of Structural Analysis of Industrial Equipment Dalian University of Technology Dalian 116023, PR CHINA
ABSTRACT A socalled ccontinuation approach is proposed for the solution of singular optima in truss topology optimization. In the approach, we start the optimization process from a relatively large value of c, apply the crelaxed formulation and obtain an optimum design of this relaxed problem; then we decrease the value of c by Ae, and choose the design from the previous optimization as initial design to begin the next optimization, and continue the process until convergence. It is shown by numerical examples that this scheme alleviates the dependence of the final solution on the choice of the initial design and increases greatly the probability of finding the singular optima from rather arbitrary choice of the initial design.
KEYWORDS structural topology optimization; numerical algorithm; singular optima; relaxed formulation;
INTRODUCTION While many efforts have been made in solving topology optimization problems of truss structures, it is wellrecognized that there are still serious difficulties in this field. One of them is the socalled singular optima phenomenon which prevents the iterative algorithms from
125
126
Guo, Xu and Cheng, Gengdong
converging to the true optimal solution. Singular optima in structural topology optimization was first shown by Sved and Gions t11.They found that in some cases, global optimal topology could not be obtained by employing the nonlinear programming techniques unless one deleted the unnecessary bars from the ground structure in advance. Since then, the phenomena of singular optima was studied intensively by many authors. Kirsch t2j investigated the singular optima in truss and grillage topology optimization problems. He suggested that singular optima corresponds to a singular point in the design space and it is very difficult or even impossible to arrive at it by numerical search algorithms. Cheng and Jiang I31pointed out that singular optima appears mainly due to the nature of local behaviour constraint. Illustrating by truss topology optimization problems subjected to stress constraints, they demonstrated that the discontinuity of the bar's stress constraint function when its cross sectional area takes zero value is the essential cause of the existence of singular optima. The authors suggested to establish a rational formulation which unifies the sizing and topology optimization and enables one to apply the sizing optimization techniques for the solution of topology optimization problems. Rozvany t41 studied the singular optima in the light of exact optimal layout theory and proposed a precise definition of singular optima. Since singular optima is a major obstacle to structural topology optimization, various methods have been suggested to overcome the difficulties. Kirsch t51 proposed a design procedure to find the optimum of singular topologies by neglecting the compatibility conditions and applying linear programming technique. His method is appropriate only when the optimal design is statically determined. Recently, Kirsch t6j has presented another twostage topology optimization procedure consisting of reduction and expansion process in order to eliminate the problems of singular optima which may appear in the conventional approach. RozvanyETJ suggested that for truss topology optimization problems subjected to stress constraints ,smooth envelope functions should be used to make the feasible set nonsingular. Recently, a erelaxed approach for the solution of singular optima has been developed by Cheng and Guo t81. In this approach, the problem is reformulated by relaxing the local behaviour constraint through the introduction of a relaxation parameter c. Under this relaxed formulation, singular optima does not exist in the design space, so that one can now apply the sizing optimization technique to solve this reformulated problem. It is proved that the solution of the relaxed problem is a good approximation of the solution of the original problem as long as c is sufficiently small. Although great successes have been achieved by applying this approach for the solutions of singular problems, however, it should be noted that this approach cannot locate the global optimal solution exactly unless the relaxation parameter e is small enough and initial design is properly chosen. Hoback t91proposed a percent method to find the global optimal design of singular problems and presented a number of numerical examples. In the present paper, based on the aforementioned crelaxed approach, a socalled ccontinuation approach is developed for the solutions of truss topology optimization. By applying this approach, in most of the test cases, global optimal solution, even if a singular
eContinuation Approach for Truss Topology Optimization
127
one, can be obtained from rather arbitrary choice of the uniform initial design. The effectiveness of this approach is illustrated by several numerical examples.
oeRELAXED A P P R O A C H For completeness, in this section, we will describe the erelaxed approach briefly. The problem of truss topology optimization subject to stress constraints can be given as follows. AI,A2,...,AN Po" to find N
Min W = ~p,l~A,
(la)
i=1
St
KUj = Pj
j = 1,..., m
O'iL   % ~~0
j = 1,..., M ;
i = 1,2,...,N
(lc)
cr~oy < 0
j = 1,..., M ;
i= 1,2,...,N
(ld)
i = 1,2,...,N
(le)
(lb)
A, >_0
As noted in [8], under crelaxed approach, the problem can be reformulated as follows: P• to find A~,A2,... ,A u N
Min W = ~ p,l, At
(2a)
i=1
St
KUj  Pj
(2b)
j = 1,...,M
(o~L  cr~) A, < c
j = 1,..., M ;
(or0  o'y ) At < e
j = 1,...,M;
i = 1,2,..., N i = 1,2,...,N
(2c) (2d)
A~L
(2e) i = 1,2,..., N where N is the total number of bars in the ground structure, M denotes the number of loading cases. At ,p, ,l~ are cross sectional area, mass density and the length of ith bar. The superscript A i ~_
U and L refer to the upper and lower bounds of the relevant quantities, respectively. Under this formulation, the stress constraints are relaxed by introducing a relaxation parameter e. It can be seen that for any given ~ > 0, the constraint (2c),(2d) can always be satisfied for sufficiently small A t . Thus, in the vicinity of the degenerated subdomain corresponding to Ai = 0, the measure of the feasible design space is nonzero. Therefore, under this formulation, optima is not singular and one can apply the conventional numerical search algorithms to arrive at the optima of the relaxed problem. With the help of the theory of point to set mapping tl01, it can be proved that the optimal solution of the relaxed problem will converge to the solution of the original problem as the relaxation parameter c tends to zero. Thus the optimal solution of the reIaxed problem provides a good approximation of the corresponding optima of the original problem, even if it is singular.
128
Guo, Xu and Cheng, Gengdong
eCONTINUATION APPROACH Although the erelaxed approach enables one to reach the singular optima by numerical optimization, it should be noted that just like most of the optimization techniques, this approach cannot locate the global optimal design unless the initial design is properly chosen. In fact, the global optimization is still a challenge in the field of structural optimization. From mathematical point of view, it is essentially unsolvable. In order to obtain the global optimal design, various methods haye been proposed. But unfortunately, these methods are either eomputationally expensive or not suitable for the solution of largescale problems. Thus for a practical problem, if global optima is needed, it seems that the best thing one can do is to extend the "attractive region" of the global optima if one is not willing to bear the enormous cost. The econtinuation approach follows this idea. Its basic steps can be described as follows. Step 1. Under erelaxed formulation, start the optimization process from a relatively large value of c and an arbitrary chosen uniform initial design; Step 2. Decrease e by Ae, and choose the optimal solution obtained from the previous optimization step as initial design to begin the next optimization step; Step 3. Continue this until a small desired value of e has been reached.
G
S
Figure 1" Feasible domain for singular problem To explain the idea behind the approach, let us consider a typical geometric representation of singular topology optimization (see Fig. 1). There is a degenerated subdomain in the design space and the global singular optima is located at point G in whose vicinity the measure of the feasible domain is zero. Under conventional formulation, numerical search algorithms terminate at nonoptimal solution point S.
Figure 2" Feasible domain under different values of e
eContinuation Approach for Truss Topology Optimization
129
Under erelaxed formulation, the design space for different values of c can be shown in Fig.2. Observing these figures, it is clear that under erelaxed formulation, the shape of the feasible domain is rather regular for a large e; with the decreasing of the value of e, the feasible domain of the relaxed problem approaches more and more that of the original problem, but becomes irregular. Since the shape of the feasible domain is rather regular for a large e , the attractive region of the global optima of the relaxed problem may also be relatively large; therefore the probability with which numerical search algorithms converge to the global optima from an arbitrary choice of the initial design is high. Moreover, based on the convergence analysis of erelaxed approach in [8] it has been concluded that the global optimal design of the relaxed problem is a continuous function o f t , so the optimal designs from two successive optimization steps may be very close if the degradation of e is not too large. Thus the optimal solution obtained from previous optimization step is a good initial design for the next optimization step. In other words, by means of this continuation approach, we can choose a proper direction to converge to the global optimal solution of the original problem. Numerical examples show that in most of the cases, the global singular optima can be obtained by means of this approach without much extra computation effort.
CONVERGENCE ANALYSIS In this section, proof of convergence for our econtinuation approach will be given. It is established under the following assumptions. A s s u m p t i o n 1 . The original optimization problem P0 has a bounded optimal solution Wgpt < +00 .
A s s u m p t i o n 2. In every optimization step, we can get the global optimal solution WS of the
relaxed problem P . Denoting the feasible set of the design variable associated with problem P0 and problem P~ as U 0 and Ue respectively, we have the following lemmas. Lemmal.
Ue, cUe, for ck _0,ct _>0.
The proof of this lemma is obvious. L e m m a 2 . I f c k 0 , C t > 0 thenW ° p ' > W °p' '


~
ek

el
"
Proof: Denoting the optimal design corresponding to problem P~ as A~p' , then from lemma 1
we
W "p' = W ( A °p`) el
cl
have 
A°P' c Uek c Ue,
if
e k
°p' • inf W(A) _<W(A~p') = We~
A ~U~t
In
opt
W °p' < Wo"p' = W(A 0 ) for all ck > 0 ~'k
m
~
*
From the above lemmas, we have the following theorem.
particular,
ck < c/ . we
Then have
130
Theorem
G u o , X u and Cheng, G e n g d o n g
1: For a positive sequence
assumptions 12, we have k.+~ lim W~'°p' ~
{ek}, ~k ~>Ek+l > 0 and k~oo limc k = 0; then under
~ W;pt"
Op` }is an above bounded Proof: By lemma 2 and assumption 1, it is clear that ,¢,Wck
monotontically increasing sequence. Therefore, from the theory of mathematical analysis, the sequence (WOp , ~k, } must have limit lim W °p' = W < Wo p' For the completeness of the proof, it is k~
6t


"
only required to show that W < Wo p' is impossible. Assuming that W < Wo p' , then we can find a real number 6 > 0 such that 0 = Wo p'  W  6 > O. Since lim W °p' = W k~
integer K o > 0, for all gk such that k > K o, ;Iwo ~, WOekpt < W  F O  W ; p t  W 
Ek
there exists an
< 0", then for all k > K o (~I W  Wgpt  (~
(3)
On the other hand, based on the continuity result of W °p' (e) with respect to c when c = 0 which was established in [8], we have for 6 > 0; there exists an integer K~ .> 0, such that for all k > K1, ,[W°p`'~  W / ' [ < 6 . So it readily follows that for all k > K 1, W °p`'` > Wo p'  6 . contradicts
This
(4a) when k is sufficiently large. Thus the assumption of W < Wo p' is invalid
and W = Wo p' holds. The proof is completed.
N U M E R I C A L EXAMPLES In this section, several numerical examples will be presented to illustrate the effectiveness of our econtinuation approach. A program which implements the Constraint Variable Metric optimization technique was used as optimizer and another structural analysis program has also been developed for the function calculations. All the design variables and constraints are scaled. Sensitivity analysis is carried out by means of backward finite difference approach. The difference step was T=0.005 for all examples. Generally speaking, the starting value of ~ may be chosen'heuristically as long as g0 is larger than a critical value %. It is clear that this value is problemdependent. Setting %=10.0 was found to be adequate for all examples shown below. The value of A~ also depends on the characteristics of the particular problem; setting ek+~=ek0.5 often works well. The termination value ~ter can be taken as 0.010.00001. If ~ter is too small, machine error will predominate; thus the effect of relaxation can not be manifested. The symbol SF,IF,RF in the following text denotes the stress formulation (1), internal force formulation which impose constraints on bar internal force, i.e., set e=0 in (2c) and (2d), and ~relaxed formulation (2), respectively.
eContinuation Approach for Truss Topology Optimization
131
3bar truss example A threebar truss subjected to three load cases (see Fig.3)is optimized. It is a classical example from [ 1]. E= 1.0,h= 1.0,p= 1.0. The allowable stresses are oia = +5.0,i = 1,3; o2a = +20.0. The three loading cases are (a) P=40,c=45°; (b) P=30, c=90 ° ; (c) P=20,c=135 °, respectively.
h
h
Figure 3" Sved and Gions' 3bar truss example As was shown t~], its global optimal design is a 2bar truss, i.e., A~ = 8.00, A2 = 1.50, A3 = 0.00. This is a singular optima because for this design lim o23 = 20.00 > 5.00, lim o'33 = 21.36 > 5.00 The optimization results obtained under different problem formulations are listed in Table 1. TABLE 1 OPTIMIZATION RESULTS UNDER DIFFERENT FORMULATIONS
Final results of different formulations from different initial designs Relaxed Formulation Initial Design Stress Formulation A1 A2 A3 W A1 A2 A3 2.4 7.024 2.135 2.767 15.97 W15.97 3.2 8.0 12.0 18.0 0.000 2.828 8.000 14.14 W=15.97 0.0 2.828 W=15.97 10.0 0.000 8.000 14.14 0.7 4.0 7.024 2.135 2.767 15.97 W=15.97 1.0 1.0 1.0 w215.97 10.0 7.024 2.135 2.767 15.97 1010 10.0 w=15.97 1.50 0.00 12.81 1.0 1.0 0.1 8.00 w=15.97 1.50 0.00 12.81 10.0 0.1 8.00 10.0 w=15.97 1.50 0.00 12.81 8.00 G1.Opt From Table 1, it can be seen that under conventional stress formulation, the singular global optimal solution cannot be obtained no matter how you choose the initial design. The numerical search algorithms always terminate at a local optima. But if the initial design is properly chosen, it is possible for crelaxed approach to arrive at the global optimal solution. It is worth noting, however, that if we start the optimization process from an arbitrary chosen initial design, crelaxed approach cannot always locate the global optima. For this example,
132
Guo, Xu and Cheng, Gengdong
the global optima cannot be achieved unless the initial value of A3 is less than 0.1. This unpleasant behavior can be attributed to the fact that when the value of the relaxation parameter 6 is very small, the shape of the feasible domain may be rather irregular in the vicinity of the singular optima; thus the attractive region of the singular optima may be very small, so the global optima can be found only when the initial design is very close to it. TABLE 2 OPTIMIZATION HISTORY AND CORRESPONDING ELEMENT STRESS
Continuation History (3' 10.0) A1 A2 A3 6.8409 0.8214 0.2406
C 5.0
4.0
7.0232 .
.
.
0.9554 .
.
.
.
.
.
Element Stress W 10.684
0.2075
11.094
.
2.0
7.4931
1.2270
0.1149
11.937
1.0
7.7425
1.3635
0.05"95
12.371
0.5
7.8703
1.4318
0.0317
12.592
0.1
7.9739 ] 1.4864
0.0061
12.768
7.9974
0.0006
12.808
0.01 .
.
.
.
.
1.4986
.
0.001
8.0000
1.4999
0.0000 [ 12.812
Gl.Opt
8.0000
i.5000
0.0000
1
2
5.735 0.891 2.028 5.570 0.694 2.130 5.267 0.327 2.326 5.129 0.158 2.416 5.062 0.082 2.454 5.013 0.015 2.492 5.001 0.002 2.499 5.000 0.000 2.500
1.680 26.087 23.756 1.309 24.187 22.150 0.616 21.630 20.085 0.298 20.733 19.399 0.154 20.319 19.080 0.029 20.067 18.905 0.003 20.007 18.860 0.000 20.001 18.857
....
,,,
3 4.056 25.197 25.,783 4.261 23.493 24.280 4.651 21.303 22.410 4.831 20.575 21.814 4.909 20.237 21.534 4.983 201052 21.397 "4.998 20.005 21.360 5.000 20.000 21.357
12.812
To illustrate the effectiveness of our cdontinuation approach, optimizations are performed by means of the propose method under different choice of the initial design. Since in practice, we have no idea about which bar will remain in the final optimal topology, it is natural to
e  C o n t i n u a t i o n A p p r o a c h for Truss T o p o l o g y Optimization
133
start the optimization process from a uniform initial design. Taking the initial design as 3* 1.0, 3*5.0, 3'10.0, 3*20.0, in all cases, ccontinuation approach gives the same optimal result. Table 2 gives the optimization results from initial design 3'1.0 with e0  5 . 0 . With reference to Table 2, it may be seen clearly that the proposed approach could achieve result very close to the global optima of the original problem when e is sufficiently small. It can also be observed that with the decreasing of e, all of the constraints can be satisfied to a high accuracy. In this example, the termination value of e is taken as 0.001; further iterations did not improve the result. Fig.4 shows the continuation history of the structural weight. W
12.5
11.5
10.5
0 o.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Figure 4" Continuation history of structural weight of 3bar truss
5bar truss example Fig.5 shows the ground structure. Geometrical and material data are E = 1.0,h = 1.0,p = 1.0, allowable stresses are 0% = +5.0,i = 1,3,4,5; a2a = +20.0. Two loading cases are considered; they are (a) Pax=0.0, P~y=0.0, Pbx=5.0, Pby=50.0; (b) P~ =5.0, P~y=50.0, Pbx=0.0 P b y = 0 . 0 ; t h e optimization results of this example are given in TABLE 3. The global optimal solution of this problem is a still a singular one because for the optimal design we have
llimcr,, I= 15.00 > 5.00, Ilimcr2,I= 7.25 > 5.00
Y
h
/ X
Figure 5" 5bar truss example
134
Guo, Xu and Cheng, Gengdong
Inspection of Table 3 reveals that the conventional formulation generates only nonoptimal design no matter from what initial design the iteration starts. For econtinuation approach, global optimum was obtained. Table 3 gives the optimization results with 5"10.0 and eo = 10.0 . Global optimal design is obtained after 7 continuation steps. TABLE 3 OPTIMIZATION RESULTS FOR 5BAR TRUSS
10.0 8.0 610 4.0
A1 0.776 1.12 1.11 1.07
1.0 0.01 0.001 GL.Opt SF
Continuation History(5* 10.0) A2 A3 A4 1.55 6.23 0.97
1.01
1.67 1.87 2.08 2.40
6.68 7.50 8.33 9.59
0.84 0.64 0.43 0.10
1.00 1.00 1.00 7.219
2.49 2.50 2.50 4.735
9.96 10.00 10.00 4.531
0.01 0.00 0.00 8.795
A5 11.06 11.51
W 26.933 28.551
12.14 12.79 13.80
30.173 32.657 33.416
14.14 14.14 14.14 7.828
33.486 33.496 33.500 39.990
Hoback's 4bar truss example The 4bar truss shown in Fig.6 was optimized. This example was discussed by Hoback for demonstrating the effectiveness of his percent method for the solution of singular problems. The structure is subjected to two load conditions, that is, P t=I 0, P2=0 for load case 1, P 1=0, P2=10 for load case 2. The allowable stresses for all members are 20 and 5 in tension and
compression respectively. 3
4 25
~75
" 100 ~  " 1 0 0 r Figure 6" Hoback's 4bar truss example Its global optimal solution is Al = 0.000,A 2 = 0.4286,A 3 = 0.4773 and A4 = 0.3476. This is a singular and static indeterminate design because for this design the limiting stress of bar 1 violates the constraint imposed on it. The results of econtinuation approach are shown in
eContinuation Approach for Truss Topology Optimization
135
Table 4. Singular optima was obtained after 7 continuation steps from the uniform initial design 4* 10.0 (see Table 4). Optimization from other uniform initial design gives the same results. Under conventional formulations, numerical algorithms iterate only to the local optimum. Compared with the percent method tgj, econtinuation approach can obtain the global optimal solution from arbitrary choice of the initial design without introducing extra design variables. And our result is lighter in weight than the one given by the percent method. TABLE 4 CONTINUATION HISTORY, FINAL DESIGNS FROM HOBACK, SF AND IF FOR 4BAR TRUSS
2.0 1.5 1.0 0.5 0.1 0.01 0.001 Hoback SF IF
A1 0.0529 0.0481 0.0372 0.0210 0.0046 0.0000 0.0000 0.0000 1.0157 0.tJ000
A2 0.1451 0.2083 0.2757 0.3480 0.4115 0.4274 0.4286 0.5600 1.0115 1.2500
A3 0.1521 0.2586 0.3361 0.4109 0.4653 0.4757 0.4773 0.3980 0.4939 0.0000
A4 0.1592 0.2021 0.2483 0.2969 0.3373 0.3471 0.3476 0.3620 0.4886 2.1479
W 79.39 105.79 132.66 159.42 180.48 185.15 185.62 191.00 410.69 500.08
,,,
10bar truss example This example is a modified version of the wellknown 10bar truss problem(see Fig. 7). The truss is subjected to two load cases, that is, P~=100, P2=0 for load case 1, P1=0, P2=200 for load case 2. E = 1.0E + 04,h = 360.0, p = 1.0. The allowable stress for all members is 20 and 5 in tension and compression respectively.
.
Pl
Figure 7:1 Obar truss example To show the superiority of our ~continuation approach to other methods; the example was optimized under different formulations.. The results are listed in Table 5. From the uniform initial design 10* 10.0, ~continuation approach gives the global optimal solution denoted by GO whereas algorithms under stress formulations can only arrive at the solution which is non
136
Guo, Xu and Cheng, Gengdong
optimal. In Table 6 we list the bar stresses under the two load cases for the two designs obtained by the econtinuation approach and the stress formulation. It shows that the solution from stress formulation is not the optimum one because for this set of cross sectional area, the stress constraint imposed on bar 9 is active though this bar does not exist in the final structure. TABLE 5 OPTIMIZATION RESULTS FOR 10BARTRUSS , m
A1 A2 A3
10.0 10.09 8.91 75.66
Continuation History (initial 6.0 1.0 0.1 10.07 9.99 10.00 9.33 9.90 9.99 77.33 79.64 79.96
design 10* 10.0) 0.01 SF 10.00 10.00 10.00 10.59 80.00 81.24
A4 A5 A6 A7 A8 A9 A10 W
1.04 0.00 3.74 12.81 2.29 1.17 51.25 7017
0.67 0.00 4.23 13.32 1.44 0.71 53.29 7160
0.00 0.00 5.00 14.14 0.00 0.00 56.56 7380
0.08 0.00 4.90 14.03 0.09 0.12 56.14 7346
0.01 0.00 4.99 14.13 0.02 0.01 56.52 7376
.......
,,,
0.02 0.00 11.30 22.54 0.04 0.00 96.62 10142
IF 11.43 8.65 74.28
GO 10.00 10.00 80.00
6.06
0.00 0.00 5.00 14.14 0.00 0.00 56.56 7380
0.11
5.40 12.47 9.18 2.14 48.94 7517
m
TABLE 6 BAR STRESS FOR FINAL DESIGNS UNDER DIFFERENT FORMULATIONS
Bar Stress C~ SF
Bar Load 1 Load 2 Load 1 Load 2
1
2
3
10.00 10.00 3.75 20.00 20.00 5.00 10.00 9.44 3.69 20.00 18.88 4.92
4
9.59 8.00 0.00 0.00
6 1.73 20.00 0.62 0.00 18.78 8.85 10.02 0.00 5
7 8 20.0016.01 20.00 12.19 12.55 0.00 12.55 0.00
9 13.57 11.31 20.00 17.38
10 2.50 5.00 1.46 2.93
RANDOM TEST OF THE ECONTINUATION APPROACH The present econtinuation approach is a heuristic algorithm. It cannot guarantee the convergence to global optimal solution from any initial design. But as we know, a good heuristic algorithm should locate the global optimal solution with high probability. To study the practical performance of our econtinuation approach for finding the singular global optimal solutions, in this section, a set of random tests are performed by applying the e
~Continuation Approach for Truss Topology Optimization
137
continuation approach together with Constraint Variable Metric optimization technique with different randomly generated initial designs and checking whether the global optimal solution of the problem is attainable.
3bar truss example Let us consider the above mentioned 3bar truss example. Using the random number generator, the sample space SNU(nonuniform initial design) and SU(uniform initial design) of the test are generated respectively. Table 7 lists the optimization results obtained from different sets of nonuniform initial designs and optimization algorithms. With 20 randomly chosen initial uniform designs,algorithms all converge to the global optima. Here symbol S denotes the fact that the singular global optima can be obtained starting from the corresponding set of initial design, whereas symbol F denotes that only local optima is achieved with this set of initial design. TABLE 7 TESTS FOR 3BAR TRUSSWITHNONUNIFORMINITIALDESIGN A°
A°
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
8.91 4.30 5.35 9.00 5.29 4.36 ..... 8.92 6.20 3.30 8.69 7.01 2.20 8.30 7.69 1.07 7.77 8.24
A°
6.85 3.59 9.65 7.76 2.35 9.22 8.53 1.08 8.62 9.15 2.19 7.86 9.60 3.44 6.97 9.89 4.64
1.09 9.46 10.06 2.32 8.62 10.56 3.71 7.64 10.88 5.04 6.51 11.00 6.29 5.28 10.92 7.44 3.96
S S S S S S S s S S S F S S F s s
18 19 20
2.06 7.10 8.65
5.96 9.99 5.76
10.64 8.45 2.58
S S S
"
c
S F F S F F F F F F F F F F F F F F F S
138
Guo, Xu and Cheng, Gengdong
Hoback's 4bar truss example In this example, 20 nonuniform and 20 uniform random generated initial designs are used to perform the test. The optimization results corresponding to each of the nonuniform initial design are listed in Table 8. Again, with 20 randomly chosen initial designs, algorithms all converge to the global optima.
TABLE 8 RESULTS OF 4BAR TRUSS WITH NONUNIFORM INITIAL DESIGN A °
A°
A°
A°
1 2 3
12.00 6.90 5.62
11.04 1.10 11.15
6.25 8.17 14.00
3.22 13.83 12.64
F S S
F F F
4
11.90
11.87
7.88
1.25
F
F
5 6
8.15 4.17
2.59 10.15
6.57 13.89
12.91 13.63
S S
F F
7 8
11.58 9.26
12.48 4.25
9.37 4.85
2.73 11.76
F F
F F
9 10 11 12 13 14 15 16
2.65 11.04 10.20 1.10 10.31 10.96 2.46 9.39
8.96 12.86 5.85 7.62 13.00 7.35 6.14 12.90
13.53 10.69 3.06 12.91 11.81 1.23 12.06 12.72
14.35 4.67 10.38 14.81 6.54 8.82 15.00 8.30
S F S S F S S S
S F F S F F S F
17 18
11.52 3.98
8.72 4156
2.60 10.99
7.10 14.90
S S
F F
19 20
8.30 11.87
12.57 9.94
13.40 4.41
9.91 5.27
S S
F F
,,,
5bar truss example For this example, the results of the performed optimizations with 20 nonuniform initial designs are given in Table 9. With 20 randomly chosen initial designs except the one A ° , A° , A° = 1.64, all of them converge to the global optima.
eContinuation Approach for Truss Topology Optimization
139
TABLE.9 RESULTS OF 5BAR TRUSS WITH NONUNIFORM INITIAL DESIGN
1
A°
A°
A°
A°
Ag
8.91
6.85
1.09
8.30
12.90
C
F
F
/
2
4.30
3.59
9.46
11.87
8.72
F
F
3
5.35
9.65
10.06
5.44
4.56
S
F
4
9.00
7.76
2.32
7.07
12.57
S
F
5
5.29
2.35
8.63
12.00
9.94
F
F
6
4.36
9.21
10.56
6.82
2.90
F
F
7
8.92
8.53
3.71
5.71
12.00
S
F
8
6.20
1.08
7.64
11.91
10.98
F
F
9
3.30
8.62
10.88
8.08
1.21
F
F
10
8.69
9.15
5.04
4.26
11.21
S
F
11
7.01
2.19
6.51
11.60
11.82
F
F
12
2.20
7.86
11.00
9.20
2.48
F
F
13
8.30
9.60
6.29
2.74
10.22
S
S
14
7.69
3.44
5.28
11.08
12.45
1~ 
F
15
1.07
6.97
10.92
10.15
4.15
S
F
16
7.77
9.89
7.44
1.19
9.04
S
S
17
8.24
4.64
3.96
10.36
12.84
F
F
18
2.06
5.96
10.64
10.92
5.75
S
F
19 20
7.10 8.65
10.00 5.76
8.45 2.59
2.36 9.45
7.71 13.00
S S
S F
With reference to the above results, the following observations are drawn: ~continuation approach has the high probability to locate the singular global optima starting from rather arbitrary uniform initial design. It is worth noting that starting the optimization process from uniform initial design is natural if we have no idea about the optimal topology in advance. Because of this feature, the method provides an ideal tool for singular topology optimization problems where the singular global optimum are usually very difficult to achieve by other methods. For some nonuniform initial design, Econtinuation approach can lead to the singular optima, though the convergence to singular global optima cannot be guaranteed. However, it is well known that global optimum remains open for nonconvex mathematical programming.
140
Guo, Xu and Cheng, Gengdong
CONCLUSIONS For making the most efficient use of resources, in engineering practice, the determination of global optima is always highly desirable. By this motivation, a econtinuation approach which is a variant version of erelaxed approach is proposed in this paper. It is based on the solution of a sequence of problems constructed by relaxing the constraints to some extent through the introduction of a relaxation parameter e. In the limit as e tends to zero, the solution of the singular problem is approached. This approach works well for finding the global optimal solution of singular problems. By means of this approach, global optimal solution can be obtained from rather arbitrary choice of the initial design; thus the initial designdependency problem of the erelaxed approach is alleviated. Though this approach cannot guarantee a global optimum from any initial design, numerical examples show that it does have the ability to avoid entrapment by local optima and locate the global optima with high probability. Although only truss topology optimization problems under stress constraints are discussed in this paper, the method presented can be applied equally to other discrete topology optimization problems with various behavior constraints. REFERENCES
1. .
.
.
.
.
7. 8. 9. 10.
Sved G. and Gions Z.(1968). Structural optimization under multiple loading, Int. J. Sci., 10,803805 Kirsch U.(1989) Optimal topologies of truss structures, Comput. Methods. Appl. Mech. Eng., 72(1),1528 Cheng G. and Jiang Z.(1992), Study on topology optimization with stress constraint, Eng. Optim., 20, 129148 Rozvany G.I.N. and Birker T.(1994), On singular topologies in exact layout optimization, Struct. Optim., 8,228235 Kirsch U.(1990), On singular topologies in optimum structural design, Struct. Optim., 2,3945 Kirsch U.(1996), Integration of reduction and expansion process in layout optimization, Struct. Optim., 11, 1318 Rozvany G.I.N.(1996), Difficulties in truss topology optimization with local buckling and system stability constraints, Struct. Optim., 11,1318 Cheng G. and Guo X.(1997), ERelaxed approach in structural topology optimization, Struct. Optim., 13(4) Hoback A.S.(1996), Optimization of singular problems, Struct. Optim., 12, 9397 Hogan W.W.(1973), Point to set maps in mathematical programming, SIAM review, 15, 591603
RECENT DEVELOPMENTS IN BASIC FINITE ELEMENT TECHNOLOGIES Carlos A. Felippa Department of Aerospace Engineering Sciences and Center for Aerospace Structures University of Colorado, Boulder, Colorado 803090429, USA
ABSTRACT High performance (HP) elements are simple finite elements constructed to deliver engineering accuracy with arbitrary coarse meshes. This paper overviews original developments based on parametrized variational principles, which provide a common foundation for various approaches to HP element construction. This work led to the template approach. Templates are parametrized algebraic forms of finite element stiffness equations. The application of templates to the construction of an optimal 3node plate bending element is outlined. The template approach holds future promise in the unification of finite element methods and connection to arbitrary grid finite differences. Essential to this unification work is the availability of powerful computer algebra systems. KEYWORDS
Finite element method, variational principles, templates, parametrization, morphing, high performance elements, optimal elements, plate bending, triangular elements, computer algebra systems. INTRODUCTION
The Finite Element Method (FEM) was first described in its presently dominant form by Turner et al. (1956). It was baptized by Clough (1960) at the beginning of an explosive growth period. The first applications book, by Zienkiewicz and Cheung (1967), appeared seven years later. The first monograph on the mathematical foundations was written by Strang and Fix (1973). The opening sentence of this book already declared the FEM an "astonishing success." And indeed the method had by then revolutionized computational structural mechanics and was in its way to impact nonstructural applications. The FEM was indeed the right idea at the right time. The key reinforcing factor was the expanding availability of digital computers. Lack of this enabling tool meant that earlier related proposals,
141
142
Felippa, C.A.
notably that of Courant (1943), had been forgotten. A second benefit was the heritage of classical structural mechanics and its reformulation in matrix form, which culminated in the elegant unification of Argyris and Kelsey (1960). A third influence was the victory of the Direct Stiffness Method (DSM) developed by Turner (1959, 1964) over the venerable Force Method, a struggle recently chronicled by Felippa (1995). Victory was sealed by the adoption of the DSM in the earlier generalpurpose FEM codes, notably NASTRAN, MARC and SAP. In the meantime the mathematical foundations were rapidly developed in the 1970s. "Astonishing success", however, carries its own dangers. By the early 1980s the FEM began to be regarded as "mature technology" by US funding agencies. By now that feeling has hardened to the point that it is virtually impossible to get significant research support for fundamental work in FEM. This viewpoint has been reinforced by major software developers, which proclaim their products as solutions to all user needs. Is this perception correct? It certainly applies to the core FEM, or orthodox FEM. This is the material taught in textbooks and which is implemented in major software products. Core FEM follows what may be called the RitzGalerkin and Direct Stiffness Method canon. Beyond the core there is an evolving FEM. This is strongly rooted on the core but goes beyond textbooks. Finally there is a frontier FEM, which makes only partial or spotty use of core knowledge. By definition core FEM is mature. As time goes, it captures segments of the evolving FEM. For example, most of the topic of FEM mesh adaptivity can be classified as evolving, but will eventually become part of the core. Frontier FEM, on the other hand, can evolve unpredictably. Some components prosper, mature and eventually join the core, some survive but never become orthodox, while others wither and die. Four brilliant contributions of Bruce Irons, all of which were certainly frontier material when first published, can be cited as examples of the three outcomes. Isoparametric elements and frontal solvers rapidly became integral part of the core technology. The patch test has not become part of the core, but survives as a useful if controversial tool for consistent element development and testing. The semiloofshell elements quietly disappeared. The words "basic FEM technologies" in the title of this article refers to construction of elementlevel models of physical problems. Frontier topics therein include interface, multiscale, and high performance elements. The material that follows deals with the last topic. HIGH PERFORMANCE ELEMENTS
An important objective of frontier FEM is the construction of high performance (HP) finite elements. These have been defined by Felippa and Militello (1989) as "simple elements that deliver engineering accuracy with arbitrary coarse meshes?' This definition requires further clarification.
Simple means the simplest geometry and freedom configuration that fits the problem and target accuracy consistent with human and computer resources. This can be summed up in one FEM modeling rule: use the simplest element that will do the job. Engineering accuracy is that generally expected in most FEM applications in Aerospace, Civil and Mechanical Engineering. Typically this is 1% in displacements and 10% in strains, stresses and derived quantities. Some applications, notably in Aerospace, require higher precision in quantities such as natural frequencies, shape tolerances, or in longtime simulations.
Recent Developments in Basic Finite Element Technologies
143
Coarse mesh is one that suffices to capture the important physics in terms of geometry, material and load properties. It does not imply few elements. For example, a coarse mesh for a fighter aircraft undergoing maneuvers may require several million elements. For simple benchmark problems such as a uniformly loaded square plate, a mesh of 4 or 16 elements may be classified as coarse. Finally, the term arbitrary mesh implies low sensitivity to skewness and distortion. This attribute is becoming important as pushbutton mesh generators gain importance, because generated meshes can be of low quality compared to those produced by an experienced analyst.
Tools for Construction of HP Elements The origins of HP finite elements may be traced to several investigators in the late 1960s and early 1970s. Notable early contributions are those of Clough, Irons, Taylor, Wilson and their coworkers. The construction techniques made use of incompatible shape functions, the patch test, reduced, selective and directional integration. These can be collectively categorized as unorthodox, and in fact were labeled as "variational crimes" at that time by Strang and Fix (1973). A more conventional development, pioneered by Pian, Tong and coworkers, made use of mixed and hybrid variational principles. They developed elements using stress or partial stress assumptions, but the end product were standard displacement elements. This approach was further refined in the 1980s. A good expository summary is provided in the book by Zienkiewicz and Taylor (1992). New innovative approaches came into existence in the 1980s. The most notable have been the Free Formulation of Bergan and Nyghrd (1984), and the Assumed Strain method pioneered by MacNeal (1978). The latter was further developed along different paths by Bathe and Dvorkin (1985), Park and Stanley (1986), and Simo and Hughes (1986).
Unification by Parameterized Variational Principles The approach taken by the author started from collaborative work with Bergan in Free Formulation (FF) high performance elements. The results of this collaboration were a membrane triangle with drilling freedoms described in Bergan and Felippa (1985) and a plate bending triangle described by Felippa and Bergan (1987). It continued with exploratory work using the Assumed Natural Strain (ANS) method of Park and Stanley (1986). Eventually FF and ANS coalesced in a variant of ANS called Assumed Natural Deviatoric Strain, or ANDES. Elements based on ANDES are described by Militello and Felippa (1991) and Felippa and Militello (1992). This unification work led naturally to a formulation of elasticity functionals containing free parameters. These were called parametrized variational principles, or PVPs in short. Setting the parameters to specific numerical values produced the classical functionals of elasticity such as Total Potential Energy, HellingerReissner and HuWashizu. For linear elasticity, 3 free parameters in a 3field functional with independently varied displacements, strains and stresses are sufficient to embed all classical functionals. Two recent survey articles with references to the original papers have been written by Felippa (1994,1996). One result from the PVP formulation is that, upon FEM discretization, free parameters appear at the element level. One thus naturally obtains families of elements. Setting the free parameters to numerical values produces specific elements. Although the PVP EulerLagrange differential equations are the same, the discrete solution produced by different elements are different. Thus an
144
Felippa, C.A.
obvious question arises: which free parameters produce the best elements? It turns out that there is no clear answer to the question because the best set of parameters depends on the element geometry. The PVP formulation led to an unexpected discovery. The configurations of elements constructed according to PVPs and the usual assumptions on displacements, stresses and strains were observed to obey certain algebraic rules. Such configurations could be parametrized directly without going through the source PVP. This observation led to a general formulation of finite elements as templates. FINITE ELEMENT TEMPLATES
A finite element template, or simply template, is an algebraic form that represents elementlevel stiffness equations, and which fulfills three conditions: (C) Consistency: the Individual Element Test (lET) form of the patch test, introduced by Bergan and Hansen (1975), is passed for any element geometry. (S) Stability: the element stiffness matrix satisfies correct rank and nonnegativity conditions. (P) Parametrization: the element stiffness equations contain free parameters. (I) Invariance: the element equations are observer invariant. In particular, they are independent of node numbering and choice of reference systems. Setting the free parameters to numeric values yields specific element instances.
Constructing Optimal Elements By making a template sufficiently general, all published finite elements for a specific configuration can be generated. This includes those derivable by orthodox assumptions and those that are not. Furthermore, an infinite number of new elements result. The same question posed above arises: Can one select the free parameters to produce an optimal element? The answer to this question is not yet known for general elements. An unresolved difficulty is the definition of unique optimality conditions at the element level. But even invoking reasonable criteria, a major technical difficulty arises: the actual construction of optimal elements poses formidable problems in symbolic matrix manipulation, because one has to carry along arbitrary geometries, materials and free parameters. Until recently those manipulations were beyond the scope of computer algebra systems (CAS) for all but the simplest elements. As personal computers and workstations gain in CPU speed and storage, it is gradually becoming possible to construct optimal twodimensional elements for plane stress and plate bending. Most threedimensional and curvedshell elements, however, still lie beyond the power of present systems. The mathematical theory behind templates is elaborate and will not be pursued in this article. Only a few basic results needed for the plate bending example presented in Section 4 are given in the next subsection. For further historical and technical details the reader is referred to a recent article by Felippa, Haugen and Militello (1995).
The Fundamental Decomposition The stiffness matrix derived through a template approach has the fundamental decomposition K  Kb(ffi) + Kh(~j)
(1)
Recent Developments in Basic Finite Element Technologies
145
Here Ko and Kh are the basic and higherorder stiffness matrices, respectively. The basic stiffness matrix Kb is constructed for consistency and mixability, whereas the higher order stiffness Kh is constructed for stability (meaning rank sufficiency and nonnegativity) and accuracy. The higher order stiffness Kh must be orthogonal to all rigidbody and constantstrain (curvature) modes. In general both matrices contain free parameters. The number of parameters O/i in the basic stiffness is typically small for simple elements. For example, in the 3node KPT elements considered below there is only one basic parameter, called c~. Its value must be the same for all elements in a mesh to ensure satisfaction of the IET. On the other hand, the number of higher order parameters flj can be in principle infinite if certain components of Kh can be represented as polynomial series of element geometrical invariants. In practice, however, such series are truncated, leading to a finite number of flj parameters. Although the flj may vary from element to element without impairing convergence, usually the same parameters are retained for all elements. A 3NODE KPT ELEMENT TEMPLATE The application of the template approach is rendered specific using a particular configuration: a 3cornernode fiat triangular element to model bending of Kirchhoff (thin) plates. The element has the conventional 3 degrees of freedom: one transverse displacement and 2 rotations at each corner. For brevity this will be referred to as a Kirchhoff Plate Triangle, or KPT, in the sequel.
Stiffness Decomposition For the KPT under study the configuration of the stiffness matrices in (1) can be shown in more detail. Assuming that the 3 × 3 momentcurvature plate constitutive matrix D is constant over the triangle, we have Kb  1 L D L T A '
T r Kh  A [BxT4DxBx4 + BxsDxBx5 + Bx6DxBx6 ]
(2)
3
Here A is the triangle area, L is the 9 × 3 force lumping matrix that transforms a constant internal moment field to node forces, Bxm are 3 × 9 matrices relating natural curvatures at triangle midpoints m  4, 5, 6 to node displacements, and DX is the plate constitutive matrix transformed to relate natural curvatures to natural moments. Parameter ot appears in L whereas parameters flj appear in Bxm. Expressions of these matrices are given in the Appendix.
The KPT136 Template A useful KPT template is ba', 0. The visible degrees of freedom of the element collected in u and the associated node forces collected in f are ur = [Uzl f r = [fzl
/~xl
Oxl
Oyl
~yl
blz2 Ox2
fz2
/~x2
Oy2
~y2
UZ3
fz3
Ox3
Oy3].
J~x3
/~y3 ] •
(5) (6)
The Cartesian components of the plate curvatures are tCxx, ICyy and 2Kxy = tCxy + tCyx, which are gathered in a 3vector ~. In the Kirchhoff model, curvatures and displacements are linked by
tCxx
=
02110 OX2 ,
Kyy

0 2W Oy 2 '
0211) 2gxy  2 ~ . OxOy
(7)
where w = w(x, y)  Uz is the plate transverse displacement. In the KPT elements considered here, however, the compatibility equations (7) must be understood in a weak sense because the assumed
154
Felippa, C.A.
curvature field is not usually integrable. The internal moment field is defined by the Cartesian components mxx, myy and mxy, which are placed in a 3vector m. Curvatures and moments are linked by the constitutive relation
[mxxI [OllO12O13][xx]
m =
myy mxy

D12 D13
D22 D23
D23 D33
I¢yy
(8)
 D~.
2tCxy
where D results from integration through the thickness in the usual way. Three dimensionless side direction coordinates zr21, zr32 are zq3 are defined as going from 0 to 1 by marching along sides 12, 23 and 31, respectively. The side coordinate 7(ji of a point not on a side is that of its projection on side ij. The second derivatives of w  Uz with respect to the dimensionless side directions will be called the natural curvatures and are denoted by Xji "" 021/3/07t'2These curvatures have dimensions of displacement. They are related to the Cartesian plate curvatures by the matrix relation
[XlI
X =
X32
O2W
[
X221 y221 x21Y21

~
X13

x~2
Y22 x32Y32
x23
Y23
x13Y13
02W
]
O2W ~X ~y
" T 1 tO,,
(9)
02w
2OxOy
the inverse of which is
02W
02W flUx ,,
02w 2 02w OxOy

1 [
Y23Y13
Y31Y21
Y12Y32
x23x13
x31x21
X12X32
Y23X31 d x32Y13 Y31Xl2nt" Xl3Y21
YlEX23 ] XElY32
02w
~
 TX.
02tO
(10) The transformation equations (9) and (10) are assumed to hold even if w(x, y) is only known in a weak sense
The Basic Stiffness Template Following Militello and Felippa (1991), the otparametrized basic stiffness is defined as the linear combination Kb  A1LDL T, L  (1  c~)L/d otLq (11) where L is a 9 x 3 forcelumping matrix that maps an internal constant moment field to node forces. Lt and Lq are called the linear and quadratic versions, respectively, of L. Matrix Lt was introduced by Bergan and Nygfird (1984) and Lq by Militello and Felippa (1991). Expressions for both matrices may be found in the latter paper.
The Higher Order Stiffness Template For an element of constant D, the higher order stiffness template is defined by
 A (B4TDxB4 + B/DxB5 + B6TDxB6)
(12)
Recent Developments in Basic Finite Element Technologies
155
where D x  T T D T is the plate constitutive relation expressed in terms of natural curvatures and moments, and Bxm are the natural curvaturedisplacement matrices evaluated at the midpoints m  4, 5, 6 opposing corners 3,1,2, respectively. These matrices are parametrized as follows. Define the geometric invariants X12X13 ~ Yl2Yl3
~'1 
X~l I y21
~'2  X23X21 ~ Y23Y21
'
X22 + y32
'
)~3 = X31X32 + Y31Y32
X23 q y23
(13)
"
which have a simple physical meaning as measures of triangle distortion (for an equilateral triangle, ~.1 = ~.2 = ~.3 = 1/2). In the following expressions, the/3derived coefficients )//and cri are selected so that the Bxm matrices are exactly orthogonal to all rigid body modes and constant curvature states. This is a requirement of the fundamental stiffness decomposition. Only the expression of Bx4 is given below. Those of Bx5 and Bx6 follow from appropriate cyclic permutations.
/~l = /~lO'qt]~ll~'3"~/~12~'l~/~13~'2 ' 1~2 ~~ /~20~/~21~'3qt/~22~1~fl23~'2 ' /~3 ~~/~30~/~31~3~~32~'1 "~'fl33~'2 ' /~4 ~ /~40 +/~41 ~3 ~ fl42/~'1 .ql_~43~,2 ' /~5  /~50 ql._/~51/~'3"~/~52~'1 ~/~53~'2 ' /~6  /~60 ~/~61 ~'3 ~ ~62~'1 ~/~63~'2 ' /~7 " /~70~'/~71~'3+/~72/~'1 +/~73~'2 ' /~8  /~80ql"f181~'aJC/~82/~'l Jf'/~83~2 ' /~9  /~90"~'/~91/~'3'~/~92~1 ~t/~93~'2' }/1 " fll + /~3' }/2 = f13' }/3 " f12 ~ f13' }/4 " /~4 @ f16' }/5 " 1~6' }/6 = f15 ~ f16'
}/7 = 137 +/~9, }/8 = 139, }/9 = 138 + 139, 0.1 = 2}/3, 0.2 = 2}/3  2Yl, 0"3 = 2}/4, 0"4 = 2}/6, 0"5 = 2}/7, 0.6 = 2}/9  2}/7, 0.7 = 2}/1, 0"8 = 2}/6  2}/4, 0"9 = 2}/9
Bx4 =
0"3 0"5
}/4Y31 + }/5Y23 }/7Y31 nt" }/8Y23 0"v }/lY31 ~ }/2Y23
I
}/4X13 ~ }/5X32 }/7X13 "nt }/8X32 }/1X13 q" }/2X32
0"4 0"9 0"1
f16Y31 '1 }/6Y23 f16X13 I }/6X32 /~9Y31 "t }/9Y23 f19X13 + }/9X32 /~3Y31 ~ }/3Y23 ~3X13 + }/3X32
0"8 f14Y31 d/~5Y23 /34X13+/35X32"] 0.6 /37Y31 +/38Ye3 ~7X13 ~ /~8X32 0.2
fllY31 + f12Y23
fllXl3 ~ /~2X32
J
(14)
This Page Intentionally Left Blank
STRUCTURAL OPTIMIZATION FOR PRACTICAL ENGINEERING: SOFTWARE DEVELOPMENT AND APPLICATIONS"
Yuanxian GU, Hongwu ZHANG, Zhan KANG, Zhenqun GUAN State Key Laboratory of Structural Analysis of Industrial Equipment Department of Engineering Mechanics. Dalian University of Technology, Dalian 116024, China
ABSTRACT The design optimization software JIFEX, with its former version MCADS, is developed with the applicationoriented concept. It is practically applicable to complex structures of general purposed engineering. The versatile structural modeling and simulation methods of JIFEX for the finite element analysis and design optimization are presented. The numerical approaches of structural sensitivity analysis for the problems of static stress, eigenvalue of vibration and buckling, dynamic frequency response, and transient dynamics and heat transfer are introduced. The software integration of the analysis, optimization, pre and postprocessing developed with advanced platform Windows 95/NT is implemented. The application examples of structural optimization has demonstrated the effectiveness of the optimization technology with applicationoriented software.
KEYWORDS Design Optimization, Sensitivity, Sot~ware
"Project supportedby the ScientificFundfor National Outstanding Youth of China (19525206)
157
158
Gu, Yuanxian et al.
INTRODUCTION With the development of modem computational technologies such as advanced computers and software techniques, finite element methods, computeraided design and engineering, and etc., the technology of design optimization is becoming more and more significant in the structural engineering. On the basis of theories and algorithms of the optimization, the numerical methods and software implementation are particularly important to the practical application of the structural optimization technology. In the practical applications of the real life engineering, the numerical methods and software implementations of design optimization have to be composed of the following components: The numerical modeling and simulation of complex structures in the analysis and the optimization. The structural analysis with general purposed capabilities. The efficient sensitivity analysis, i.e. the sensitivity derivatives calculation for structural behavior functions. The efficient and robust solution algorithms of optimization problems. The reliable and applicable software with largescale numerical computing capability and user friendly operations. In the past three decades, there was much research on the numerical methods of structural optimization and some programs developed. However, most of these kinds existed programs are research oriented and insufficient in the capabilities such as the types of elements, loads, and variables, functions of analysis, constraints as well as sensitivity calculations, and computational scale. This has limited the applications of the structural design optimization technology in practical engineering. The problem is that the linkage between the theoretical research and the software development, and the requirement of sottware being application oriented. Due to the complexity of practical engineering, the abovementioned facilities are necessary for structural optimization softwares. This paper presents the development and applications of some numerical methods of structural design optimization implemented in the software system JIFEX t~l and its former version MCADS t2]. The development of JIFEX software is practical application oriented and based on the advanced computer techniques. First of all, the JIFEX and MCADS were developed on the basis of general purpose software for the finite element analysis of complex structures. Its general purpose analysis for the static, dynamic, and buckling problems is supported by the versatile modeling capability to complex structures with various elements, loads, and boundary conditions. On the basis of these structural analysis and modeling capabilities, the numerical methods of sensitivity analysis have been studied for the constraints of static stress and deformation, eigenvalues of vibration frequencies and
Structural Optimization for Practical Engineering
159
buckling loads, dynamic responses in the frequency domain and time domain, as well as the transient heat transfer in structures. This way, the design optimization model of JIFEX is implemented for various elements, loads, variables, constraints, and objectives. The JIFEX system is developed on the platforms of Pentium PC and MS Windows 95/NT, with the 32bit C/C++ programming tool. The preprocessing of the finite element modeling is developed on the CAD package AutoCAD, and implemented the interactive geometric modeling and automatically generation of finite element model data. The computing visualization under the Windows 95/NT environment is available. The advanced technologies of computer sottware and numerical computing make JIFEX to be new generation software of structural analysis and design optimization. STRUCTURAL MODELING OF ANALYSIS AND OPTIMIZATION Firstly, JIFEX software is capable of simulating complex structures by means of various kinds of elements, loads, and boundary conditions. Particularly, it simulates boundary support and component connection conditions with displacement constraint on each freedom of node in arbitrary local coordinate system. And the displacement constraints include status of fixed, preassigned displacement value, and particularly, multilevel masterslave relation. The socalled masterslave relation means that the displacement status of a freedom of slave nodes can be controlled by (i.e. with the same value of) one master node and this control relation can be defined level by level. These methocts make the simulation of complex structures convenient and easy. The finite element analysis of JIFEX generally applicable to the design optimization of complex structures includes the following fundamental requirements: • Structural static strength analysis, • Free vibration frequency analysis, • Dynamic response analysis both in frequency and time domains, • Global buckling stability analysis, • Transient heat transfer and thermal stress analysis, In accordance to the versatile finite element modeling, the design optimization model of JIFEX has developed following three types of design variables: • Size design variable. It includes the geometric sizes of cross section of various elements, such as bar, beam, membrane, plate, and shell elements. • Shape design variable. It includes the coordinates of special nodes and geometric parameters of boundary shape interpolations. Then, the boundary shape of continuum structures, the configuration of frame structures, the locations of stiffeners and linkages can be optimized. The shape optimization modeling is cooperated with the parameterized design concept of the CAD. • Composites design variable. It includes the design parameters of composite laminate and honeycomb sandwich plates, such as ply orientation, layer thickness, core height, and material parameters of special composite components.
160
Gu, Yuanxian et al.
The design constraint functions of JIFEX95 cover the following structural behaviors: structural weight, node displacement, stresses, vibration frequencies, dynamic displacement and stress responses, and buckling loads. The design objective function can be selected from any one constraint function or combined with several constraint functions. This way, the objective of design optimization can be = Minimizing structural weight, • Reducing structural stress and/or deformation, • Improving structural stiffness, • Increasing structural fundamental frequency or adjusting the distribution of a group of vibration frequencies, • Minimizing structural dynamic responses, = Increasing structural buckling loads, • Multiple objectives optimization. The basic feature of the modeling and simulation facilities of JIFEX is its general applicability to complex structures. This applicability is presented from the finite element model, analysis functions, design variables, design constraints, and optimization objectives. JIFEX is also capable of largescale engineering computing. The finite element model with 100000 nodes can be computed with Pentium PC. Some analysis applications of practical structures with this scale have been completed.
STRUCTURAL SENSITIVITY ANALYSIS The optimization solution algorithms of JIFEX are sequential linear programming and sequential quadratic programming. The solution efficiency and convergence stability of these algorithms are based on the sensitivity analysis on the constraint and objective functions. The difficulty of sensitivity analysis for the general purpose optimization software such as JIFEX is that it must be suitable to various elements, variables, and structural responses. The numerical methods of sensitivity analysis have been studied for (1) static stress and deformation, (2) eigenvalue problems of free vibration and buckling stability, (3) frequency response problem, and (4) transient problems of dynamic response in time domain and heat transfer in structures. Static Strength Problem
The sensitivity derivatives of displacement vector U and stresses o are calculated with Eq.(2) derived from Eq.(1), the finite element formulation. The K and P are structural stiffness matrix and load vector respectively, S is element stress matrix, subscript e denotes element. KU=P,
cr=SeU ,
x u ' = p '  xz , d  s'u, + s y ;
(1)
(2)
Structural Optimization for Practical Engineering
161
Eigenvalue Problems The eigenvalue problem can be expressed as below K~b 3M~ = 0
(3)
Two kinds of eigenvalue problems have been considered: (a) The freevibration, M is mass matrix, and ;t is the square of vibration frequency; (b) The buckling problem, M is geometric stiffness matrix, 2 is critical load factor. The derivatives of the eigenvalue are obtained by Eq.(4). 2.' ~br (K'AM')~b qbrMq~ = 1
(4)
The following Rayleigh quotient can be used to improve the approximation accuracy of the eigenvalues.
~j(x): uj(x)/~(x) O , ( x )  v,(Xo)+ ~=1
(5)
ev, (x,  Xo,), r,.(x) .. . er,. :,,, = r,.(Xo) + 2=~r.,
~=l v' x,
 Xo,)
where Uj and Tj are modal strain energy and modal kinetic energy of the jth eigenvector. For the buckling problem, Eq.(5) is approximate since M depends on prebuckling stress field. The accurate derivatives are
,~' = #~ ( K '  aM')~ + ~A~K'U~
(6)
where Us is prebuckling displacement, and the adjoint field A is obtained by K A : {D(#rMqk)} r cTJB
(7)
Dynamic Frequency Response Problem
The dynamic frequency response problem is formulated as below Mii + Cft + Ku = f
(8)
The C is damping matrix. The harmonic exciting is expressed asf=foCos(0t)orf=foSin(0t). Two kind of loads fo are dealt with: (a) the nodeforce load exciting is fo = P, (b) the accelerate movement exciting is fo = MAUo. u o is the rigid displacement vector with unit
162
Gu, Yuanxian et al.
movement, A is accelerate scope. The displacement response of Eq.(8) is: u  s sin(0 t) + c cos(0 t)
(9)
With the modal superposition method, the derivatives of eigenvectors are required. It is difficult in the case of repeated eigenvalues. Following subspace iteration method and direct eigenvector space method are implemented in JIFEX. S u b s p a c e iteration m e t h o d
Starting from trial vectors qj(0), compute eigenvector derivatives by following subspace iteration procedures with k=0,1,2,. .....
(a) Inverseiteration: ( K + o M ) V (k+°  f f + Mq/(k)
(10)
w
F  M'~bp  ( X ' + ~ ' ) W (ll) W
~ppl,
p = Ap + ¢rI
Ap and ~p are p eigenpairs, o is shift value. (b) Projecting into trial vector space for eigenvector derivatives: /(D (k+°  A~/D(k+')Ap = A~A].
+ A~'(k+0pAp

~,(k+,)p
(12)
~,(k+l) _ W r K , W + W r K V ~k+° + (W rKV(k+°) r l~l'(k+') = W r M ~ V + WrMV (k+' + (WrMV(km) r  Wr KW,
(13)
I~t  W r M W
(d) Updating the trial vectors and convergence checking: T(k+~) _ V(k+Op + W D (k+~)
(14)
The Nelson's method is used to solve Eq.(12) for repeated eigenvalues. E i g e n v e c t o r space m e t h o d
Let C be independent to variable, excitingf=foCos(0t). Differentiating Eq.(8) gives Mii' + Cit' + K u '  f s sin(0 f s  M ' s 0 2  K's,
t) + .[c cos(0 t)
fc  M'02c  K'c + fo
(15)
Structural Optimization for Practical Engineering
163
Solving Eq.(15) in eigenvector space twice with respect exciting fsSin(0t) and fcCos(Ot) respectively, we have s
Ss sin(0 t) + c s cos(0 t)
(16)
U'c = sc sin(0 t) + c c cos(0 t)
Then, the derivatives of u can be obtained {u}' : {u~}' + {u~}' : ({S~} + {S~})Sin0 t + ({C~ } + {Co })CosO t
(17)
Both of the methods have advantages of easy to be implemented on the basis of dynamic analysis, computational efficiency and accuracy, but without difficulty of eigenvector sensitivity calculation with repeated eigenvalues Transient Problems
The motion equations and criticalpoint constraint of transient problems can be written in a general firstorder form Air = f (u,a,t), g ( u , a )  f ' p ( u , a , t ) d t < O,
(18)
u(a,O) = u o p ( u , a , t )  g(u,a,t)fi(ttm~)
(19)
where u is response vector, a is design variable, g is a general form constraint function, t denotes time. The sensitivity analysis with the direct method and the adjoint method are implememed with Eq.(20) and Eq.(21) respectively. A dit _ j _ _du dot da
dA i, + f da Oa'
AT2 + (jT + AT ))t, " 0 ,
du(a,O____.__~)= 0, J, j. = ~ da ' &tj
AT2(t,,, ) = .c3g(tm, ) \ &l
(20) (21)
For the transient heat transfer problem, the finite element formulation is MT+KT=Q=P+~
(22)
where M is the heat capacity matrix; T is the vector of node temperature; Pb is reduced by the given temperature on the boundary. With the implicit 0difference method, the temperature T"~' of time t,+~ is computed from T' with
v (a,,~l + K)T "+'  (a,,M  b K ) T " +0 a,=l/(OAt,),
b=(1O)/O,
00
d) calculate U~+,~t.:,
Ut+,~t,:
4) Output the computational results if required, then turn to 3) for the next step, else stop at the last step. It is noteworthy to point out that the parallel efficiency of this parallel algorithm on the substructure level is dependent on the structure decomposition. Thus, when decomposing and dispensing the substructures, we should pay attention to two points: 1) Balance load among different processors as possible as we can get. Otherwise, the parallel efficiency will be decreased because of the heaviest loading processor following the Amdahl principle. 2) Interfaces between the substructures (also the number of the correlative elements of each substructure) should be minimized to decrease the communication between processors and the computation repeatedly.
PARALLELIZATION OF NEWMARK INTEGRATION METHOD For Eqn. 1, Newmark integration method employs the following timemarching scheme: (11)
U~+,~tU~+ U~At + I ( 1  a ) ~ )
' + aU,+z/lAt2
(12)
By satisfying the Eqn. 1 at time t + At, we obtain
Rg),+~,  k,+,,, Where
(13)
Parallel Integration Algorithms for Dynamic Analysis of Structures System
[(  M + 6 AtC + teat 2K
173
(14)
(15) K
U, +UtAt+(2
Computing time for the Newmark integration is mostly spent on the assembling of global matrices, evaluating effective load vector R, formation, factorization and back/forward substitutions of effective stiffness matrix /f through the each time step. Assembling
Since the effective stiffness matrix /f is generally of large scale in the FEM analysis of engineering structures, the assembling, factorization and back/forward substitution, which have to be segmented into blocks even in the modern computers, and the incore outcore swapping techniques have to be employed. A routine in the available PFEM provides the function of automatic segmentation of global matrices, which can be easily used for the parallelization of assembling. Suppose that m processors are available and t( has been separated into il blocks. Then the kth processor cyclically processes block k, block re+k, etc., and store them in the local hard disk. Computing effective load vector
After loading element stiffness matrices, the correlative nodes and the load vectors distributed on slave processors are updated simultaneously and are further transferred among the processors to fulfill the final load vector updating. Triangular f actorization
The global stiffness matrix can be factored: R  Lz)/5
(~ 6)
where L is the lower triangular matrix and D is the diagonal matrix. The global matrices are stored in onedimensional array by employing the variable bandwidth technique. The largest bandwidth is denoted as d. Denote m i as the row number of the first nonzero entry of column j, for which Smj,j  Kmj,j
(17)
For other elements of column j i1
 max(mj , ran) r=mn
174
Cheng, J.G. et al.
where L 0  S , j / d , ( j  2,3,.., d ; i 
(19)
1,2,.., n)
i1
di,  K ,  ~ L,.,r d r,.Lr,r
(i . . . .1,2,
, n)
(20)
rm!
Block 2 ~
......
Figure 1: Blocking diagram of the variable bandwidth stiffness matrix The block partition of A" is illustrated in Figure 1. Triangular blocks and rectangular blocks are denoted as T1, T2 ..... etc. and R1, R2, ..., etc., respectively. The triangular factorization of T1 is the same as the original. The factorizations of T2 and R1 are affected by the factorization of T1. Generally the process of the following blocks is affected by its predecessors, and the blocks had to be swapped between the incore memory and the secondary memory. Denote the triangular block size as n, and partition it into b=[n/m] subblocks. Thus processor Cn processes column Cn of each block: (1) decompose the first subblock in every processor (2) processor 1 decomposes column 1 of subblock 2 and transfers the decomposed element into the other processors. (3) processor 2 decomposes column 2 of subblock 2 and transfers the decomposed element into the follower processors. (4) processors 3, 4 ..... m carry out the same processing as processors 1 and 2 (5) processor m transfers column m of block 2 into the former processors (6) processor (ml) transfers column m of block 2 into the tbrmer processors (7) processors (m2), (m3) ..... 2 complete the message passing as processors m and (ml) Repeat the above procedures till current block is factorized.
Parallel Integration Algorithms for Dynamic Analysis of Structures System
175
Each rectangular block should be partitioned into jk=[n/2m] subblocks. Every subblock is decomposed simultaneously as follows: (1) processor 1 decomposes subblocks 1 and jk (2) processor 2 decomposes subblock 2 block jk1 (3) processors 3, 4 ..... m complete the parallel work as processors 1 and 2 The results are transferred from each other after decomposing the last rectangular block in front of the triangular block.
Back substitution Considering Eqn. 1, Eqn. 13 can be solved by the following two steps, i.e.
Ly  k DLT ~f t+At : Y
(21)
(22)
The loading manipulations from outcore memory are intensive for large problems, so the parallelization of the loading manipulation is as important as that of algorithms. Here we take Eqn.21 to illustrate the procedure: (1) processor 1 loads block 1, updates the corresponding part of y , and transfers results to the following processors. (2) processors 2, 3, 4 . . . . , m load blocks 2, 3, 4 . . . . . m respectively and update the corresponding part of y with the triangular part of correlative blocks. Then these processors receive the results transferred from the former processors, update y, transfer the updated parts to other processors, and so on. For blocks m+ 1, m+2, ..., the same solution procedures can be used.
NUMERICAL EXAMPLES
Examples of the parallel central difference integration algorithm The computations are performed on three pyramidshaped space DOFs 128, 288 and 512. The implementations are carried out in the speed hubclustered PCs. The pre and postprocessing are finished FEM package to which the parallel algorithms are incorporated. analysis are fully coincident with those of the serial analysis.
flames, respectively, of PVM platform on highwith an available serial The results of parallel
The effectiveness of a parallel algorithm can be assessed by the speedup ratio Sp(n) and the parallel efficiency Ep(n), where Sp(n)=ts(n)/tp(n), Ep(n)=Sp(n)/n. ts(n) and b(n) are, respectively, the computer time used on the single processor with the best serial program and that on the n processor system with the parallel program.
176
Cheng, J.G. et al.
The achieved parallel efficiencies are shown in Table 1. These parallel efficiencies increase with the number of time steps and the size of the problems.
TABLE 1 PARALLEL EFFICIENCY ON PVM Examples
1
2
3 i
Number of processors
Time steps of load
2
3
4
5
500
0.92
0.86
0.76
0.72
800
0.93
0.88
0.77
0.69
1000
0.92
0.84
0.77
0.71
500
0.96
0.91
0.83
0.77
800
0.94
0.89
0.84
0.76
1000
0.95
0.89
0.82
0.75
500
0.95
0.90
0.84
0.76
800
0.94
0.90
0.82
0.75
1000
0.93
0.90
0.82
0.76
'
Examples of parallel Newmark algorithm Computations are carried out on a plane stress problem and a 3D bending plate problem, shown in Figure 2. The material and geometrical parameters are as follows: E = 2.1e + 11N / m', v  0.3, p  7.84e + 3 kg / m 3, l  100m, h  20m , b  0.1m, co19.4224,a: =  1 0 0 0 N , a , 
100
100
100
N,a0.25
80.50,02.00,At=0.04.
Speedup Sp and parallel efficiency Ep under different DOFs n of the structure, bandwidth d, number of blocks il, and time steps are depicted in Table 2. F(t)  a a sin co t
/ / / / /
F(t)  a 2 sin co t
Epv , ~
l
Figure 2 Structures and parameters
F ( t )   a 2 s i n cot
Parallel Integration Algorithms for Dynamic Analysis of Structures System
177
These results show that the parallel e~ciency increases with the DOFs and the bandwidth of the problem, and with the time steps to be solved. For the available network system and the present examples, time spent on parallel computation is relatively smaller than that on communication, so the parallel efficiency of triangular factorization is lower comparatively. Since the stiffness matrix in Example 1 is partitioned into 6 blocks, the parallel efficiency is on the low side because of the load unbalance when the number of the processors is 4 and 5. The super linear speedups are observed in Example 2 due to the relatively larger buffers for swapping when the number of processors increased. TABLE 2 PARALLEL EFFICIENCIES 1N THE EXAMPLES Example 1
Example 2
N
3262
7579
Number of
D
318
726
Processors
I1
6
32
2
3
4
5
Steps
20
50
100
20
50
100
Sp
1.67
1.83
1.88
1 . 8 1 1.89
1.96
Ep
0.84
0.92
0.94
0.91
0.95
0.98
Sp
2.30
2.62
2.65
2.61
2.75
2.91
Ep
0.77
0.87
0.88
0.87
0.92
0.97
Sp
2.75
2.69
2.82
3.34
3.52
3.78
Ep
0.69
0.68
0.71
0.84
0.88
0.95
Sp
2.79
2.75
2.88
3.56
4.04
4.49
Ep
0.56
0.55
0.58
0.71
0.81
0.90
CONCLUSIONS In this paper, the central difference method and the Newmark method are parallelized on a clustered network system and are also incorporated into a generalpurpose finite element package  PFEM, which results in a serialparallel mixed FEM package. Satisfactory parallel efficiencies are obtained for the examples presented. In both algorithms, the parallel efficiencies are found to increase with the size of the problems. In the Newmark algorithm, super linear speedups have also been observed, the reason for which can be ascribed to the increase of the swapping buffers with the number of processors and the swapping intensive nature of the devised algorithm. Due to limitations of the hardware and software, the scale of the examples presented is rather small; it will be increased in further work.
178
Cheng, J.G. et al.
REFERENCES
Cheng J.G. and Li M.R. (1995). A serialparallel mixed finite element analysis soitware PFEM. Computer and Mechanics Research. Chengdu: Chengdu Science and Technology University Press, China. Cheng J.G., Yao Z.H., Li M.R. and Huang W.B. (1996). Parallel algorithms for explicit integration of dynamics analysis of structures and their implementations. Journal of Tsinghua University, Science and Technology, 36:10, 8085. Chiang K.N. and Fulton R.E.(1990).Structural dynamics methods for concurrent processing computers. Computer & Structure 36:6, 10311037. China Architecture Academy, Zhejiang University (1991), JGJ 791. People's Republic of China jobs normalframe structures design and described programs for implementation. China Architecture Industry Press, China. Hajjar J.F. and Abel J.F.(1989). Parallel processing of central difference transient analysis for threedimensional nonlinear framed structures. Communications in Applied Numerical Methods 5:1, 3946. Noor A.K.(1979). Finite element dynamics analysis on CDC Star 100 computer. Computer & Structure 10:1, 719. Ou R. and Fulton R.E.(1986). An investigation of parallel numerical integration methods for nonlinear dynamics. Communications in App#edNumerical Methods 30:1/2, 403409. Wang X.C., Shao M.(1988). The fundamental principle of the finite element method and numerical methods, Tsinghua University Press.
M A X I M U M E N T R O P Y PRINCIPLE A N D TOPOLOGICAL OFITMIZATION OF T R U S S S T R U C T U R E S
B. Y. Duan, Y. Z. Zhao and H. Liu School of Mechanical and Electronic Eng. Xidian University, X f an 710071, China
ABSTRACT This paper presents an entropy  based topological optLrnization method for truss structures. A group of variables meaning the density distribution of energy is defined first, so that a bridge between energy distribution and topological optimization is built. After the energy distribution is found, optimal cross  sectional areas are f o u n d . Then the above two steps are executed in turn till the convergence is reached. Next, the stiffness and strength reliability corresponding to the final topology are computed with reasonable results.
KEYWORDS maximum entropy principle, topology optimization, structural reliability, linear programming
INTRODUCTION Maximum entropy principle(MEP) [1] has been applied in most information systems. Since engineering analysis and optimum design can be viewed as an information system, ~ has been applied in this area too. Simoes and Templeman E21applied IVIFP in a synthetic problem of pretensioned steel net; Li and Templeman [sl[s] utilized MEP to solve structural optimum design; Erlander [6] used MEP in problems of distribution and assignment of traffic; Tiku and Templeman [sl applied MEP to the problem of analysis and optimum design of underground water net in a city. Recently, Abraham I. Beltzer [141 applied ~ to finite element with some interesting
179
180
Duan, B.Y. et al.
results. Study on topological optimization has continued several decades since the first paper Es] was published in 1964. So far, the two main methods are Ground Structural Approach(GSA) and Homogeneous Method(HM) [1°][11]. Both methods view topological optimization from a viewpoint of either force transmitted or material distribution. Whether it is the force transmitted or material distribution, both can be classified as a problem of information treatment in engineering. That is to say, topological optimization is also information treatment which could be considered from an information entropy concept. It can be anticipated that a new area of study might be opened if topological optimization can be studied from the angle of information treatment by MEP. This paper attempts to make such a study. MATHEMATICAL MODEL OF ENTROPY  BASED METHOD FOR TOPOLOGICAL OPTIMIZATION OF TRUSS STRUCTURES I n t r o d u c t i o n o f D e s i g n Variables a n d Objective F u n c t i o n s
Considering a truss from the viewpoint of energy distribution reveals that different topological forms correspond to different energy distribution. Suppose the total structural energy is Q0 and energy of the/th bar is q~, both Q0 and q~ will satisfy Q0 =
~]q~
(1)
i1
in which n is the total number of elements. Now a set of variables {~}~"   (~vl, J~2, " " , ;~
q' 
Q0
~)T
(i= '
is defined as 1
2 '
... '
n)
(2)
'
Variable & is of obvious physical sense, its value stands for the percent of energy stored in the/th element. If & equals zero, then the/th bar can be eliminated from the ground structure. As a resuit, all those elements with variable & (i  1 , 2 , . . . , n) being non  zero construct a topological form under consideration of a certain sense. It is obvious that variables & (i  1 , 2 , . . . ,n) satisfy equations (3) i1
and
&~0,
g iE
F
(4)
where F is a set consisting of all connectable elements for the given fixed points. Because variables & (i = 1 , 2 , . . . , n) satisfy non  negative and normality conditions required by entropy, the corresponding structural entropy can be described in the form of
Maximum Entropy Principle and Topological Optimization of Truss Structures
S =
O~~& • ln&
181
(5)
i=1
where O is a large positive constant and S is structural entropy. In adch'tion, an e n e r g y  related function can be described in
f ~,.
IF, I
(6)
i1
in which Ti __
"]l;,['_.~.q,...~
~q~~ is the permissible stress of the/th bar ~ L~ is the length of the/th bar 2E E is the modulus of material.
It can be noted from Eqn. 6 that volume of f is directly proportional to structural weight for the case of fully stressed design. In the general case ( n o n  fully stressed design), this will not be valid. However, what we are interested in is if the ith bar exists and what is the concrete response in f. Actually, the situation whether the ith bar exists or not can be mathematically shown by if [F,[ ~ 0. According to the sense of topological optimization, [F~ I ~ 0 means that the ith bar will tend to be eliminated from ground structure. That is to say, the/th barfs weight will tend to be eliminated from the whole structural weight, which means the decrement of structural weight. Seen from this point of view, both problems of minimizing f and minimizing structural weight are equivalent. Up to this point, two parameters have been defined. One is the structural entropy expressed in Eqn. 5 and the other is the function shown in Eqn. 6. It can be seen from the following discussion that both functions can be expressed as the functions of variables {Z }r = (Zl, Zz, . . . , £,)r
Description of Behaviour
Constraints
Substituting energy definition of the ith element into physical equation yields 1
q, =  ~ g , .
A,~
(7)
where K~ = EA~/L~ is the ith barfs stiffness. A,, A and F~ are cross  sectional area, elongation and member force respectively. Solving A from Eqn. 7 and considering Eqn. 2 leads to A=+_
~2O0&K,
,
(i=
1, 2, . . . , n)
(8)
Similarly, the member force F~ can also be written as
F, = + J2K, o0
,
= 1,2,...,,)
(9)
Substituting Eqn. 9 into equilibrium equation [ N J { F } = {P} results in
(~0)
Duan, B.Y. et al.
182
+ ~/2K1Qo& rN]
.

{P}
(11)
I ~ / 2 K . O o & in which I N ] is a projection  geometry matrix; {P } is a loading vector to which structure is subjected. From geometry equation {A} = [N]~{6}
(12)
and considering boundary condition, we have {6} = ED]{A}
(13)
where {A} is a vector of elongations; {6} is a vector of displacements ~ matrix [D] 
EEN]EN]~']I[N] . ff the allowable values of m displacement constraints are given as {6}  (a~, a2, ..., a.)~ , displacement constraints can be, by substituting Eqn. 8 into Eqn. 13, described as, + ~ 2O0& KI
ED]
:
{a}
(14)
~2O0&
+
K,
As for stress constraint, it can be written in the form of
~/2K, Oo& ~ A,[a,],
(i
1, 2, " " , n)
(15)
where [~r~] is the permissible stress of the ith element. Up to this point, all constraints such as equilibrium equation, displacement and stress have been described as the functions of variables {£}  (&, ~2, " " , &)~".
Mathematical Model of Topological Optimization Synthesizing Eqns. 3 , 6 , 1 1 , 1 4 and 15, and considering minimum weight and maximum entropy design simultaneously, yields the following mathematical model. (P I )
find variables
{~}  (~1, &, " " , £,)~ in space
minimize
f  ~ T, . F, i=1
maximize
~ ~, . ~/2K, Oo&
(16)
i=1
S =  O ~ & • In& i=1
(5)
Maximum Entropy Principle and Topological Optimization of Truss Structures
183
+ J2KIQoX1
EN]
subject to
.
 {P }
(11 )
! ~/2K,,Qo~
:J:
~2Q0~1 K~ i
~_ {~}
EDJ
(14)
+__ ~/2OO~K~
~/2K~qo& ~_ A, Ea~], &
=
(i = 1, 2, . . . , n)
1
(15) (3)
i=l
>_ 0
V i E r
,
(4)
Actually non  negative condition Eqn. 4 can be satisfied automatically due to the logarithmic form of entropy function. In order to simplify problem (P I ) , to eliminate the absolute symbol and to have a standard LP model, let = 2Q0
'
(i=1
'
2
'
...
'
n)
=1.2.
(17) .
and adding slacking variables {Y}  (yl, y2, " " , y~)~ and {Z}  (zl, z2, . . . , z.)~to the left hand side of Eqns. 14 and 15, respectively, and introducing two weight coefficients ~1 and ~)2, yields the following problem (PII)
find variables
({~}~, {~}~)~ = (~,, ~2, ..., ~ , ~,, ~2, ..., h . ) ' i n space E 2"
n ~ ~
¢ =
subject to
~1~,(~. + ~) + i=1
I[DK] L [I..]
IEDK'] [I..]
~ O ~i= 1 ~
  ~>" " In 200
(~ 200   ~>"
[I,,,,,,] [0,~]0 { ~ ~} EO..] [L.]J z}J L{~>J
(19)
(20)
{a} , {a} , {Y} , {Z} ~ 0
(21)
2
(22)
i=1
(~. _ ~.)2 = 200
184
Duan, B.Y. et al.
0 (23)
°o °

1
[OK] = Eo]
(24)
9 Oq
i [~,]
=
A![~,]
(25)
[ I ~ ] is a unit matrix. [0~.], [0,,] and [0..] are zero matrixes. In problem (P 11 ) , the expressions except for entropy and normality condition Eqn. 22 are all linear functions, provided cross sectional areas are given. Considering the characteristics of problem (P 1I ) , SLP is utilized to solve it without considering the normality condition first. METHOD OF FINDING OUT ENERGY DISTRIBUTION The normality condition of Eqn. 22 is dealt with in the following method. For the known variables ({a (k) }~, {a (k) }~)~ after the kth iteration, let 2Q0
,
( i   1, 2, ..., n)
(27)
Then parameters {A(k)} are normalized by
~where
~= ~
,
(ii,
2, ..., n)
(28)
)~(')
i=1
Next, update variables
a~k) = ~(k) _ ~(k)
,
&(') =_+ ~/2Q0~(~)
(i  1, 2, " " , n) as,
,
(i
1, 2, ..., ~)
symbol ± of ~(k) should be the same as that of ~(k). Finally, ~k) is considered as the present design point ; go to the next iteration again. METHOD OF FINDING OUT CROSS  SECTIONAL AREAS
(29)
Maximum Entropy Principle and Topological Optimization of Truss Structures
185
The previous discussion is merely concerned with seeking {,~}  (,~1, Z2,  " , ,~)~ under the known area {A }  (A1, A2, ".., A,) r. As a matter of fact, topological form will be influenced by {A } . Meanwhile, what is necessary in engineering is to know the cross  s e c t i o n a l area. In order to find optimum topology and cross  sectional areas as well, both parameters of {,%} and {A } should be viewed as variables. ff the work of finding {,%} under the given {A } is called step ( I ) , process of seeking {A } is called step ( lI ). The final result of {Z } and {A} can be obtained by carrying out two steps in turn.
Problem of seeking cross  s e c t i o n a l area under the given energy distribution can be described as follows. ( p 11I ) find variables
{¢}  (¢1, ¢2, " " ,
¢.)~in space
pL
T
minimize
W
maximize
S =  O ~
,=~
(30) ,8,g, • lnfl,¢,
(31 )
i=1
subject to
"~, 1;
are the piecewise linear forces of elastic components and damping, respectively. To provide the theoretical background for the design, it is essential to study the nonlinear primary resonance of the system when maxly] >1. If the parameters (0,4"1,/~ and f are small, the primary resonance can be approximated, through the use of average approach, as
(5)
y(r) = a(r) cos[2r + (p(r)], where a(r) and (,o(r) yield
a  q ( a )  f~coscp, 22
ip = 1 ~2 + p(a) 22
2a
+
f 22a
(6)
s i n co,
I p(a)  pa ~ (2(,oo  sin2cp o ), (7) [q(a) = a[4:07rg + ~:l (2(f10  s i n 2c,o0) ],
(Po = arcc°s(1), a
a > 1.
From Eqn.6, one has the relationship between the amplitude of the steady state resonance and the excitation frequency [ p ( a ) +  ~a ( 1  22 )] 2 + ,~,2q2( a )  ( ~ f) 2 = 0 .
(8)
To classify the types of the primary resonance, one can focus on the case of a >> 1 and let z = 1 / a 0 , which is the natural consequence of (0 > 0 and (1 > 0 in practice. Figure 2 shows a typical transition set of these two unfolding parameters. The transition set divides the right half plane of (a I , a2) into 4 regions as shown in Figure 2. In each region, the amplitudefrequency curve of the primary resonance looks the same qualitatively. This figure, thus, enables one to choose an appropriate combination of unfolding parameters so that the vibration isolation system possesses the desired qualitative behavior in primary resonance. It should be emphasized at the end of this section that though the above analysis is made for a piecewise linear system on the assumption of weak nonlinearity, i.e., the parameters ~ and (l are small, the numerical simulations in Wang and Hu (1997) showed that the results were valid even when /~ and (1 w e r e not small.
DESIGN
APPROACH
The basic idea of present approach is to design a linear vibration isolation system with very small damping first, and then an elastic stop with large damping. In the working frequency range, the slightly damped vibration isolation system has required vibration transmissibility. Once the vibration isolation system undergoes the primary resonance somehow, both the elastic stop and the large damping reduce the vibration amplitude and remove the jumping phenomenon that may occur for a harmonically forced nonlinear oscillator. 3.1 Design of primary system The vibration isolation system without any stop is referred to as the primary system hereinafter for brevity. The vibration transmissibility of the primary system yields T= I.
1+(2~:oA)2 (1  2 2 )2 _+. (2~:o~)2
'
(15)
where only two dimensionless parameters 2 and (0 are to be designed. For a linear vibration isolation system in traditional sense, the vibration transmissibility in resonance can only be attenuated by increasing the damping ratio (0. For the vibration isolation system with an elastic stop, however, the task of attenuating the vibration transmissibility in the case of resonance can be left to the damped stop. Hence, a very small damping ratio 4"0 can be chosen in the design of primary system in order to avoid the system impacting the stop when the system starts running, see Hu (1996). In the case of 4"0 < 0.1 and 2 > 2, the stiffness of the main elastic component can be determined by using the following approximation of Eqn. 15
A~_41+1/T .
(16)
236
Hu, H.Y. and Wang, F.X.
3.2 Design of Damped Elastic Stop 3.2.1 Preliminary Design Given the dimensionless excitation frequency 2, the damping ratios (0 and (~, the system parameters to be designed are only /~ and f . These two parameters appear in the expressions of unfolding parameters a 1 and a 2, and hence, follow the selection of the two unfolding parameters. In principle, any parameter combination of a 1 and a 2 in region IV in Figure 2 makes sure that the frequencyamplitude curve of primary resonance does not have jumping. Thus, an arbitrary combination of (a 1,c~2) in region IV can be chosen to determine the corresponding parameters /.t and f , provided that the vibration transmissibility is acceptable. A great number of numerical simulations showed that the stiffness ratio y should not be too large. As shown in Figure 3, the function of stop is very obvious in the lower frequency range. The response amplitude goes down very rapidly in the beginning of the increase of /1, and then changes not verymuch later until very complicated dynamics happens. For the sinusoidal excitation of high frequency, the response amplitude has a peak as shown in Figure 4. It is smaller than the initial value only when ¢z is very large. As a result, an excessive stiffness ratio ~ is harmful. 1.6
3.4 2.85
1.5
Ymax
Y max
1.4.
1.3
2.3
~
1.75 j
0
. . . . . . . . . . . . . . . .
2
4
8
8
10
Figure 3" Maximal displacement versus stiffness ratio when )~ = 0.5
1.2
0
.......... 2
4
8
8
10
Figure 4' Maximal vibration amplitude versus stiffness ratio when )~ = 1.4
Moreover, the dimensionless excitation amplitude f
defined in Eqn.1
is inversely
proportional to the clearance 6 when the excitation amplitude F is fixed. So, the clearance can be determined from c~2. If the clearance is too large, the stop can not be in function, tf too small, the vibration may become nonlinear and then undergoes a subharmonic resonance in working frequency range. In summary, the stiffness of the stop should not be very large and the clearance should be appropriate. So, it is necessary to optimize these two parameters, or namely two unfolding parameters in region IV in Figure 2.
Computer Aided Design for Vibration Isolation Systems with Damped Elastic Stops
237
3.2.2 Optimization for Parameters Now f , the dimensionless excitation amplitude, is taken as the design variable to look for the minimal stiffness ratio of an elastic stop such that the requirement for vibration transmissibility in the primary resonance is met. The constraint conditions for this problem are as follows: (a) the parameter boundary where the forced vibration of the system is linear in the working frequency range; (b) the minimal stiffness ratio of the stop for given resonance transmissibility at different excitation amplitudes (or clearances); (c) the hysteresis set H that guarantees no jumping and no hysteresis in the primary resonance. It is obvious that condition (c) has been given in Eqn. 13. Hence, only the first two conditions are discussed hereafter. By eliminating /2 in the expressions of a l and G~2 in Eqn. 10, one obtains 64 ' ~2' ,2
O~2
71;2f 2
2
471;2;]2~,~(l + O.)2 a , .
(17)
The forced vibration of the system in the working frequency range is linear if the following inequality holds true "~/(1Aw2)2  (2g02w) + 5, (18) where 2 w is the ratio of working frequency to the natural frequency of the primary system. Given fmax, the critical value of Eqn.18, a parabola denoted by fma~ in Figure 5(b) can be f
< fmax =
determined from Eqn. 17. This is the parameter boundary of condition (a). 2
a
2
2
a
0
\H
2 0.0
0
Bl
_
0.2
2
0.4 a
0.6 l
0.8
1.0
2 0.0
0.2
0.4
0.6
a
0.8
1
(a) Transition set of unfolding parameters (b) Optimal parameter region Figure 5' Transition set and optimal parameter region in (a~, a 2) plane when(0=0.01,
4"1=0.2, Z w=2.75 and a a=1.2 at f  0 . 5 2
1.0
238
Hu, H.Y. and Wang, F.X.
Inthe fourth quadrant ofthe ( ~ l , a 2 ) plane, i.e., a 1 > 0 and a 2 < 0, a contour chart shown in Figure 5~b) can be obtained from Eqn. 17 if f = f = const., i = 1,2,... are taken respectively. It is easy to see from Eqn.10 that (al, a2) + (0, 0) when /a + +oo. Thus, the intersection of a contour curve with the hysteresis set // in Figure 5(b) gives the maximal stiffness ratio of stop, which guarantees no jumping and no hysteresis in the frequencyamplitude curve of the primary resonance. In addition, Eqn.10 implies that a 2 ~   c o when f + +oo or 6"+ 0 equivalently. To derive condition (b), it is necessary consider the maximal amplitude of the primary da resonance. By differentiating Eqn.8 with respect to 2 and imposing ~ = 0, one has
/]2
=
1+2p(a)4q(a)2. a
(19)
a
Substituting Eqn. 19 into Eqn.8 yields 2 a
[l+p(a)
1 2 q(a) a"
2
]q(a)  (
f2
) =0.
(20)
Solving Eqn.20 for the stiffness ratio, one obtains _~Tr f ]2 [q(a)]2 ~/11/a 2 = t[ + " 1}[arccos( 1 ) ]1. a2 2 2q(a) a a
(21)
If the acceptable dimensionless deformation of the vibration isolator is specified as a a when f = fa" Substituting the vibration amplitude a in Eqn.21 with a a , one obtains the minimal stiffness ratio of the stop that guarantees the vibration transmissibility for given f , . Let /~mm be the stiffness ratio. Then,/Jmi n results in a pair of (a 1, a 2 ) from Eqn.10. For other f , one can make use of the fact that the maximal amplitude of the primary resonance is approximate to the amplitude in Eqn.8 when p(a) + a(1  22) / 2 = 0, and proportional to f approximately. Thus, one obtains the change of //rain with variation of f from Eqn.21 by substituting a = a a f / f~ into Eqn.21, and then has a curve denoted by /dmin shown in Figure 5(b). As shown in Figure 5b, the region of optimal design parameters is the shaded one surrounded by the boundary of linear vibration, the boundary of minimal stiffness ratio of stop and the hysteresis set H. In summary, the parameters f and /a should be chosen such that the system does not have the excessive vibration transmissibility and jumping in the resonance frequency range. Hence, the design of dynamic characteristics of the system is independent, and will not affect the vibration transmissibility of the primary system designed according to the linear theory of vibration isolation in the working frequency range.
Computer Aided Design for Vibration Isolation Systems with Damped Elastic Stops
239
3.3 Numerical Verification The above design for the parameters of the damped elastic stop is based on the analysis of the primary resonance. Very often, a harmonically forced nonlinear system undergoes the subharmonic or superharmonic resonance, or vibrates chaotically when the system parameters are slightly perturbed or the initial state of the system does not fall into the basin of attraction of the designed state. Thus, it is necessary to verify the system design numerically from the viewpoint of both local and global nonlinear dynamics. The numerical verification includes the following steps. 3.3.1 Accurate Computation of Primary Resonance At this step, the excitation frequency is taken as the control parameter and the periodic vibration designed is computed by using the continuation technique developed in Wang and Hu (1996). If the analytic result greatly deviates from the numerical one, a new parameter combination in the optimal parameter region should be taken to check the results. If the new result is still not good, the constraint damping ratio (~ can be modified and the design in subsection 3.2 should be repeated. 3.3.2 Analysis of Global Dynamics One of the most important features of a nonlinear dynamic system is the coexistence of multiple steady state motions evolving from different initial states, say, three coexisting periodic motions or two coexisting periodic motions and a chaotic motion of large amplitude. It is obvious that the continuation technique is not able to determine all these motions. So, it is necessary to examine the effect of the initial states on the system dynamics when the system parameters are fixed. This can be accomplished by using the technique of interpolated celltocell mapping or the technique developed in Wang and Hu (1998). For the sake of reliability, it is better to make the analysis of celltocell mapping in a large region of the Poincare section. The size of convergence criterion in the celltocell mapping should not be too large in order that the results are reliable. If coexisting steady state vibrations are found in the celltocell mapping, the amplitude of each vibration should be examined. Once the expected vibration is not tolerable, the design has to be modified. Otherwise, the possibility of jumping phenomenon should be examined in order that the isolated equipment will not undergo dangerous shock due to the jumping when the excitation frequency varies. If there is a wide frequency range wherein multiple steady state vibrations exist, the robustness of the expected vibration should be examined. The robustness includes the stability of expected vibration, the stability redundancy of the vibration against the perturbation of system parameters and initial states. For this purpose, the following concepts will be used in the evaluation of a design. (a) The stability redundancy index is defined as the distance between the largest module of the eigenvalue of the linearized Poincare mapping and 1, i.e.,
240
Hu, H.Y. and Wang, F.X.
R, IIv[11
(22)
lOO%,
where v is the eigenvalue with the largest module. (b) The redundancy index against the perturbation of a system parameter p is defined as
Rp = ( 1 
=lI I lll>100% zXp
(23)
where v I and v 2 are the eigenvalues with the largest modules before and after the variation ofp. (c) The redundancy index of a basin of attraction is defined as
qy,. 11,' x 100%,
(24)
where A is the area of the Poincare section of concern,
YF is the fixed point, Yb is the
Ra =
A
point which is on the boundary of basin of attraction and most close to the fixed point. The value of Ra reflects the robustness of a basin of attraction. 3.3.3
Analysis of Bifurcated Periodic Motions
The result of celltocell mapping provides a set of initial fixed points for the continuation of periodic motions. At this step, larger meshes can be used so as to find the number and locations of the periodic motions efficiently. The technique developed in Wang and Hu (1998) is suggested to determine the fixed points of nodesaddle type since they can hardly be determined by using celltocell mapping. Once the initial fixed points are given, several kinds of continuation techniques can be used. Among them, the method suggested in Foal and Thompson (1991 ) is relatively simple. The subharmonic number of the resonance of a piecewise linear system depends on /a, the stiffness ratio of the elastic stop, but the occurrence of a subharmonic resonance depends mainly on damping. In the traditional design, small damping is usually used to guarantee low vibration transmissibility in high frequency range, and no damping is artificially arranged in the elastic stop. Thus, the subharmonic resonance is likely to occur. In the present design, large damping will be arranged in the elastic stop to avoid the subharmonic resonance effectively.
A NUMERICAL EXAMPLE For simplicity, an example is discussed here in the form of dimensionless parameters. The damping ratio of the primary system of concern is 4"0 = 0.01. It is required that the vibration
Computer Aided Design for Vibration Isolation Systems with Damped Elastic Stops
241
transmissibility at 2 w = 2.75 should be 0.2, the acceptable vibration amplitude should be 1.2 for f = 0.52, and the vibration transmissibility in primary resonance should be T < 5.0. The region of the optimal parameters for this system is shown in Figure 5b when the damping ratio of stop was set t o ( 1 = 0.2. If tWO unfolding parameters were chosen a s a I = 0.l 6 and a 2 = 1.05, there followed f = 2.73 and Iz = 2.06. The corresponding frequencyamplitude curve of the primary resonance is shown in Figure 6. It is easy to verify from Figure 6 that T ~ 2.81 < 5.0. In the continuation of the periodic motion, two turning points were observed when 2 e[0.5045, 0.5054]. In this frequency range, there exist three periodic motions. As shown in Figures 7, two of these periodic motions are stable, while the other is unstable. The multiplicity of these periodic motions exists in a very narrow frequency range. Even though there is a jumping between the stable periodic motions, the variation of the vibration amplitude in jumping is very small. The global dynamics of the system in the Poincare section of [4,0] x [2.4,2.4] was analyzed by using celltocell mapping approach at 2 = 0.5053. The result is shown in Figure 8, where two fixed points represent the abovementioned stable periodic vibrations. The redundancy indexes of these two stable fixed points are as following. The first fixed point corresponds to Rs = 72.5%, Re, 64.13%, Rf = 53.56% and R a =0.32%, and so does the second fixed point to Rs = 39.1%, R~ = 95.31%, Rf = 94.41% and Ra = 5.58%. It is easy to see that all the redundancy indexes of the second fixed point are relatively larger. Compared with the first fixed point, the second fixed point is more robust. This assertion was verified in the continuation of these periodic motions. Finally, the celltocell mapping was made for the system when the excitation frequency was fixed at 2 = 0.5 , 0 . 7 5 , 1 . 0 , 1.25, respectively. As expected, the only periodic motion found is the motion of period 1. In addition, no subharmonic resonance was found when the excitation frequency was set as multiplication of natural frequency of the linearized system. o.9 0.7
5
0.5
Yrr~
0.3
2:7.~:~_..~ x_ Aaalytica 0
0.5
1
0.1 1.5
2
Figure 6" Frequencyamplitude curves of primary resonance when (l = 0.2
0.5
0.505
0.51
0.515
0.52
Figure 7: The Poincare velocity versus excitation frequency curve when ( 1   0 . 2
242
Hu, H.Y. and Wang, F.X.
239
239
, , , , i , , , , i , , , , I , , , , i , , , , i , , , , i , , , , i ~ , , , i , , , , i , , ; ,
C3
Figure 8" The Poincare section of basins of attraction of two periodic vibrations when 2 = 0.5053 If a large stiffness ratio of the stop, say, ~ = 15.0 was chosen, a chaotic vibration was observed in the celltocell mapping. The Poincare section of corresponding strange attractor is shown in Figure 9. In this case, the system lost the function of vibration isolation. Hence, it is very dangerous to increase the stiffness of the elastic stop intuitively in order to limit the vibration amplitude. However, the vibration amplitude can be greatly reduced if the damping in the elastic stop is increased. For instance, when the damping in the above system was increased to (1 = 0.4, the maximal amplitude was greatly reduced as shown in Figure 10. Moreover, only a period 1 motion was found in the celltocell mapping. 3.5
0.8
i f
Nt.mi
2.5 ))
0.4
Y, rax
0.2
1.5 ,
.2
0.8
0.4
0
0.4
Figure 9: The Poincare section of a strange attractor when /~ = 15.0 and (1 = 0.2
0
D.5
,
,
,
,
1
,
,
i
,
1.5
,
\ 2
Figure 1O" Frequencyamplitude curves when ~t = 2.06 and £'1 = 0.4
Computer Aided Design for Vibration Isolation Systems with Damped Elastic Stops
243
CONCLUSIONS
A systematic design approach is suggested, by integrating several analytical and numerical approaches of nonlinear dynamics together, for the vibration isolation system with a damped elastic stop. This design approach has satisfied accuracy and reliabilit~ so that a great number of tests can be avoided. Both the numerical simulation and real test showed that the design of the vibration isolation system can greatly reduce the vibration transmissibility in the frequency range of resonance through the use of damping in the stop and keep very low vibration transmissibility in the working frequency range.
ACKNOWLEDGMENTS
This work was supported in part by the Natiofial Natural Science Foundation of China under Grant No. 59625511. The authors wish to thank Dr. Dongping Jin and Miss Wenfeng Zhang for their patient assistance during the research.
REFERENCES
1.
2. 3.
4. 5.
6.
7.
Foal S. and Thompson J.M.T. (1991). Geometrical Concepts and Computational Techniques of Nonlinear Dynamics. Computer Methods in Applied Mechanics and Engineering 89:381394. Hu H.Y. (1995a). Simulation Complexities in the Dynamics of a Continuously PiecewiseLinear Oscillator. Chaos, Solitons and Fractals 5:11,22012212. Hu H.Y. (1995b). Detection of Grazing Orbits and Incident Bifurcations of a Forced Continuous, PiecewiseLinear Oscillator. Journal of Sound and Vibration 187"3, 485493. Hu H.Y. (1996). Design of Elastic Constraints from Viewpoint of Nonlinear Dynamics. Chinese Journal of Mechanical Engineering 9:2, 135142. Wang F.X. and Hu H.Y. (1996). Calculating Accurate Periodic Solution of PiecewiseLinear Systems. Transactions of Nanjing University of Aeronautics and Astronautics 13:2, 142146. Wang F.X. and Hu H.Y. (1997). Bifurcations of Primary Resonance of Harmonically Forced Symmetrically PiecewiseLinear Systems. Journal of Nanjing University of Aeronautics and Astronautics 29"3, 283288. Wang F.X. and Hu H.Y. (1998). An Algorithm for Locating Coexisting Periodic Motions of Nonlinear Oscillators. Chinese Journal of Applied Mechanics 15:1, 105108.
This Page Intentionally Left Blank
CALCULATION OF THIN PLATES ON STATISTICAL N O N   U N I F O R M FOUNDATIONS
Yi
Huang,
Yuming Men and Guansheng Yin
Department of Basic Sciences,Xi'an University of Architecture and Technology, Xi' an 710055, P. R. China
ABSTRACT
In this paper,the calculation methods of reliability of structure (plate)  medium (soil) interaction are presented. Based on reliability theory, calculation methods of reliability of thin plates on statistical nonhomogeneous foundation are studied. Calculation formulas of reliability of elastic thin plates on Winkler, elastic halfspace and twoparameter foundation are developed. The calculation procedure is further demonstrated by examples of plates with four free edges and the plates with foursimply supported edges. KEYWORDS
Statistical nonhomogeneous foundation, thin plates, reliability foundation, elastic halfspace, twoparameter foundation
calculation,
Winkler
INTRODUCTION
Thin plates on elastic foundations are widely employed in engineering. Many applications,such as raft foundation, road pavement ,airport runway, etc. can be calculated by reducing them to thin plates on an elastic foundation~ that is the typical problem of structuremedium (soil) interaction. Concerning the calculation of plate on an elastic foundation, the certainty analysis method has been used, i.e. the physical parameters, geometrical dimensions of plate and foundation properties are considered as determinate factors, and actual variabilities of them are considered through socalled "safety factor". Actually, each kind of parameter has a relatively large variability because the foundation (soil) is highly dispersed, so it should be considered as a statistical nonhomogeneous medium, that is a key to reliability calculation of structuremedium interation. On the other hand, there are also variabilities in varying degrees for material properties and
245
Huang, Yi et al.
246
geometrical dimension of plates. This leads to uncertainty of the actual parameters that is difficult to represent by the classical "safety factor". Due to the lack of quantitative analysis of variability effect of all kinds of parameters on the deformation and strength of thin plates, unreasonable errors are caused in the design of foundation plate structures. Therefore ,in order to reasonably design the foundation plates,it is very important to study the effects of all kinds of uncertain parameters on the deformation and internal force of the thin plates. So it is inevitable to introduce reliability calculation of s t r u c t u r e  soil interaction. Recently, the fast development of reliability theory provides a powerful method to study the variable effects of the parameters on the structural internal force and deformation, and the reliability theory has been used in structural design. But up till the present moment there are not many papers dealing with reliability calculation of structuresoil interaction at home and abroad. In this paper, the reliability theory of platesoil interaction is discussed. Based on reliability theory, the reliability calculation of thin plates on elastic foundation is studied by considering the foundation as a statistical nonhomogeneous medium. The calculation procedure is further demonstrated by examples of a plate with four free edges and a plate with four simply supported edges. In fact , this paper presents reliability calculation method of structuremedium (soil) interaction. RELIABILITY CALCULATION HOMOGENOUS FOUNDATION
OF
THIN
PLATES
ON
STATISTICAL NON
Reliability calculation o f thin plates on Winkler foundation Winkler model is one of the simplest linear elastic foundation models. Although it has some theoretical drawbacks , it is still used in current engineering because of its simplicity. According to Winkler model, the deflection surface differential equation of elastic thin plates is
a,w + bk w 

(1)
b q
where A4(:) is symbol of the double Laplacian, W is deflection of the thin plate, k is
Et 3
modulus of the foundation, q is load of the plates, D  12(1 /~2) , h is thickness of the plate, 9 is Poisson's ratio of the plate material, and E is modulus of elasticity of the plate material. Once the deflection equation of plates is solved, the bending moment of plates can be obtained by the following equations: azW
M~ =    D ( ~
azW
+/~ Oyq)
azw azw My  D( 3yT +/14~_z) (TiA;
(2)
Calculation of Thin Plate on Statistical NonUniform Foundations
247
In the strength design of elastic foundation plates, ultimate moment of resistance is considered as the critical point. Thus, the reliability of elastic foundation plates is the probability that the maximum moment is less than the ultimate moment of resistance. Under normal design, construction and usage,the probability of reliability ps can be expressed as (3)
Ps = P ( M s > Mma~)
where M j is ultimate moment of resistance Mm~ is actual maximum bending moment. In the reliablity calculation of thin plates on Winkler foundation, the value of the modulus of subgrade is an important parameter, which is also the parameter of most variability. So, from a statistical nonhomogeneous viewpoint, it must be considered as a random variable in calculation. In addition, load p , modulus of elasticity E and plate thickness t , which should be considered as random variables,have a considerable effect on reliability. The ultimate moment of resistance, which depends on material properties and construction conditions, also has variability and should therefore be considered as a random variable. Poisson ratio of the plates has less variability and thus has less effect on calculation results. For this reason it can be considered as a constant in order to simplify the calculation. Other random variables can be regarded as independent of each other. Thus,the probability of reliability of thin plates on Winkler foundation can be approximately expressed in a function of the random variables, M s, E , t, q and k ~ that is (4)
Ps = f ( M s , E , t, q , k )
To simplify the expression, the above random variables can be written as Xi~ that is, X1 M s , X z = E , X3 = t, X4 = q, Xs = k. The mean and standard deviation of each variable
is/~x, and ax, respectively. Hence the limit state equation of thin plates is Z = g ( X 1 , X 2 , X 3 , X~, Xs)  M s  [ M ( X z , X 3 , X 4 , Xs)[ ....
(5)
where Ms is ultimate moment of resistance, which represent the resistance o~ the plates, and M is design bending moment, which represents the effects of actions of loads. When every random variable is of normal distribution, the reliability index equation can be obtained by using checking points method of firstorder second moment. 5
g(X;
, X ; , ... , X ~ ) 
___
/~
~
~.7..
(X~ 
~)
i1
og
~x
where p" is design checking points. Substituting Eqns. 2 and 5 into Eqn. 6 and noting the expression
(6)
248
Huang, Yi et al.
g(X(
, X ; , ... , X ; ) = 0
(7)
at the checking points, the formula of reliability index can be written in the following form fl=
{ g ( X ( )  { ( X [
Et3 1  t~x,)  12(1/~2){ME (X2"  / % )  + 
3
+ t(Xa" /~x,)] +
~xT(X;/~x,)},.}}/A
OM /=2
(8)
where M=
A 
O2W
a2W
{2 + [ ax,
(9M 2
(or
Eta
12(1 3 aM 2)
/j2)
a2W
o2W
(8')
?21( M2 1 8M2} 2 ~ + E ~2 ax,_ aM 2} },/2 i= 2
~
GXi
(8")
P"
Eqn. 8 is the reliability index equation of thin plates on Winkler foundation. The direction cosine of the normal op" to coordinate vector is. cosOx, =  
~ax,,
cOSOx, = 12(1
cosOx, =   12(1 
/~2)A t
~
~)A E + ~
ax,
ax,
12(1 ])A
~
ax,
(9)
E6 (aM) c°sOA = 12(1  /~2)A ~ ax~
Coordinates Xi" of the checking points are X;
= gx, + flax cosOx,
(10)
Solving Eqns. 7,9 and 10 simultaneously, the fl value of plate on elastic foundation can be obtained. In actual calculation, fl can be obtained only by iteration, of which the procedures are as,shown in Figure 1.
Calculation of Thin Plate on Statistical NonUniform Foundations
[Assume X,* (px, can
249
beuSed . .as. initialvalue)] . .
no l
[Calculate/3 by using limittstate equation g ( . )   0
i
¥
~tAre X;' and fl equal to the previously Calculated value or within acceptable error? h yes [/3 is the calculated reliability index, and .Xi* is the design . checking point. I" ] Figure 1
Procedures of iteration
In the above analysis, basic variables X~ are assumed of to be normal distribution. If X~ is not a normal random variable, it should be converted to an equivalent normal random variable before calculation. The calculation equations of mean/~x' and standard diviation ax' can be found elsewhere [4].
Reliability calculation o f thin plates on elastic halfspace foundation According to halfspace foundation model, the foundation is assumed to be a continuous, homogeneous, isotropic, completely elastic halfspace body. The mechanical properties of the foundation are functions of modulus of elasticity E, and Poisson's ratio/z, of soils. Although the solution of this model is much more complicated than that of Winkler foundation model, Eqn. 2 can still be used to describe the relationship between deflection and internal forces. Thus, the limit state equation of this model is similar to Eqn. 5. Among all parameters, M s, E, t and q are still considered as random variables. Although both modulus of elasticity E, and Poisson ratio/~, of soil are actually random variables, only E, is considered as a random variable during calculation in order to simplify calculation, and /~s is regarded as a determinate variable because it has less variability. Negl6cting the actual variability of/~, does not very affect the results of reliability calculation. When each random variable Xi is of normal distribution, Eqns. 7, 8, 9 and 10 can still be used as reliability calculation for the plates on the elastic halfspace foundation provided, replacing k in the above equations by E,.
Reliability calculation of thin plates on twoparameter foundation Two independent elastic constants are used to describe soil mechanical properties in the twoparameter model. According to the difference in the parameters used, this type of model can be further classified as FilonekoBorodich model, Hetenyi model, Pasternak model, Vlazov model, etc. In practice, the two parameters that characterize mechanical properties of soil are
250
Huang, Yi et al.
always variable with the discreteness of soil properties, no matter which model is used. So they are all actually random variables. However, whether both parameters are used as random variables in calculation depends upon practical circumstances. If it is complicated or timeconsuming to determine both parameters, only the parameter that represents the main soil properties is considered as a random variable, while the other one can be considered as a constant. For instance, subgrade coefficient k is considered as a random variable in FilonenkoBorodich model, while membrane tension T can be considered as a constant. The reliability calculation equations are thus exactly the same as those of Winkler model. However, It should to be noted that the subgrade coefficient here, k ,is different from that in Winkler model when 7" is introduced. If both parameters are relatively easy to determine, they can all be considered as random variables. For example, k and G in Pasternak model can be determined by means of ground testing,a relatively easy test and calculation. Both parameters are then considered as random variables. Let Xs = k , X 6 = G. Then Eqn. 8, the reliability index equation, becomes
Et 3 1 fl = {g(Xx)  {(X~"  /ax)  12(1 /~z) {M[~(X~  /~x,) + 3 ~ 4 y ( X ;  lax,)] +
aM  ~ ( X ~"  t~x, ) } p " } } / A
(11)
i=2
where
{ z + [ Eta ]z{{ Mz 1 3 M Z } z +_ A = ax, 12(1  !1z) EY + E ~2 ax,_ t ( 9MZ 3 OM z OM z 1/z i= z ffxiaX'
(11')
and the expression of M is the same as Eqn. 8' A new item is added in the equation of cosOx,
Et a 3M £1 cosux~ = 12 (1  ~z) A 3X6 ax~
( 12 )
g ( X [ , X2" , .. ,X~" ) = 0
(13)
and
Calculation procedures are the same as above. Obviously, the reliability index equations of thin plates on twoparameter foundation are almost the same as those on Winkler foundation except for the addition of the X6 item. T h u s , we conclude that for elastic thin plates on threeparameter foundation, the reliability
~(Xi*
index equations can be obtained only if items i=2
/~x,)
~ax,
and i=2
in the
Calculation of Thin Plate on Statistical NonUniform Foundations
3M
above expression are changed to
(X;'/~x,)
OM
and
ax,
251
2
, respectively,
iffiffi2
i=2
when three parameters are considered as random variables. Also, in the expression of Et 3 3M cosOxi , a new item cosOx, = 12 (1   / j 2 ) A OX7 ax, should be added . When reliability index/? is obtained, the probability of reliability p, can be calculated by p ,   f(fl)
(14)
where, p, is the probability of reliability of the plate, and ~o(fl) is the standard normal distribution function. If there are complex limit state equations or many random variables in practical problems, this computation procedure is very complicated to solve reliability index by using checking point method of firstorder second momnent. This can be avoid by using the iteration method of Fiesslur (1980). The principles can be found in reference [7] and the procedure is as follows: (1) Determine an expression for g ( X ) ( 2 ) Express the limit state equation g ( X ) as a function of yi by introducing X i  Yx, standardization variable y; = (i = 1,2, ...... ,n) ax (3) Determine expressions for all first derivatives of h ( y ) , h'i (4) S e t y = 0 a n d f l = 0 (5) Evaluate all h'; values (6) Evaluate h (y) (7) Evaluate standard deviation of Z from
(8)
Evaluate new values for y from
Y =
h', [fl + h ( y ) ] o; ~z
(9) Evaluate ~ = ~ / E y ~ (10) Repeat steps (5) to (9) until values converge. Fiesslur's iterative procedure can be used together with variables. This method is very convenient for multivariable problems or problems with complex relationships, because of its relative simplicity. EXAMPLES
As an first example, we discuss a reliability calculation of a fourfreeedge rectangular
252
Huang, Yi et al.
plate on t w o   p a r a m e t e r foundation ESl. Length of the plate is 2l = 8m and width of the plate is 2 b = 4m. The plate is acted upon by uniform load p . When soil thickness H ~c~, the two parameters, K and Gp , of the foundation are Eo7 K = 2(1  ~0)
Eot Gp = 4)'(1 + ~0)
and E, Eo  1   ~
v, Vo = l   v ,
where 7 is a measurement of nonlinear variation o{ normal strain of foundation,and E,, u, are elastic constants of soil material. Two parameters , K and G p , can then be expressed by E, and v,. Let Y and v0 be equal to 1. 5 and 0. 4 respectively; the relative rigidity between foundation and plate system is K , = 1.0. Solution of the deflection equation of plate is/8~ W(xy)
zra: zt'y = p0L(1  ~0) X 10 3 [408 + 2 . 3 7 c o s (  ~ ) + 17.4cos(f[) +E0
7rx
Try
ZO
"~t
(a)
+ O. 4cos ( ~ ) cos ( =v )
The values of mean and standard deviation of each parameter are given in Table 1. TABLE 1 VALUES OF MEAN AND STANDARD DEVIATION OF PARAMETERS Variable
Symbol
Mean(/%)
Standard deviation ( a x )
Unit
Ms
X1
500
50
KN. m / m
X2
2.25 X 107
3.6X106
KN/m 2
X3
0.35
Po
X4
300
30
KN/m 2
E$
Xs
635
190.5
KN/m 2
,
m
Calculation of Thin Plate on Statistical NonUniform Foundations
253
Poisson ratio of thin plates is/~ = 1/6 ,Poisson ratio of soil is u, 0. 2857. Both of them are constants; then
M
o~W .,r + oy
nl a,.Tg
o2W
x=O y=O
= 93 • 4913 X 10 4 Po E] = 93. 4913 X 10 ~ X~ Xs
Therefore, the limit state equation of the plate is Z = g(Xi)
= Ms
IDM
....
I = X1  O. 3436 X l O  ' X z X , / X s
(b)
The reliability index equation can be expressed in the following form by using checking points method of firstorder second moment as
t9
x;
{g(X;*)  {(X[  p t x )   O. 3436 X 1 0  4 [ ( ~ T ) ( X z " x~"
/*x,) 
x
/~x,) q O. 3436 X lO'(X2*z Xs" 4. ) (Xs* /~x~)
O. 3436 X I O  4 ( w ) ( X { ,~x 5
3} }/A
where
A 
2
2
X; , ix;x: 27 a~,+[ x;'
2 4
GX z
{ GX1
2 2 ]}I/2 GX s
At the checking point
g(X?, x;,
= X?
•. . , X d )

0. 3436
X
10'
IX;X;
[ x;
0
that is X~*Xs*  o. 3436 x lO4Xa * X ; = o
(c)
Each direction cosine is COS~x t
=
COSOx z = COSOx,~
123.696(X;/X~
)/A
=0
COSOx 4 = COS~x 5
50/A (d) 10.308
× 10 ~ ( x ; / x ~ ) / A
=  65.4558 ( X ; X ; / X ~
Coordinate of each checking point is
a)/A J
Huang, Yi et al.
254
X~ = 500 + 50flcosOx, X j = 2.25 X 107 + 3.6 X 106flcosO G Xa~ = 0.35 = 300 + 30pcosOx,
(e)
X;
Xs* = 635 + 190.5flcosOx s Using the system of Eqns. (d) (e) and (b) ,the results can be obtained by iterations ~the procedures are given in Table 2. TABLE 2 PROCEDURES AND RESULTS OF ITERATIONS
Iteration order
X;'
X;
5,00
2.25 X 107
0.35
300
635
0. 8032
485. 5304
2. 3718 X 107
0.35
306.342
514. 1684
0. 7922
490. 4025
2. 3519 X 107
0.35
305. 4781
503. 3878
0. 7922
490. 6520
2. 3510 × 107
0.35
305.4020
502. 8143
0. 7922
490. 6679
2. 3510 X 107
0.35
305. 3973
502. 7842
0. 7922
X;
Final result is ~3=0. 7922; the probability of reliability of the plate can then be obtained p, = ¢ ( P ) = 7 8 . 5 8 %
and the probability of failure is Pf=
1P,=
21.42%
It is evident that the probability of failure is larger. If we use the traditional safety factor method to calculate its safety factor,we obtain
K=
5OO = 1.37 365. 2441
Calculation of Thin Plate on Statistical NonUniform Foundations
255
It follows that the plate has enough safety capacity. That the variabilities of all kind of calculation parameters cannot be considered in safety factor is the main reason for different conclusions between the reliability and safety factor methods. Therefore unreasonable results arise. If the standard deviation of E is decreased to ax~  2.25 X 10SkN/m 2 , and that of E, to ax~ = 7 6 . 2 k N / m 2, then probability of failure of the plate will be decreased to
py  7 . 0 8 % As our second example, we discuss reliability of a fouredge simply supported rectangular thin plate on Winkler foundation. This plate supports uniform load q and has 1 length a   6 m , width b = 4 m and Poisson ratio/1 = ~ . Other variables are given in Table 3. TABI.E 3 VALUES OF MEAN AND STANDARD DEVIATION OF PARAMETERS
Variable
Symbol
Mean (ff~,)
Standard deviation( G )
Unit
Mj
Xl
5O0
5O
KN. m / m
E
X2
2. 304 X 107
2.88 X 106
KN/m 2
Xa
0.3
0.02
m
X4
400
40
KN/m 2
X~
30000
7500
KN/m a
q0
T h u s , the solution of a fouredge simply supported rectangular thin plate under uniform load q0 is Eg~
W=
16q0 zr2
oo
~
~~
~
,,,=,,a,,~,...,,,=,,a,s,...
sin(mn)x sin(
a m2
nTr
bn2 ) 2
)y
rnn[~r,D(_.a_~ + _~ +
(f) kI
This series converges fast. Here , only the first item is taken to illustrate the calculation procedure.
Huang, Yi et al.
256
Substituting Eqn. f into Eqn. 8, we find
=
M
i
~W
~W
I Oy~
+~~1
7r)s
= 16q0 ( b ) 2 + t ~ a zr__.y V,o ( + ~1 ~ +k ,sin ~rXsin a b [
i
(g)
Putting in known data into the formula and expressing it as a basic variable, the result is 1074.1X4 ~rx . Try M = 68.05X~X] + l O00Xs sin ~sm 4
X2X]X, So, Z = g ( X , ) = Mj
IMI ....
=
Mj
IDMIm,, ~  
Xi
Introducing the standardization variable y; = ~
,
ax,
X
~x,
O. 7392XzX] + 10. 8618Xs
1
(i = 1 , 2 , " ' , 5 ) , the limit state
equation can be expressed as
Z = h ( y ) = (y,ax, + px,) (yzax, + px~)(ysax~ + IJx~)3(yiax, + px,) O. 7392(y2ax~ + Zx2)(y3ax~ + Zx~) 3 + 10. 8618(ysax, + Zx~) Forming the derivative of this expression with respect to y~(i= 1 ~ 5 ) , w e obtain
3h (y) ~Yl
 ax, 10. 8618(ysax 5 + [Jxs)(y3ax, + tlx3)3(Y4ax, + Fx )ax, [0. 7392(yzax~ + ~x,)(y3ax, + [~x~)3 + 10. 8618(ysax, + ~x )] 2
~h(y) ~Y2 3h(y)
..
~Y3
(Yzax 2 + Px,)(Y3ax 3 + Px3)ax , O. 7392(yzax2 + ~x~)(y3ax~ + t%) 3 + 10. 8618(ysaxs + ~x,)
~h(y) ~Y4 Oh(y) ~Y5
+ /2x,)(ysax~ +/2x)Z(y4ax, + llx )(ysax~ + t2x )ax 3 [0. 7392(yzax2 + ~x,) (y3ax3 + ~x~) 3 + 10. 8618(ysax, + t~x )] z
3 X 10. 8618(y2ax,
..__
10. 8618(y2ax 2 + px2)(y3ax ~ + lax ) ( y , ax, + px )ax s [0. 7392(y2ax~ + [lx,)(yaax~ + / 2 x ) 3 + 10. 8618(ysax, + ~x )] 2
From Fiesslurrs method,the iteration results are given in Table 4.
Calculation of Thin Plate on Statistical NonUniform Foundations
257
TABLE 4 PROCEDURES AND RESULTS OF ITERATIONS Iteration order
where yi 
yl
y2
y3
y4
ys
0.0
0.0
0.0
0.0
0.0
1. 6540
0. 5431
0.8690
1. 0477
1. 0863
2.4624
1. 4704
0. 4236
0. 6841
1. 1301
1. 2420
2.3726
1. 4207
0. 4172
0.6721
1. 0937
1. 2741
2.3375
1. 4250
0. 4152
0. 6700
1. 1000
1. 2838
2.3474
1. 4242
0. 4155
O. 6693
1. 1000
1. 2872
2.3486
1. 4237
0. 4151
O. 6688
1.1000
1. 2876
2.3483
1. 4236
0. 4150
O. 6687
1. 0999
1. 2877
2.3482
Xi*
~
ax,
~x i
0.0
is the conversion variable.
Reliability index is ~2. 3482 ~therefore, the probability of reliability of the thin piate is p, = ~ f l )  99.0567 500 According to the safetyfactor method,safety factor K = 316.70= 1.58 The above result is approximate because only the first item of the series is taken to demonstrate the calculation procedure. CONCLUSIONS The reliability calculation method of thin plates on Winkler, elastic halfspace and twoparameter foundations is discussed in this paper. The equation of reliability calculation is established for three different models of foundation. In the reliability calculation, the ultimate moment of resistance M i , modulus of elasticity E , thickness t and load q0 are considered as random variables. Other random variables which represent soil properties should be selected according to the concrete model of the foundation. Furthermore~ the selection of random variables suggested in this paper is not unvaried~ they depend on specific conditions of practical engineering. For instance, high construction quality of foundation plates leads to little variability of structure dimensions ~the thickness t of the plate can then be treated as a constant. Otherwise, it should be treated as a random variable.
258
Huang, Yi et al.
The reliability index equation in this paper is only suitable for the calculation of thin plate structure on elastic foundation. As to the probability of reliability for soilfoundation structure interaction ,this topic will be discussed in future papers. REFERENCES
]1 ] Jinzhang Zai. (1994). Foundation Analysis and Design o f TaU Buildings, SoilStructure Interaction Theory and Their Application, China Construction Publishing House, Beijing. [2]Selvadurai, A. P. S. (1979). Elastic Analysis o f SoilFoundation Interaction, Elsevier Scientific Publishing Co. [3]Zhaohong Zhu, Binggang Wang and Dazhi Guo. (1985). Mechanics Calculation for Pavement ,Peoples Traffic Publishing House, Beijing. [4] The State Standard of the People's Republic of China. (1992). Unified Design Standard for Reliability o f Engineering Structures (GB5015392) ,China Planning Publishing House, Beijing. [5]The State Standard of The People's Republic of China. (1994). Unified Design Standard for Reliability o f Railway Engineering Structures (GB5021694) . China Planning Publishing House,Beijing. [6] Tianyi Zou. (1991). Reliability o f Bridge Structures, People's Traffic Publishing House, Beijing. [7]G. N. Smith. (1986). Probability and Statistics in Civil Engineering. Nichols Publishing Company, New York. [8]Yao Sheng and Yi Huang. (1987). A Free Rectangular Plate on The T w o  Parameter Elastic Foundation, Applied Mathematics and Mechanics 8..4, 317 329. [9]S. Timoshenko and S. WoinowskyKrieger. (1959). Theory o f Plates and Shells. M cGra w  Hill Book Co. [10]Zhiming Tan and Yao Zukang. (1991). Reliability Analysis o f Concrete Pavement, Journal of Tongji University 19 : 2, 157 165. [11]Dajun Ding and Zhongde Liu. (1986). Calculation Theory and Method o f Beam on Elastic Foundation, Nanjing Technology Institute Publishing House, Nanjing. [12]Zhilun Xu. (1979). Mechanics o f Elasticity, Peoplers Education Publishing House, Beijing. [ 13 ] Shaozhong Xi, Shiying Zheng. ( 1981 ). Applied Mechanics o f Elasticity. China Railway Publishing House, Beijing.
COMPUTER SIMULATION OF STRUCTURAL ANALYSIS IN CIVIL ENGINEERING Jiang JianJing, Guo WenJun and Hua Bin
Department of Civil Engineering, Tsinghua University, Beijing
ABSTRACT In this paper, structural analysis and failure process in civil engineering are discussed. The main contents are: (1) The philosophy of computer simulation in structural engineering; (2) The mathematical model for engineering problem; (3) The visualization of numerical results. Some simulation examples are presented. KEYWORDS Computer Simulation, Structure, Civil Engineering, Disaster Risk.
INTRODUCTION Prototype Test or Site Observation
Constitutive Law (Physical Model)
d Simulation Analysis ,,
~
Numerical Method (mathematical Model)
No
Graph&Image System (Visual Equipment)
~J
Application ]
Modify Constitutive Law or Mathematical Model Improve Numerical method of Graphic System
Figure 1" Philosophy of computer simulation in structural analysis
259
260
JianJing, Jiang
et al.
Computers are used widely to simulate the objective world, including natural phenomenon, system engineering, kinematics principles and even the human brain. Though civil engineering is a traditional trade, computer simulation has been applied successfully, especially in structural analysis. Three prerequisites are needed to perform structural analysis: (1) Constitutive law of specific material, which can be obtained by smallscale test; (2) Effective numerical method, such as finite element method (FEM), direct integration, etc.; (3) Graph display software and visual system. Figure 1 shows the philosophy of computer simulation in structural analysis. The following parts give a comprehensive explanation of several aspects.
SIMULATION OF STRUCTURAL FAILURE Structural behavior under various loading conditions and environment is of great importance, especially its collapse procedure and ultimate loading capacity. When the structure form is very special, we usually resort to model experiment in order to determine the characteristic of the structure. Yet the model usually is small because of the constraints of space and equipment, thus can not reflect the behavior of real structure. If we want to study the influence of a parameter, series analogous experiments must be done, which is very timeconsuming and expensive.
(a) ~   ~
(b) ~ 3 0 . 3 I ~
(c) P=33.6KN
Figure 2: Development of micro fracture
Computer Simulation of Structural Analysis in Civil Engineering
261
Taking advantage of computer simulation, we can carry out "fullscale experiments" without worrying about time and budget. Simply changing several input parameters, we can obtain the influence of these parameters. Sometimes it is impossible to do real experiments, then computer simulation is of more significance. For example, it is obvious that we cannot repeat the procedure of accident of nuclear reactor safety shell. Using computer simulation, we can perform reversal analysis of the accidents to determine the cause and unfolding procedure of accidents. Other examples are highspeed collisions of car and structure collapse under earthquake, which can be performed only by computer. Under high speed loading, structural response is very quick, therefore what we get is the final result. An extremely contrary situation is the evolution of the earth's crust, which may take millions of years. Both can be simulated by computer by changing time rate according to request. Figures 2 [41 and 3 [11 show the development of microfracture and dynamic response of a frame, respectively.
1 t=7.86
s
(a) (b)
Figure 3" Dynamic response of flame
APPLICATION IN HAZARDS PREVENTION The history of human being is a history of hard battle with various hazards, such as flood, fire, earthquake and so on. Because hazards are unreplicable for experiments, computer simulation has become more and more popular and of great value. Many simulation systems have been developed successfully. For instance, software has been developed and preloaded with landform, topography and surface features of a flood area. Given flood standard and specific location of burst, computers can demonstrate submerged areas at different time on screen, which are calculated according to water quantity, speed and area. People can view the gradual inundation process and, in turn, work out flood prevention and personnel dispersion programs. Fire prevention software is another successful example. Using this software, we could simulate spreading of fire in forest and building, which give guidelines to firefighting. Figure 4 [31 shows
262
JianJing,
et al.
Jiang
the distribution of debris after an earthquake.
giNCK t~trt~il/ l~lltgl~:8
Figure 4" Distribution of debris after earthquake
1~:'.':':5:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: iii:i:i::::!, iiii~!iiiiiiiiiiiiiiiiiiii~iiiii~!~iii~i~i!i~iii~ii~iiiii~iii~i!iiii~iiiii~iiiii~!i!i!iiiiiiii!!iii!i! iiii!iii!~i~iiii!!iii!iiiii i.lj.:]"!iiii!i~iiiiiiiiiii :~,~.A~,, :'.0~~ ,t'.. ' ~ b~ ................ , ................." .................""..................................................................................... ii:".................. :~:~:~:~:~:~:~:;:~:~:~:~:~:~:~:~:~:~:~:~:~:~:~:~:~;:~:;:.:;~:~:~:~:~:~:~:;:~:~i~......... :~:~:~:~; , . ,~,~i~.;." ..............!i~::~!~::!::~::!!i::~::~::~iiii::ii~::i!~i~::~i~!~!::~::~ii!::~::~::~::~iii~: :~!:~
.....•""'''"""""""""""""""'."'."::::::::::::::::::::::::: .'.'.'.'.'.'.'.'.'.'.'.'.'.'.'.' ~.~.~i"~ ~ ~ ~g "" ~ ::::::::::::::::::::::::::::::::::::::::::::::::::::::: ~~:::::::::::::::::::::::::::::::::::::::::::::::::::::: i i ii~ii::!~ ==================.... ===============~!=~,, ~'~ ~¢J," ~ ~ ~o~_.l
~i~i~:~:.~i~i~:~i~.~i~i~i~~:~~:~.:~.~...... !!:.'.'i~!ii:.::~ii::i:.i'.ii::iiiiii':iiii':i!::i::ii!:.::!::i:.::
~i:':"
".'.'.'.'.'.'.'.'.'.''.'.'.'.'.'.'.'.'.'.'.':'.'.'.'.'."" ........ ""~." "'"'"":'.... ::': "::~:i:!:~:!:i:!:!:!:i:~:!:!:i:~:~:~:~:~:~:~:~:~:~:! ====================================================================================i iiiiiii~" :::::::::::5:::: i~i~
~ .......
:'i~!!:i!~"ii!!i!i!~iiii!............ ~i~ ..'..~..'~.~ .......
~.}.'i!/iii,.
":. . . . . . . " ~
1~
~ ~
' : .~
.
.
.
.
.
.
.
.
.
.
.
.
.
',
iiiiiiii!iiiili:N ii!iiiii!1i1!iiii i !i!iiiiii!i!iiiiii?iii!!!iiii!iiiiiiiiii!ii!iiii!!iiiii!i!iiiii Figure 5: Submerged area
Computer Simulation of Structural Analysis in Civil Engineering
263
APPLICATION IN R O C K AND SOIL ENGINEERING Construction in rock engineering is underground, hence it cannot be observed directly. Computer simulation which can reveal its inner procedure, is of great value. For example, during the excavation of underground construction, there is always collapse, which can be solved by thorough geological survey to find the strike of faults, crevices and joints. Through smallscale experiments we can determine the mechanics of rock body and joints, which can be stored in computer for later use. Besides finite element method (FEM), there is discrete element method (DEM). The behavior of elements of DEM is similar to FEM ones in equilibrium, while elements of DEM will move under external force and gravity when equilibrium is lost until they get new equilibrium. In the analysis of underground rock structure and stability of slope, the structure is divided into elements along crevice and joint faults In simulation of excavation, the upper and side part elements may lose equilibrium and fall down; this procedure can be shown on screen and thus we can obtain the cavein area and provide reliable instruction to anchor design Computer simulation is also applied in research fields such as seepage of underground water, deposition of fiver silt, settlement of foundation and so on For example, a simulation so,ware of fivermouth deposition has been developed. When given the condition of fiver mouth, it can show the deposition rate of different size of silt and accumulating thickness, which give instruction to port design and fiver course dredging Figure 5t21 shows a simulation of bridge pier collapse
~°°°'°° . . . . . . . . . i .. B ~ . " :
: : ' . % . . / Y .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
°°~
°
~a× 01 ~ p l e c ~ m e n t
O. ~ 7 5 0 + 0 2
Figure 6: Collapse of bridge
264
JianJing, Jiang et al.
SIMULATION SYSTEM OF TEACHING EXPERIMENTS During the teaching of reinforced concrete element for students of civil engineering, demonstrating experiments is a key part, which can make the knowledge easier to understand and strengthen perceptual knowledge of students. Yet, the real element failure test is very tedious, timeconsuming and expensive, not affordable by university teaching budget. With the aid of computer simulation method, we can build a graphic environment to simulate the experiment. After input geometric and physical data, the students can observe the procedure of element failure, details of the inner process, and other changes. Compared with teaching test, this method can initiate students' activity, giving them the opportunity to participate in the experiment as well as saving large quantity of work, material and time. Three aims can be achieved through simulation teaching: 1)The student can get a clear and vivid knowledge of the phenomenon of element failure and its characteristics. Loading...
50
[~ ]
I~
I
II,
1
/,( Loading...
+'
) i
'
100
'o
[ 9f; )
( Figure 7" Demonstration of test 2)They can select different parameters of elements, such as section dimension, concrete strength, reinforced ratio, and the location of force and the influence of these parameters on failure shape, ultimate loading capacity. 3)When doing the simulation, there are instructions on the screen, as if there is a teacher. In the simulation instruction system developed by the author, there are two types of experimentthe example experiment and free experiment. The former takes the role of a real instruction experiment usually used nowadays in teaching, and the students' job is to select what type of failure they want to see. Then everything is done by computer automatically.
Computer Simulation of Structural Analysis in Civil Engineering
265
With the selfexperiment, students have an opportunity to experiment by themselves. They design their own elements; the system automatically analyses element failure and demonstrates the failure procedure. By selecting different values, the students can acquaint themselves with the influence of these parameters.
CONCLUSIONS Computer simulation has achieved great success in many fields, including structural analysis. Through the above description, three conclusions can be reached: (1) Along with the rapid improvement of CPU speed and update of hardware, the computer has been used not only as a tool of scientific calculation but also in structural analysis, hazard prevention, construction management and failure simulation. (2) Much hard work, expense and time can be saved by using computer simulation. A combination of computer simulation and experiments will be the main tools of engineers. (3) Computers can give miscellaneous dull data a vivid and lifelike form, they will play a more important role in teaching and management. (4) Compared with the high level of hardware, simulation systems lag behind. More mature and businesslike simulation systems need to be developed.
REFERENCES 1 Jiang J.J. (1991). The computer simulation system for R. C structure process. Computational mechanics and application (in Chinese), (2), 3238. 2 Wei Q.(1990). The fundamental principle of distinct element method in geotechnics, Ph.D. dissertation, Tsinghua University. 3 Jiang J.J. (1990). Application of computer simulation to civil engineering. Computational technique and computer application (in Chinese), (2), 1619. 4 Liu G.T. and Wang Z.M.(1996). Numerical simulation study of fracture of concrete materials using random aggregate model. Journal of Tsinghua University (Science and technology), (1), 8489.
This Page Intentionally Left Blank
A MIXED FINITE E L E M E N T F O R L O C A L AND N O N L O C A L PLASTICITY Xikui Li National Laboratory for Structural Analysis of Industrial Equipment, Dalian University of Technology, 116023, P.R.China
ABSTRACT This paper develops a mixed finite element with one point quadrature and hourglass control in local and nonlocal (gradient) plasticity for pressure dependent and independent materials at large strains, which allows for its application to the modelling of strain hardening and softening (strain localization) behaviors. The two yield strength parameters of strain hardening/softening materials not only depend on the internal state variable but also on its Laplacian. The evaluation of the Laplacian is based on a least square polynomial approximation of the internal state variable around each integration point. To derive the consistent element formulations in pressure dependent elastoplasticity, a natural coordinate system in the stress space and a new definition of internal state variable are introduced. Numerical examples are given to demonstrate the performance of the mixed element, particularly in preserving ellipticity as strain softening behaviour is incorporated into the computational model.
KEYWORDS mixed finite element, local plasticity, gradient plasticity, pressure dependent materials, natural coordinates, internal state variable, strain localization, consistent formulations
INTRODUCTION Owing to the advantages of the one point quadrature mixed finite element with hourglass control in both accuracy and efficiency, it has been widely implemented in finite element codes and used to solve a variety of practical engineering problems with success. Many efforts have been devoted to develop and improve this type of mixed elements since it was first published
267
268
Li, Xikui
by Flanagan and Belytschko (1981). Jetteur and Cescotto (1991) developed their one point quadrature mixed element for the vonMises elastoplasticity with strain hardening at large strain and reported excellent results obtained with cheap cost in computation. Localization of deformation into narrow bands of intense strain caused by strain softening is a characteristic feature of plastic deformation and is a common occurrence in many pressure independent and dependent materials. Numerous attempts to simulate the localization behaviour by using classical (local) plastic continuum theories have been unsatisfactory. It has been realized that the effective and radical measure to remedy this situation is to introduce the regularization mechanism into the continuum to preserve ellipticity or hyperbolicity of the governing field equation for the quasistatic problems or the dynamic problems respectively. Among diverse approaches to introduce the regularization mechanism are: nonlocal continuum theory (Bazant et a1,,1984, de Borst, 1992), Cosserat continuum theory (Muhlhaus et al., 1989) and rate dependent continua (Needleman, 1988). The basic idea of the nonlocal gradient plasticity is to include higher order spatial gradients of the effective plastic strain in the yield condition, de Borst et al. (1992) and Sluys et a1.(1993) presented the formulations and algorithms for the gradient plastic continuum in a finite element context. A finite element method for gradient elastoplastic continuum was presented by Xikui Li et al,(1996), in which the Laplacian of the effective plastic strain at a quadrature point is evaluated on the basis of a least square polynomial approximation by using the values of the effective plastic strains at neighbouring quadrature points. This nonlocal approach allows to satisfy exactly the nonlocal consistency condition of the yield function at each quadrature point, whereas the consistency condition is only enforced in a weak form and is not satisfied at each iteration but only at the end of a load step in de Borst et a1.(1992) and Sluys et a1.(1993). The objective of the present work is to develop a nonlinear version of the element for both gradient and local plasticity of pressure dependent material models, such as the DruckerPrager and the modified vonMises models, at large strain. The element is formulated not only for the plane strain and the axisymmetric solid, but also for the plane stress state. To formulate the nonlocal consistent compliance matrix and the nonlocal consistent integration algorithm in a concise and numerically efficient manner for closedform implementation, the socalled natural coordinate system in the stress (strain) space (Duxbury & Xikui Li et a1.,1996) is introduced. The separation of the plastic strain into its deviatoric and hydrostatic components due to the introduction of the natural coordinates provides an opportunity to define a new internal state variable that is capable of capturing different postyield curves in tension and compression simultaneously. To analyze the geometrically nonlinear problem by utili~ng the present mixed element, the corotational formulation (Jetteur and Cescotto, 1991) is adopted and the corotational Cauchy stress tensor and its energetically conjugated strain measure are employed.
TIlE MIXED FINITE ELEMENT Let us start with the variation of the fimctional H for the HuWashizu principle in the form
6n= j'[6~(c~ o)+~d(~  Vu)+5(Vu) ~o]dA, Ae
(1)
Mixed Finite Element for Local and Nonlocal Plasticity
269
where the stress , the strain ~, the displacement gradient Vu and the elastic modulus matrix C are referred to the original coordinates. For the twodimensional case, we have er =[~, ~;y g~ ~;~y]; Vu r =[u,,., uy,y u,,, u:,,y +u,,~]
6 r =[6~ 6, 6~ X,y];
(2)
The natural coordinate system in the stress (strain) space is introduced in such a way that the deviatoric and the hydrostatic components of the normal stress (strain) components are split and 6, e and Vu are transformed into their counterparts 6, ~ and Vfi referred to the natural coordinates as follows 8=T6;
Vfi=TVu
~ =T~;
with~yr=[% % 6mX,y]; ~ r = [ e , engm~y ]
(3)
where the transformation matrix T is defined by
 /4g
T =
1/4
(4)
l 1/xf~O 1/xf30 1/X~O OlJ
According to the displacement field equivalent to the classical bilinear function of isoparametric coordinates ~ and 11 , U i   a 0;  ~  a j
Ae x i Jra 3
(i=l, 2)
n
the
(5)
where A~ is the area of the element, the displacement gradient can be derived as n
(6)
Vfi = T ( B + h bF ) q = T B q
where qr = [u r v r] is the nodal displacement vector of the element and
Br
=
[b 0 0hi ,x
0
,y
by
0
b~
1
br,~ = ~ [ Y 2 4
00h ]
r
"
'" h b =
h ..v 0 h~o
" F=
'
1 Y31 Y42 Y13] ,"
b,yT : ~ _ ~ _ IX4 2 X13 X24 X31 ] Ig
2A e = x31Y42 + x24Y31 , X/j = Xi  Xj ; y~ = yi  yj A, c?h 3h 1 h "y = 3 y ; y = ~ [ h  ( h r x j ) b ° ] h = 4 ~ r / , • h~,  3 x '• h r [1 1 1 1];
X T1 . _ [X 1 X2 X3 X4] ;
X2T = [Yl
Y2 Y3 Y4]
The stress ^ and the strain ~ are chosen in the 'optimal incompressible' modes as
(7)
270
Li, Xikui
=¢y+h~& ; ~ = ~ + h ~
; with h T = h.y
h,y
o Ool
0
(8)
where ~r = [~i ~i1 F'mF_,xy], ~xT = [~1 ~2] and ~T = l% % ~m xxy], &r = [~ ~ ] . The constant field ( ~ , ~) is uncoupled with the antihourglass mode ( & , U ) . With substitution of expressions (6) and (9) into equation (1) we may integrate each term in the equation and obtain ~)1I= A e[5~ T( ~ _ ~) + 5~T (~ _ TBq) + ~SqrBTTr~] + ~SUr [4GH~ _ 2H~¢ ]
+~&r [2H~;~  H*Fq]+SqTFTH *r& Here C
c R 4×4
(9)
referred to the natural coordinates for isotropic elasticity can be given by
C= TCT r =diag(2G, 2G, 3K, G)
H=LH~,
H.
;
(10)
h~ ha dA , ; H * =
Ho =
hrThb
dA e
(11)
where G and K are the elastic shear and bulk moduli. From the stationary condition of (9) and the arbitrariness of the variations 6~,5~,~q,~ x and ~&, we obtain the constitutive laws, the element straindisplacement relations, the internal nodal forces F and the linear stiffness matrix of the element ~= TBq ; 2 ~ ~' = H1H*Fq ; &=2G~'; K = A e B T C B +GFTH*TH1H*F F = A e B r T r ~ + FTH *r~x ;
~=C~
(12) (13)
THE GRADIENT DEPENDENT PLASTICITY Two pressure dependent elastoplastic models for the gradient plasticity are particularly considered. The yield functions for the gradient DruckerPrager (GDP) and the modified vonMises (GMVM) models can be given in the form (Duxbury and Li, 1996) F = q + A(~P,V2~P~y m + B(~P,V2~p) F = q2 + A ( ~ v , V 2 ~ p ) % + B(~v,V2~V)
(14)
(15)
where the effective deviatoric stress q and the hydrostatic stress ~m are defined by 1
q = (~I~T e{~) 1/2 ", {~m=~(l~x+~+l~z) ", with
P = diag( 3, 3, 0, 6)
(16)
A and B are the current material parameters defined as 2Sind~ (~ P, V2~ p) A = ~(3Sind~(~v,V2~V))
_6c h (~ v, V2~ p ) Cosd~ (~ p, v2~ p) '
B=
,f3(3Sind~(~v,V2~V))
(17)
Mixed Finite Element for Local and Nonlocal Plasticity
271
for the DruckerPrager criterion, where c h is the cohesion and ~ the internal frictional angle. For the modified vonMises criterion, c p A : 45(~;(~ ,v~,) I 9 . 4 ! 1 7 . 8 ~ml 16.1~ 14.5~ 12.9~ 11.3~
8.0 96"_74 1 4.8 3.2
W~
@
i~iiii!:i}iiiiiiiii:::
i}i}i il}i i !ili!i i i! ~,~!iii!i!iiiii!!!i!i)iiiiiiiiii:' ~ ........~:~:iiiiiii!ii~i!!i}~ii:;i~iiii)iiiiii~il;!?iii;!
iliiiiii}iiiiiiii!}!i}iiiililiiiii)iiii!!iiiiiiiiii ,,,, 1 6
o  ! _5 3.1 4.6 6.2 7.7 9_3  ! o.lq 12.4 I 3.9 15.5 17.1 I
Z I4 W
500
U.
o
o
6
1o
15
20
25
3o
35
MEASLREMENT /
40
45
50
55
60
65
Figure 10. Monitoring data for a cable stay from the Tabor Bridge in Slovak Republic
Monitoring of Cable Forces Using MagnetoElasticSensors
347
REFERENCES
1. Bartolli, G., Chiarugi, A. and Gusella, V., "Monitoring Systems on Historic Buildings: The Brtmelleschi Dome", ASCE J. Struct. Engrg., pp. 663673, June (1996). 2. Bozorth, R. M., Ferromagnetism, D. Van Nostrand Co., NY, USA, (1951). 3. Jarosevic, A., Fabo, P., Chandoga, M., and Begg, D. W., "Elastomagnetic method of force measurement in prestressing steel", lnzinerske ls Stavby, v. 7., pp. 262267, Bratisala, Slovakia, (1996). 4. Krauss, J. D., Electromagnetics, McGraw Hill Inc., (1992). 5. Kvasnica, B., and Fabo, P., "Highly precise noncontact instrumentation for measurement of mechanical stress in low carbon steel wires", Meas. Sci. & Tech., pp. 763767, (1996). 6. Kwun, H., and Teller, C. M., " Detection of fractured wires in steel cables using magnetostrictive sensors", Materials Eval., pp. 503507, April (1994). 7. McCurie, R., Ferromagnetic materials: Structure and Properties, Academic Press, (1994). 8. Moon, F., Magnetosolid Mechanics, Wiley InterScience, (1984). 9. Sablik, M., and Jiles, D., "Coupled magnetoelastic theory of magnetic and magnetostrictive hysteresis", 1EEE Trans. On Magnetics, v. 29, no. 3, pp. 21132122, (1993). 10. Shahawy M., and Arockiasamy, M., "Field Instrumentation to study the Time Dependent Behavior in Sunshine Skyway Bridge. I", ASCE J. Bridge Engrg., pp. 7686, (1996). 11. Shahawy M., and Arockiasamy, M., "Field Instrumentation to study the Time Dependent Behavior in Sunshine Skyway Bridge. II", ASCE J. Bridge Engrg., pp. 8797, (1996). 12. Stanley, R. K., "Simple explanation of the theory of total magnetic flux method for the measurement of ferromagnetic coss sections", Materials Eval., Jan (1995).
This Page Intentionally Left Blank
FRACTURE ESTIMATION: BOUND THEOREM AND NUMERICAL STRATEGY ChangChun Wu
QiZhi Xiao
ZiRan Li
Dept. of Modern Mechanics, Univ. of Science and Technology of China Hefei 230026 China
ABSTRACT A bound analysis of fracture parameters is proposed. It is found that the lower bound for Jintegral can be obtained by a compatible displacement finite element method. On the other hand, the upper bound of the/*integral, as the dual of the Jintegral can be obtained by an equilibrium finite element method. To avoid the difficulty of designing equilibrium finite element models, a hybrid stress model is modified by incorporating a penalty equilibrium constraint. Moreover, a relative error measure formula for J and I* is suggested. Numerical examples on different crack and loading configurations are presented to verify the validity of the bound theorems.
KEYWORDS Fracture, Pathindependent integral, Upper/lower bound, Hybrid finite element
INTRODUCTION The Jintegral had been proven to be equivalent to the release rate of the strain energy rI(lt,) with respect to the crack area (Rice,1968). Hence, the bound of the numerical solutions for the Jintegral, if it exists, may possibly be established using the assumed displacement finite element method which is founded on H(lti). On the other hand, the/*integral proposed as the dual counterpart of Jintegral (by Wu et al.) has also been shown to be the release rate of the
349
350
Wu, ChangChun
et al.
complementary energy He(a0. ) with respect to the crack area. A natural conjecture is that the bound of the numerical solutions for the/*integral, if it exists, may possibly be established using the assumed stress finite element method which is founded on II¢(~j). In computational fracture mechanics, the estimation of upper/lower bounds for fracture parameters becomes a matter of great significance as one cannot obtain an accurate or reliable solution no matter what experimental or numerical method is used due to the complexity of fracture problems. The bound problem consists of two aspects: (1) theoretically, the existence of an approximate upper/lower bound for a certain path integral, and (2) if it exists, numerically, the evaluation approach. In addition, the error measure should be considered once the bound solutions are obtained.
DUAL PATHINDEPENDENT INTEGRAL
For a given plane crack system with actual states of stress, strains and displacement (ajj, e~j, u~ ), the Jintegral (Rice, 1968) can be defined as:
j_
dH
 Jaa
=
Ir [W ( e 'j ) dx 2 
a ,j n j (
~ 1 )ds]
(1)
To find a dual integral of J, we introduce the Legendre transformation
W(u~ ) + B(o'~j ) = a~je~j
(2)
into (1); then an alternative path integral can be derived
_,, (: J) : I,..
 B(o,.,.
 o,.,.,.,.N~'
ds}
(3) It is easy to verify that the present I* is a pathindependent integral. The I*integral can also be defined as a complementary energy release rate(Xiao,1996):
I*

dH ~
&r ,j
d a  Sr [B(a,j )dx 2 + u,  ~ l njds +
8 c3cj
(u,a,2)dx
Y
l
(4)
Fracture Estimation: Bound Theorem and Numerical Strategy
For an equilibrium stress field,
351
crij,j = O, then it can be verified that (3) is equivalent to (4). In
numerical calculations, the expression (3) is especially recommended due to its simplicity and the absence of derivative of stresses.
BOUND THEOREMS
Corresponding to the dual integrals J and 1', the following bound theorems can be established for a certain linear or nonlinear elasticity crack system with homogeneous displacement boundary constraint, i.e. Lower Bound Theorem for J: For the given cracked system, if u, and ~ are respectively the exact displacement cmd the approximate one based on the mh#mum potential energyprinciple, the approximate Jintegral will take a lower bound of its exact one, Le. _< J(,,,)
(s)
Therefore, the lower bound o l d can be obtained by displacement compatible elements. Upper Bound Theorem for I*: For the given cracked system, if o"o. and ~j are, respectively, the exact stresses and the approximate one based on the minimum complementary energy principle, the approximate I*integral will take the upper bound of its exact one: I*(o'o) > I'(cru)
(6)
Therefore, the upper bound of I* ccm be obtained by stress equilibrium elements.
In the case of linear elasticity, a proof for above theorems was presented by Wu, Xiao and Yagawa(1998). Furthermore, in the case of nonlinear elasticity, including deformation theorybased plasticity, the theoremsstill hold. As an illustrating example, the lower bound theorem for J is proved here. For the given nonlinear crack system with homogenous displacement boundary constraints, it can be verified that for the actual solutions 06, or0.) :
n(,,,) =
= v ( ¢ o.)
(7)
ku(,,,)
(8)
352
Wu, ChangChun
et al.
where V(crq) and U(u~) are the complementary energy and the strain energy, and the finite constant k = V ( % ) / U(u~ ) > o
For linear elasticity, k = 1. Let
(9)
u~ = u~ + 8u i , then we have
n(a, ) = n(u, ) + 8rl + 8:n(8,, )
(10)
As the displacement finite element method is an implementation of the potential energy principle, the first and second variations are respectively
817= 0 and 62 FI(6u~ )  Iv A(Su, )dV Hence, Eqn.(1 O) becomes H(~'~ )  H(u~ ) = Iv A(fu~ )dr In accordance with the definition of J,
n(.,)) _  ~ a( . ~ A(Su,)av
j ( ~ , )  j ( u , )   a,aS
(a)
Considering an actual status of the given system, the Jintegral must be positive, i.e. J(u~)>__O, and Eqs.(8) and (9) must be satisfied. Thus we have d j ( u . )   T a n ( . . ) 
d
d
ku(.,. )] 
k f A(..)dv
>_o
(b)
Observing that the employed compatible displacement elements can keep the strain energy to be positive definite, the comparison of (a) and (b) results in
"L
 ~
A(a.,)dV _< 0
Thus the inequality(5) holds.
(1~)
Fracture Estimation: Bound Theorem and Numerical Strategy
353
NUMERICAL STRATEGY As for J, its lower bound can easily be obtained by using conventional isoparametric elements. For I °, however, its upper bound should be estimated by stress equilibrium elements. Unfortunately it is hard to get a reliable equilibrium model for 2D and 3D problems because of numerical difficulties, such as rank deficiency, displacement indeterminacy, etc. cannot be avoided .We face the problem of how to implement the upper bound theorem for I °. It is observed that the stress equilibrium element is based on the complementary energy formulation 1Ic(o), while the stress hybrid element based on the Reissner formulation l'I R(o,u). However, I[R(o, u) is identical to FI c (o). Thus a hybrid model may degenerate into equilibrium model when the stress equilibrium equations are enforced to the hybrid element. For an individual element, let
•
HR a =
. or(Du)_~o So
?V
(12)
A generalized functional can be created in the manner of ( Wu and Cheung, 1995).
/7~6 = //~ _ 2
, (D r o)r (D r o)dV
(13)
In Eqn.(13) the penalty factor ct >0 is taken to be a large constant such that the homogeneous equilibrium condition D r o  0 is enforced to the element in a leastsquares sense. Recalling the 4node plane hybrid stress element, termed as PS, proposed by Pian and Sumihara(1984), the assumed element stress trial functions can be expressed, in terms of the element coordinates (~, rl), as
[a,:
t
(14)
0 0 1 a,b,rI a3b34d[fls
where the coefficients a, and manner:
o ,lft
b~ are related to the element nodal coordinates (xi, yi) in the
354
Wu, ChangChun et al.
a1
b1
a:
b:
1 1 l =;
1
1
1
1
1 1
1
I"
a3
b3
_1][ x,
Y,
_ l/x
lj[x x4
(15) Y4
By substituting the stress (13) and the bilinear displacements uq  N(~, r/)q into the functional (12), we have e riRG(p,q)=flTGq_
fiT( H + a H )fl E P
(16)
After condensing/3, the element stiffness matrix is now
K e  G r ( H +  Ha
E
In Eqn.(17) the matrices penalty matrix
P
) 1 G G and H
(17) are identical to those
of PS element, while the
Hp  ~v" (Dr~)T (DT~) dV
(18)
In such a way, PS hybrid element, is developed into a penaltyequilibrating model, termed as PS(ot), in which the stress equilibrium equation is imposed by the penalty function method.
ERROR MEASURE Let 6u,  u ,  u s be the displacement error induced by using the assumed displacement finite elements. Then, in accordance with the lower bound theorem (5), the relative error for the Jintegral must be J(Su, )  J ( ~ )  J(u, ) < 0
(19)
[J(Su~ )] J(u, )  J(ff~ )
(20)
and the absolute error
Let
8a~j  cr'~j a~j be the stress error induced by using the assumed stress finite elements; in
accordance with the upper bound theorem (6), the relative error for I*integral must be
Fracture Estimation: Bound Theorem and Numerical Strategy
I* (8o,j) I* (5,j) I*
o
355
(21)
and
II * (8a~j)[ = I * ( ~ j ) 
I * (o~j)
(22)
Corresponding to (20) and (22), the relative errors can respectively be
[J(fu~)[/J(u~)
and
[I*(6o,j)[/I*(o~j)
expressed as (23)
The strength estimation for a certain structure can be carried out once the upper and lower bounds are obtained, and the approximate strength is usually taken to be a combination of the bound solutions. In order to measure the error of the fracture parameters given by finite element methods, a relative dual error for J and I* is defined as:
Aj_I. =
(24)
J(u, ) + I * (o~j )
In Eqn.(24), the sum of reference solutions is (25) Observing the small quantities
J(Su,) < 0, whereas
I * (6o~j) > 0, the last term in (24) can be
ignored, we then have
J(u, ) + I * (o~j ) ~ J(K, ) + I * (~j )
(26)
Substitute (20), (22) and (26) into (23), and note that the actual J and I* are identical in value; we finally obtain .~
_I*
Aj_,. 
(aiJ)  J(u~) I * (au) + J(u'~ )
(27)
The above relative dual error formula only depends on the approximate solutions of J and I*, so as to be easily used in nonlinear fracture estimations. Obviously, the error will vanish when the adopted finite element meshes become more and more fine.
356
Wu, ChangChun et
al.
NUMERICAL TESTS In numerical calculations, the well known 4node isoparametric element Q4 and the present penaltyequilibrium element PS(ot) are employed to estimate J and I* respectively. The center cracked panel CCP with uniform stretching load o'®(Fig. la) and the single edge cracked panel SECP with uniform stretching load ~®(Fig.2b) are calculated. Only a quarter or half of a specimen needs to be considered due to the symmetry. Three finite dement meshes and two integral paths are shown in Fig.2a,b,c. The specimen material consists of Young's modulus E =1.0 and Poisson ratio v =0.3. The distributed stretching load or® = 1.0. For linear elastic crack problems, all the solutions of J or I* will be transferred to the stress intensity factor K+ for convenient comparisons. The reference solutions of K+ are offered by H.L. Ewalds and R.J.H. Wanhill (1984). To inspect the convergence behavior of the solutions of J and I*, three meshes with different densities and two independent integral paths are simultaneously considered for each specimen. From the results shown in Fig.3a,b and Fig.4a,b, it can be seen that the solutions of J by Q4 always converge to the exact one from a certain lower bound. On the contrary, the solutions of I* by PS(ct) always converge to the exact one from a certain upper bound. All the numerical solutions demonstrate the bound theorems presented in the paper. The error formula A+_+. in (27) is implemented to measure the relative error of the bound solutions for CCP. The results are listed in Table 12. The results given by the formula As_I. in (24) are also shown in parentheses in the tables for comparisons. It is found that, independent of the selection of meshes and paths, both A+_~. and A+_+. always offered almost the same results. These numerical tests exhibit the efficiency of the present approximate error formula (27).
CONCLUSIONS Lower and upper bound theorems have been established for J and I* respectively such that the estimation of fracture parameters can be carried out by means of
J ( ~ )