- Email: [email protected]
Adaptive finite element analysis of structures under transient dynamic loading using modal superposition
Finite Elements in Analysis and Design 31 (1999) 255—272
Adaptive finite element analysis of structures under transient dynamic loading using modal superposition A. Dutta , C.V. Ramakrishnan *, P. Mahajan Regional Engineering College, Silchar 788010, India Applied Mechanics Department, Indian Institute of Technology, Delhi 110016, India
Abstract In this paper, the error estimation and adaptive strategy developed for the linear elastodynamic problem under transient dynamic loading based on the Z—Z criterion is utilized for 2D and plate bending problems. An automatic mesh generator based on “growth meshing” is utilized effectively for adaptive mesh refinement. Optimal meshes are obtained iteratively corresponding to the prescribed domain discretization error limit and for a chosen number of basis modes satisfying modal truncation errors. Numerous examples show the effectiveness of the integrated approach in achieving the target accuracy in finite element transient dynamic analysis. 1999 Elsevier Science B.V. All rights reserved. Keywords: Finite element; Adaptive; Error estimation; Mesh generation; Stress recovery
1. Introduction Estimates of errors and adaptive analyses of practical finite element problems are subjects of great importance if confidence in results is needed for engineering applications. Without the assessment of the reliability of the results it is hardly reasonable to use numerical methods like the finite element method in such safety-sensitive areas as shape optimization of machine parts, construction of aeroplanes or dimensioning of nuclear power plants. Therefore, the a posteriori error estimation is now considered to be nearly as important as the finite element analysis itself. Using finite elements, the linear dynamic transient response is usually solved either by mode superposition or by direct integration schemes. The direct integration method is very highly useful for solving nonlinear problems. Procedures describing adaptive time stepping are available in [1,2] and also spatial mesh adaptation [3]. Wiberg and Li [4] have proposed a postprocessed type of a posteriori estimates in space and also in time when direct integration is used for dynamic
* Corresponding author. 0168-874X/99/$ — see front matter 1999 Elsevier Science B.V. All rights reserved PII: S 0 1 6 8 - 8 7 4 X ( 9 8 ) 0 0 0 6 1 - 4
256
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
response evaluation. It updates the spatial mesh and time step so that the discretization errors are controlled. This process is not only time consuming but also gives rise to new errors as nodal values need to be interpolated from the previous mesh to the newly generated mesh whenever a mesh is changed. But for large structural problems involving a large time domain, the mode superposition is widely used. The accuracy obtained by the mode superposition is usually dominated by the accuracy of the orthogonal modes being used and proper representation of spatially distributed loads by the number of modes selected for the modal analysis. Thus, we find that the following errors are introduced. These are (a) domain discretization error, (b) modal truncation error, (c) numerical error in eigen modes and eigen values, and (d) truncation error in numerical integration. While a proper value of time step eliminates truncation error in numerical integration to a great extent and an appropriate tolerance limit evaluates eigenmodes and eigenvalues accurately, domain discretization error and modal truncation error need elaborate study. The work reported so far on discretization error estimation under dynamic loading and using the modal superposition scheme are very few. Wilson and Joo [5] have arrived at the final mesh using Ritz vector as the basis of transformation. In their investigation, the authors have made use of modal participation and amplification factors and obtained error estimates based on Babuska’s criterion using amplified Ritz modes. This procedure is a bit cumbersome and does not give a measure of the actual error in the transient response. Cook and Avrashi [6] have discussed the procedure for estimating the discretization error of the finite element modelling as applied to the calculation of natural frequency of vibrations. Meshes are obtained corresponding to each mode. Dutta and Ramakrishnan [7] proposed a measure for discretization which is a logical extension of Zienkiewicz and Zhu’s [8] error criterion by involving time integration to consider the variation of response with time. Using this error measure, an adaptive mesh refinement strategy is proposed which yields good control over the discretization errors in transient dynamic analysis. An optimal mesh is achieved iteratively wherein meshes are refined on the basis of error indicators so that the discretization errors are within the prescribed limit for a chosen number of basis modes satisfying the modal cut-off criteria, when the modal superposition method is used. In this paper, we have adopted the adaptive strategy given by Dutta and Ramakrishnan [7] for arriving at an optimal mesh for a prescribed domain discretization error limit and for a chosen number of basis modes satisfying modal cut-off criteria for accurate evaluation of response under transient dynamic loading. A 2-D quadrilateral automatic mesh generator [9] is used for carrying out adaptive refinement of the finite element mesh. The adaptive strategy is extended to estimation of errors and adaptive analysis for plates and shells where the well-known degenerate shell element by Ahmad et al. [10] is ued for evaluation of response under dynamic loads.
2. Error measure for domain discretization The dynamic response of a linear structural dynamic system that is subjected to general dynamic loads is given Mz(#Kz#Cz"q(t). The initial conditions are z(0)"z and z(0)"z .
(1)
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
257
Here z is the vector of nodal displacements, M is assembled consistent mass matrix" L Me where n is the number of elements, C L K is assembled stiffness matrix" Ke, C C is assembled damping matrix" L Ce, (2) C q(t) is external excitation (time dependent). FE stresses are calculated as p"DBz,
(3)
where D is the elasticity matrix and B is the strain displacement matrix. The approximate solution p containing discretization errors differs from the accurate solution pH and the difference is the error e . N Thus, e "pH!p, (4) N Zienkiewicz and Zhu [8] suggested a simple approach to get better approximation of p by using a projection tecnique or nodal averaging procedure. In the projection technique, the use of weighted residual requirement for equality between pH and p leads to the computation of an accurate solution of pH which is identical to the least-squares smoothing technique of Hinton and Campbell [11]. Nodal averaging, though very simple, gives very good convergence for most of the problems. A recent and elegant technique of better stress recovery is done using superconvergent patch recovery method by Zienkiewicz and Zhu [12]. In the recovery process, it is assumed that the accurate nodal values pN H belong to a polynomial expansion pH of the same complete order ‘p’ as C that present in the basis function N and which is valid over an element patch surrounding the particular assembly node considered. Such a ‘patch’ represents a union of elements containing this vertex node (Fig. 1). This polynomial expansion will be used for each component of pH and C we get pH"Pa, (5) C where P contains the appropriate polynomial terms and a is a set of unknown parameters. For a general plate and shell problem, the polynomial can be expressed as a function of x, y and z. However, since only plate problems are solved in this paper, x and y are considered only as the parameters in the polynomial expression. Thus, it can be written as P"[1, x, y, x,2]
(6)
and a"[a , a , a , a ,2]2. (7) The determination of the unknown parameters a of the expression given in Eq. (5) is best made by ensuring a least-squares fit of this to the set of superconvergent points existing in the patch
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
258
Fig. 1. Computation of superconvergent nodal values: (䉭) Gauss point, (䢇) nodal values determined by recovery procedure, (䉺) patch assembly point.
considered. To do this we minimize L F(a)" (p (x , y )!pH(x , y )) C G G C G G G L " (p (x , y )!P(x , y )a), (8) C G G G G G where (x , y ) are the coordinates of a group of sampling points, n"mk is the total number of G G sampling points and k is the number of the sampling points on each element m (m "1,2,2,m) of H H the element patch. The minimization condition of F(a) implies that a satisfies L L P2(x , y )P(x , y )a" P2(x , y )p (x , y ). G G G G G G C G G G G This can be solved in matrix form as a"A\b,
(9)
(10)
where L L A" P2(x , y )P(x , y ) and b" P2(x , y )p (x , y ). (11) G G G G G G C G G G G Once the parameters a are determined, the recovered nodal values of pN H are simply calculated by insertion of appropriate coordinates into the expression for pH. The stresses in the nodes inside the C patch are evaluated using a from Eq. (5). The stresses at boundary nodes can be determined using interior patches.
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
259
For the time-dependent problem, we propose #e # ;100%, g" #u # where g is the overall discretization error with 2#e# dt . #e# " ¹
(12)
(13)
The integration is carried out numerically using Simpson’s rule. Also 2#u# dt #u# " ¹
(14)
where ¹ is the duration of response, #e# is the ‘energy norm’ of the error for the whole domain
L , #e#" #e# G G where
(15)
#e# " G
e2D\e dX N X N
,
(16)
the energy norm of the error in the element i, #u#"W #uM ###e# X ,
(17)
where
L #uM #" #uM # G G and
#uM # " G
X
p2D\p dX
,
(18)
the energy norm of the FE solution. The complete procedure of discretization error measure is given below in an algorithmic form: Table 1
3. Estimation of error due to modal truncation Modal truncation errors are introduced in the response calculation because of the introduction of a reduced system. Since for a specified number of modes, discretization errors have now been quantified and an appropriate mesh satisfying the prescribed tolerance limit of gN can be properly designed, the measure for control of error due to modal truncation can be carried out effectively. Inclusion of higher modes alone without taking care of domain discretization error leads to
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
260
Table 1 Algorithm for domain discretization error (1) Solve MVX"KV where V is a matrix of eigen vectors and X is matrix of eigenvalues (2) Calculate dynamic response for t"0,¹ +(2.1) For i"1, m nos. of eigenvector solve 1 + y (t)" R f (q)e\DGSGR\O sin u (t!q) dq G u G z t " y (t), G G where u is the ith damped natural frequency, f is the damping ratio corresponding to mode G G i and f is the load vector corresponding to mode i and can be written as f "2q(t). G G G (2.2) for e"1, n elements + Retrieve ze t from z t where z is the displacement vector corresponding to the element level CR for j"1, number of Gaussian point + calculate pH "DBHze t CR calculate smooth stress pHH, CR , for e"1, n elements + for j"1, number of Gaussian point + eHe "pHH!pH N R CR CR #e# " X(eHe )2D\ eH dX CR N R e NeR
#uM # " X(pHe )2D\ pH dX CR N R e NeR , for the whole struture #e# "[ #e# ] R CR
#uM # "[ #uM # ] , R CR , (3) Overall domain discretization error 2#e# dt #e# " R ¹ 2#u# dt R #u# " ¹ #u# "W #uM # ##e#X R R R #e# ;100% g" #u# (4) Discretization error at element level For i"1, n elements #e# + m " GR is a refinement indicator at element level G eN K
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
261
Table 1 Continued where
eN "gN K
#uM # ##e# R R n
is the allowable element error g is the permissible overall discretization error 2m dt m " G G ¹ If g)gN and If m (1 for i"1, n STOP G h Else h " G G (m ) N G h is the size of element i of the previous mesh, G h is the size of element i of the new mesh G Endif, Go to step 1
erroneous results [7], since there is no guarantee that the chosen mesh is suitable for the inclusion of higher modes. Given a finite element mesh, it is possible to control the modal truncation errors in finite element analysis using the ratios of MPF, modal participation factor and its maximum value (MPF ) and MSF, modal significance factor and its maximum value (MSF ) as demonstrated
by Dutta and Ramakrishnan [7]. The modal participation factor for the ith mode can be written as 2p/u. The ratio of MPF and MPF is b and b (d where d is the permissible limit. The G
modal significance factor for the ith mode can be written as 1 uN 2p ; where b " G u u ((1!b)#2f b ) G G G G G and u is the exciting frequency and the ratio of MSF and MSF is b (d.
4. Automated meshing algorithm A general formulation developed by Sinha and Ramakrishnan [9], termed ‘growth meshing’ which is characterized by rule based local selection of elements resulting in an inward propagation of the mesh and includes several existing methods of unstructured gridding, is utilized for adaptive refinement of the domain. The 2-D quadrilateral meshes are generated mainly following the strategy adopted by Zhu et al. [13], where two triangles which share a common side are combined. The main feature of the quadrilateral generation is that first two triangles which have a common side are generated and then the triangles are combined to form a quadrilateral. The process
262
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
continues until the domain is fully covered by non-overlapping quadrilaterals. Furthermore, it is rule driven and offers flexibility through the incorporated constraints along with the possibility of adaptive modification. It is very effective in achieving a rapid gradation effect over irregular and complicated domains. However, one drawback of this type of quadrilateral generation by combining triangles is that the quality of the raw quadrilateral mesh produced is usually not good and mesh quality enhancement procedures are adopted following Zhu et al. [13] and Lee and Lo [14]. In order to capture the rapid variation of element size which commonly occurs in adaptive refinement meshes, the mesh generator employs the previous quadrilateral mesh as the background mesh. Nodal spacing values of the previous mesh obtained through error estimation is utilized for calculating node spacing at a point of the new mesh to be generated. The nodal spacing at a point is calculated as [15] L (1/r)d (19) d" G G G, L 1/r G G where r is the distance between the node to be generated and the nodal point of the previous mesh, G d is the nodal spacing at the nodal point of the previous mesh and ‘n’ is the total number of nodes. G The nodal points of the previous mesh which are neighbouring close to the point under consideration influence in arriving at the value of the nodal spacing at that point. For the generation of an entire quadrilateral element mesh over the domain, the boundary of the domain is transformed into a polygon with an even number of sides. A set of uniformly distributed sample points is first placed along the boundary line and node spacing on each sample point is calculated by interpolating its value obtained from the background mesh. The boundary nodes are then generated following a simple formula given by Lo [16] which allows the position of the new nodes to be calculated one by one along the boundary. The interior nodes are generated following the strategy given by Peraire et al. [17]. Starting from an active front chosen from the boundary segments, the nodal spacing values at the background mesh are utilized such that the elements with acceptable aspect ratio and minimum distortion are generated.
5. Computation of optimal mesh using modal superposition An overall strategy for accurate computation of dynamic response is given by Dutta and Ramakrishnan [7] for plane stress/plane strain problems. Following the same strategy, the number of basis modes are first obtained for a coarse mesh to begin with, which will satisfy the cut-off criteria based on b Od or b Od adopted. Discretization errors at element level m (i"1, number G of elements) are determined for those number of modes. Average nodal spacing values are calculated using m values and mesh is then adaptively refined based on those nodal spacing values. G The number of modes required for the refined mesh are determined satisfying the modal cut-off criteria. The overall domain discretization error g and the discretization error at element level m are G determined for chosen number of modes. Iteration is carried out until the mesh is an optimal one satisfying the prescribed domain discretization error limit gN corresponding to the number of modes chosen on b or b limit.
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
263
5.1. Numerical example A number of plane stress and plate bending problems have been solved to demonstrate the adaptive strategy and use of the automatic mesh generator in obtaining the optimal meshes. Example 1 (Plane stress plate with a hole). A square plane stress plate with a central hole subjected to harmonic loading of exciting frequency uN "5 rad/s is considered. Because of symmetry, only a quadrant could be considered for analysis and Fig. 2a is assumed as the starting mesh. The limits of MSF ratios and MPF ratios are chosen as 5e-2. The limit of domain discretization error gN is considered as 1%. The geometric and material properties of the plate are as below. Plate: 100 in;100 in;1 in with a central hole of radius of 5, E"0.3e5 lbf/in, l"0.3; o"1.0 lbm/in. For the problem with uN "5 rad/s, 15 eigenmodes are selected to satisfy the modal cut-off criteria corresponding to the initial mesh (Fig. 2a). The final mesh is obtained after one iteration. The automatic mesh generator described earlier is used for adaptive refinement. The exciting frequency is very low for this example and lies between first and second natural frequencies. The results are presented in Table 2. The final mesh is shown in Fig. 2b. NIGEN is the number of eigenmodes. A similar exercise is repeated with exciting frequency uN "25 rad/s and Fig. 3a is assumed as the starting mesh. The problem with higher exciting frequency uN "25 rad/s requires 19 basis modes corresponding to the initial mesh (Fig. 3a) in order to satisfy the modal cut-off requirement of b and b O5.e-2.
Fig. 2. Plane stress plate subjected to harmonic loading (u"5): (a) Mesh 1, 12 elements, g"5.38%; (b) Mesh 2, 127 elements, g"1.12%.
264
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
Table 2 Optimal mesh for the problem of a plate with a hole subjected to harmonic loading of uN "5 rad/s Mesh
Limit of b &b
NIGEN
g (%)
Fig. 2a Fig. 2a
5.e-2 5.e-2
15 17
5.383 1.122
Fig. 3. Plane stress plate subjected to harmonic loading (u"25): (a) Mesh 1, 12 elements, g"17.84%; (b) Mesh 2, 348 elements, g"2.04%; (c) Mesh 3, 393 elements, g"1.07%.
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
265
Table 3 Optimal mesh for the problem of a plate with a hole subjected to harmonic loading of uN "25 rad/s Mesh
Limit of b &b
NIGEN
g (%)
Fig. 3a Fig. 3b Fig. 3c
5.e-2 5.e-2 5.e-2
19 23 25
17.842 2.037 1.075
The final mesh is obtained after two iterations and results are shown in Table 3. The exciting frequency lies between the seventh and eighth natural frequency. The final mesh (Fig. 3c) corresponding to this problem is quite different from the one (Fig. 3b) corresponding to uN "5 rad/s because of the integration of different numbers of basis modes in the analysis. Example 2 (Plate problems). A number of plate examples are chosen with different boundary conditions, loading pattern and geometry of the plates so that the error estimation and adaptive analysis can reflect the variation in mesh gradation. A plate clampled on all the four edges and subjected to suddenly applied pressure load is considered. A quadrant is considered because of symmetry. The initial mesh (Fig. 4a) has nine elements. Modal truncation error is controlled by keeping the limit of b and b as 1.e-3 and the overall domain discretization error tolerance gN is kept as 4%. The geometric and material properties of the plate are as below. Plate: 300 in;300 in;3 in; E"3.e7 lbf/in; l"0.316; o"7.324 e-4 lbm/in. The optimal mesh (Fig. 4c) for the clampled plate subjected to suddenly applied pressure load is obtained iteratively following the strategy described in Section 3.8 and the results are shown in Table 4. Two iterations are required to arrive at the mesh for which g is less than the prescribed tolerance limit for overall domain discretization errors and the number of modes selected for the analysis satisfy the modal cut-off criteria. Intermediate mesh and the number of basis modes required at all the stages are mentioned in Table 4. The total number of elements corresponding to the final mesh (Fig. 4c) are 234 and the mesh clearly shows the refinement requirement near the clamped boundaries where the stress gradient is high. Similarly, the clamped plate is loaded with a suddenly applied point load at the centre of the plate. The discretization pattern for such a problem under static loading is quite clear; but a study of this problem under transient dynamic loading is made to understand as to how the mesh gradation vary. The optimal mesh is determined using the prescribed limit as mentioned above. Initial mesh used has nine elements as shown in Fig. 5a. The optimal mesh is determined using the prescribed limit as mentioned above. The optimal mesh (Fig. 5c) for the clamped plate subjected to a suddenly applied point load at the centre of the plate is also obtained iteratively and the results are shown in Table 5. Table 5 Two iterations are required to arrive at the optimal mesh satisfying prescribed error limits. The intermediate mesh (Fig. 5b) and the final mesh (Fig. 5c) show that the refinement requirement is
266
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
Fig. 4. Symmetric quadrant of clamped square plate (a/t"100) with uniformly distributed suddenly applied load: (a) Mesh 1, 9 elements, g"31.62%; (b) Mesh 2, 89 elements, g"6.65%; (c) Mesh 3, 234 elements, g"3.45%.
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
267
Table 4 Optimal mesh for a clamped plate subjected to a suddenly applied pressure load Mesh
Limit of b
Fig. 4a Fig. 4b Fig. 4c
1.e-3 1.e-3 1.e-3
&b
NIGEN
No. of elements
g (%)
7 9 10
9 89 234
31.62 6.65 3.45
more near the point of application of the load which is quite obvious. The effective stresses are calculated along the centroidal line parallel to the ½-axis for all the three meshes and these stresses are plotted as shown in Fig. 6 at some specified time steps. The stress plots clearly show that the stress distribution is affected by the finite element mesh refinement and stresses stabilize to the correct value when domain discretization errors g corresponding to a finite element mesh is low. Table 5 shows that the Mesh 1 has quite a high value of g when compared to Meshes 2 and 3. Thus, the differences in stresses corresponding to Meshes 2 and 3 are much lower when compared to that of Mesh 1 except near the centre of the plate. The refinement requirement near the centre of the plate will be more stringent as the concentrated load is acting at that point and hence the stress picture near that zone can be further improved accordingly. The study shows that the adaptive mesh refinement can yield a better stress distribution over the entire problem domain and the process is fully automatic. A simply supported skew plate subjected to a suddenly applied pressure load is also considered. Material properties are the same as mentioned above. The overall domain discretization error tolerance gN is kept as 9% and b O1.e-3. The analysis is carried out for the full plate of size 300;300;3 in 3. The initial mesh has nine elements and the mesh is shown in Fig. 7a. This problem is considered since it is expected that stress concentrations will occur near the corners where the edges meet at an obtuse angle and hence such a problem is very interesting from the point of view of adaptive analysis. The optimal mesh (Fig. 7c) for the simply supported skew plate subjected to a suddenly applied pressure load is obtained following the same strategy adopted above and the results are shown in Table 6. For proper representation of the boundary conditions along the skew edges, the following considerations are made. Simply supported boundary conditions along an edge are maintained by allowing only the rotation about that edge and restraining all other degrees of freedom. At a node ‘i’ (Fig. 8) along the skew edge, the rotations about normal to that skew edge can be written as b cos h!b sin h. Using the penalty approach, this amounts to adding a term to the potential energy which can be written as C(b cos h!b sin h), where C is a large number.
268
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
Fig. 5. Symmetric quadrant of clamped square plate (a/t"100) with suddenly applied point load at the centre of the plate: (a) Mesh 1, 9 elements, g"32.37%; (b) Mesh 2, 167 elements, g"5.01%; (c) Mesh 3, 313 elements, g"3.21%.
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
269
Table 5 Optimal mesh for a clamped plate subjected to a suddenly applied point load at the centre of the plate Mesh
Limit of b
Fig. 5a Fig. 5b Fig. 5c
1.e-3 1.e-3 1.e-3
&b
NIGEN
No. of elements
g (%)
10 12 13
9 167 313
32.37 5.01 3.21
Fig. 6. Effective stress (lbf/in) vs. distance along the centroidal line: (a) At 1/8th cycle; (b) at 1/4th cycle; (c) at 3/8th cycle; (d) at 1/2th cycle.
270
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
Fig. 7. Simply supported skew plate (a/t"100) with uniformly distributed suddenly applied load: (a) Mesh 1, 9 elements, g"39.57%; (b) Mesh 2, 175 elements, g"13.69%; (c) Mesh 3, 317 elements, g"8.82%.
Table 6 Optimal mesh for a simply supported skewed plate subjected to a suddenly applied pressure load Mesh
Limit of b
Fig. 7a Fig. 7b Fig. 7c
1.e-3 1.e-3 1.e-3
&b
NIGEN
No. of elements
g (%)
8 12 15
9 175 317
39.59 13.69 8.82
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
271
Fig. 8. Node ‘i’ along skew edge.
This can be rearranged as
C cos h 1 [b b ] 2 !C sin h cos h
!C sin h cos h
b
C sin h
b
. The terms C cos h, !C sin h cos h and c sin h get added to the stiffness matrix, for every node on the incline and the new stiffness matrix is solved for the unknown variables. Similar considerations are made for translational degrees of freedom. Two iterations are required to arrive at the final mesh. The intermediate mesh (Fig. 7b) and the final mesh (Fig. 7c) show that the refinement requirement is more near the corners where the two adjoining edges meet at an obtuse angle.
6. Conclusion The error estimation scheme for plane stress/plane strain and general plate and shell problems subjected to transient dynamic loads is computationally efficient and easily implementable. Local error indicators combined with adaptive meshing is very useful, leading to solutions of specified accuracy. The overall strategy is very effective in achieving high level of accuracy introducing only the required additional degrees of freedom during mesh refinement.
References [1] O.C. Zienkiewicz, Y.M. Xie, A simple error estimator and adaptive time stepping procedure for dynamic analysis, Earthquake Eng. Struct. Dyn. 20 (1991) 871—887. [2] L.F. Zeng, N.-E. Wiberg, X.D. Li, Y.M. Xie, A posteriori local error estimation and adaptive time-stepping for Newmark Integration in dynamic analysis, Earthquake Eng. Struct. Dyn. 21 (1992) 555—571. [3] L.F. Zeng, N.-E. Wiberg, Spatial mesh adaptation in semidiscrete finite element analysis of linear elastodynamic problems, Comput. Mech. 9 (1992) 315—332. [4] N.E. Wiberg, X.D. Li, A post-processed error estimate and an adaptive procedure for the semidiscrete finite element method in dynamic analysis, Int. J. Numer. Meth. Eng. 37 (1994) 3585—3603. [5] K.J. Joo, E.L. Wilson, An adaptive finite element technique for structural dynamic analysis, Comput. Struct. 30 (6) (1988) 1319—1339.
272
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
[6] R.D. Cook, J. Avrashi, Error estimation and adaptive meshing for vibration problems, Comput. Struct. 44 (3) (1992) 619—626. [7] A. Dutta, C.V. Ramakrishnan, Error estimation in finite element transient dynamic analysis using modal superposition method, Eng. Comput. 14 (1) 1997. [8] O.C. Zienkiewicz, J.Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, Int. J. Numer. Meth. Eng. 24 (1987) 337—357. [9] A. Sinha, C.V. Ramakrishnan, Growth methods and binary space partitioning for efficient finite element mesh generation, Internal Report, Applied Mechanics Dept., IIT, Delhi (1998). [10] S. Ahmad, B.M. Irons, O.C. Zienkiewicz, Analysis of thick and thin shell structures by curved finite elements, Int. J. Num. Meth. Eng. 2 (1970) 419—451. [11] E. Hinton, J.S. Campbell, Local and global smoothing of discontinuous finite element functions using a least squares method, Int. J. Num. Meth. Eng. 8 (1974) 461—480. [12] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique, Int. J. Num. Meth. Eng. 33 (1992) 1331—1364. [13] J.Z. Zhu, O.C. Zienkiewicz, E. Hinton, J. Wu, A new approach to the development of automatic quadrilateral mesh generation, Int. J. Num. Meth. Eng. 32 (1991) 849—866. [14] C.K. Lee, S.H. Lo, A new scheme for the generation of a graded quadrilateral mesh, Comput. Struct. 52 (5) (1994) 847—857. [15] R.J. Cass, S.E. Benzley, R.J. Meyers, T.D. Blacker, Generalized 3-D paving: An automated quadrilateral surface mesh generation algorithm, Int. J. Num. Meth. Eng. 39 (1996) 1475—1489. [16] S.H. Lo, Automatic mesh generation and adaptation by using contours, Int. J. Num. Meth. Eng. 31 (1991) 689—707. [17] J. Peraire, M. Vahdati, K. Morgan, O.C. Zienkiewicz, Adaptive remeshing for compressible flow computations, J. Comp. Phys. 72 (2) (1987) 449—466.
Adaptive finite element analysis of structures under transient dynamic loading using modal superposition A. Dutta , C.V. Ramakrishnan *, P. Mahajan Regional Engineering College, Silchar 788010, India Applied Mechanics Department, Indian Institute of Technology, Delhi 110016, India
Abstract In this paper, the error estimation and adaptive strategy developed for the linear elastodynamic problem under transient dynamic loading based on the Z—Z criterion is utilized for 2D and plate bending problems. An automatic mesh generator based on “growth meshing” is utilized effectively for adaptive mesh refinement. Optimal meshes are obtained iteratively corresponding to the prescribed domain discretization error limit and for a chosen number of basis modes satisfying modal truncation errors. Numerous examples show the effectiveness of the integrated approach in achieving the target accuracy in finite element transient dynamic analysis. 1999 Elsevier Science B.V. All rights reserved. Keywords: Finite element; Adaptive; Error estimation; Mesh generation; Stress recovery
1. Introduction Estimates of errors and adaptive analyses of practical finite element problems are subjects of great importance if confidence in results is needed for engineering applications. Without the assessment of the reliability of the results it is hardly reasonable to use numerical methods like the finite element method in such safety-sensitive areas as shape optimization of machine parts, construction of aeroplanes or dimensioning of nuclear power plants. Therefore, the a posteriori error estimation is now considered to be nearly as important as the finite element analysis itself. Using finite elements, the linear dynamic transient response is usually solved either by mode superposition or by direct integration schemes. The direct integration method is very highly useful for solving nonlinear problems. Procedures describing adaptive time stepping are available in [1,2] and also spatial mesh adaptation [3]. Wiberg and Li [4] have proposed a postprocessed type of a posteriori estimates in space and also in time when direct integration is used for dynamic
* Corresponding author. 0168-874X/99/$ — see front matter 1999 Elsevier Science B.V. All rights reserved PII: S 0 1 6 8 - 8 7 4 X ( 9 8 ) 0 0 0 6 1 - 4
256
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
response evaluation. It updates the spatial mesh and time step so that the discretization errors are controlled. This process is not only time consuming but also gives rise to new errors as nodal values need to be interpolated from the previous mesh to the newly generated mesh whenever a mesh is changed. But for large structural problems involving a large time domain, the mode superposition is widely used. The accuracy obtained by the mode superposition is usually dominated by the accuracy of the orthogonal modes being used and proper representation of spatially distributed loads by the number of modes selected for the modal analysis. Thus, we find that the following errors are introduced. These are (a) domain discretization error, (b) modal truncation error, (c) numerical error in eigen modes and eigen values, and (d) truncation error in numerical integration. While a proper value of time step eliminates truncation error in numerical integration to a great extent and an appropriate tolerance limit evaluates eigenmodes and eigenvalues accurately, domain discretization error and modal truncation error need elaborate study. The work reported so far on discretization error estimation under dynamic loading and using the modal superposition scheme are very few. Wilson and Joo [5] have arrived at the final mesh using Ritz vector as the basis of transformation. In their investigation, the authors have made use of modal participation and amplification factors and obtained error estimates based on Babuska’s criterion using amplified Ritz modes. This procedure is a bit cumbersome and does not give a measure of the actual error in the transient response. Cook and Avrashi [6] have discussed the procedure for estimating the discretization error of the finite element modelling as applied to the calculation of natural frequency of vibrations. Meshes are obtained corresponding to each mode. Dutta and Ramakrishnan [7] proposed a measure for discretization which is a logical extension of Zienkiewicz and Zhu’s [8] error criterion by involving time integration to consider the variation of response with time. Using this error measure, an adaptive mesh refinement strategy is proposed which yields good control over the discretization errors in transient dynamic analysis. An optimal mesh is achieved iteratively wherein meshes are refined on the basis of error indicators so that the discretization errors are within the prescribed limit for a chosen number of basis modes satisfying the modal cut-off criteria, when the modal superposition method is used. In this paper, we have adopted the adaptive strategy given by Dutta and Ramakrishnan [7] for arriving at an optimal mesh for a prescribed domain discretization error limit and for a chosen number of basis modes satisfying modal cut-off criteria for accurate evaluation of response under transient dynamic loading. A 2-D quadrilateral automatic mesh generator [9] is used for carrying out adaptive refinement of the finite element mesh. The adaptive strategy is extended to estimation of errors and adaptive analysis for plates and shells where the well-known degenerate shell element by Ahmad et al. [10] is ued for evaluation of response under dynamic loads.
2. Error measure for domain discretization The dynamic response of a linear structural dynamic system that is subjected to general dynamic loads is given Mz(#Kz#Cz"q(t). The initial conditions are z(0)"z and z(0)"z .
(1)
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
257
Here z is the vector of nodal displacements, M is assembled consistent mass matrix" L Me where n is the number of elements, C L K is assembled stiffness matrix" Ke, C C is assembled damping matrix" L Ce, (2) C q(t) is external excitation (time dependent). FE stresses are calculated as p"DBz,
(3)
where D is the elasticity matrix and B is the strain displacement matrix. The approximate solution p containing discretization errors differs from the accurate solution pH and the difference is the error e . N Thus, e "pH!p, (4) N Zienkiewicz and Zhu [8] suggested a simple approach to get better approximation of p by using a projection tecnique or nodal averaging procedure. In the projection technique, the use of weighted residual requirement for equality between pH and p leads to the computation of an accurate solution of pH which is identical to the least-squares smoothing technique of Hinton and Campbell [11]. Nodal averaging, though very simple, gives very good convergence for most of the problems. A recent and elegant technique of better stress recovery is done using superconvergent patch recovery method by Zienkiewicz and Zhu [12]. In the recovery process, it is assumed that the accurate nodal values pN H belong to a polynomial expansion pH of the same complete order ‘p’ as C that present in the basis function N and which is valid over an element patch surrounding the particular assembly node considered. Such a ‘patch’ represents a union of elements containing this vertex node (Fig. 1). This polynomial expansion will be used for each component of pH and C we get pH"Pa, (5) C where P contains the appropriate polynomial terms and a is a set of unknown parameters. For a general plate and shell problem, the polynomial can be expressed as a function of x, y and z. However, since only plate problems are solved in this paper, x and y are considered only as the parameters in the polynomial expression. Thus, it can be written as P"[1, x, y, x,2]
(6)
and a"[a , a , a , a ,2]2. (7) The determination of the unknown parameters a of the expression given in Eq. (5) is best made by ensuring a least-squares fit of this to the set of superconvergent points existing in the patch
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
258
Fig. 1. Computation of superconvergent nodal values: (䉭) Gauss point, (䢇) nodal values determined by recovery procedure, (䉺) patch assembly point.
considered. To do this we minimize L F(a)" (p (x , y )!pH(x , y )) C G G C G G G L " (p (x , y )!P(x , y )a), (8) C G G G G G where (x , y ) are the coordinates of a group of sampling points, n"mk is the total number of G G sampling points and k is the number of the sampling points on each element m (m "1,2,2,m) of H H the element patch. The minimization condition of F(a) implies that a satisfies L L P2(x , y )P(x , y )a" P2(x , y )p (x , y ). G G G G G G C G G G G This can be solved in matrix form as a"A\b,
(9)
(10)
where L L A" P2(x , y )P(x , y ) and b" P2(x , y )p (x , y ). (11) G G G G G G C G G G G Once the parameters a are determined, the recovered nodal values of pN H are simply calculated by insertion of appropriate coordinates into the expression for pH. The stresses in the nodes inside the C patch are evaluated using a from Eq. (5). The stresses at boundary nodes can be determined using interior patches.
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
259
For the time-dependent problem, we propose #e # ;100%, g" #u # where g is the overall discretization error with 2#e# dt . #e# " ¹
(12)
(13)
The integration is carried out numerically using Simpson’s rule. Also 2#u# dt #u# " ¹
(14)
where ¹ is the duration of response, #e# is the ‘energy norm’ of the error for the whole domain
L , #e#" #e# G G where
(15)
#e# " G
e2D\e dX N X N
,
(16)
the energy norm of the error in the element i, #u#"W #uM ###e# X ,
(17)
where
L #uM #" #uM # G G and
#uM # " G
X
p2D\p dX
,
(18)
the energy norm of the FE solution. The complete procedure of discretization error measure is given below in an algorithmic form: Table 1
3. Estimation of error due to modal truncation Modal truncation errors are introduced in the response calculation because of the introduction of a reduced system. Since for a specified number of modes, discretization errors have now been quantified and an appropriate mesh satisfying the prescribed tolerance limit of gN can be properly designed, the measure for control of error due to modal truncation can be carried out effectively. Inclusion of higher modes alone without taking care of domain discretization error leads to
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
260
Table 1 Algorithm for domain discretization error (1) Solve MVX"KV where V is a matrix of eigen vectors and X is matrix of eigenvalues (2) Calculate dynamic response for t"0,¹ +(2.1) For i"1, m nos. of eigenvector solve 1 + y (t)" R f (q)e\DGSGR\O sin u (t!q) dq G u G z t " y (t), G G where u is the ith damped natural frequency, f is the damping ratio corresponding to mode G G i and f is the load vector corresponding to mode i and can be written as f "2q(t). G G G (2.2) for e"1, n elements + Retrieve ze t from z t where z is the displacement vector corresponding to the element level CR for j"1, number of Gaussian point + calculate pH "DBHze t CR calculate smooth stress pHH, CR , for e"1, n elements + for j"1, number of Gaussian point + eHe "pHH!pH N R CR CR #e# " X(eHe )2D\ eH dX CR N R e NeR
#uM # " X(pHe )2D\ pH dX CR N R e NeR , for the whole struture #e# "[ #e# ] R CR
#uM # "[ #uM # ] , R CR , (3) Overall domain discretization error 2#e# dt #e# " R ¹ 2#u# dt R #u# " ¹ #u# "W #uM # ##e#X R R R #e# ;100% g" #u# (4) Discretization error at element level For i"1, n elements #e# + m " GR is a refinement indicator at element level G eN K
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
261
Table 1 Continued where
eN "gN K
#uM # ##e# R R n
is the allowable element error g is the permissible overall discretization error 2m dt m " G G ¹ If g)gN and If m (1 for i"1, n STOP G h Else h " G G (m ) N G h is the size of element i of the previous mesh, G h is the size of element i of the new mesh G Endif, Go to step 1
erroneous results [7], since there is no guarantee that the chosen mesh is suitable for the inclusion of higher modes. Given a finite element mesh, it is possible to control the modal truncation errors in finite element analysis using the ratios of MPF, modal participation factor and its maximum value (MPF ) and MSF, modal significance factor and its maximum value (MSF ) as demonstrated
by Dutta and Ramakrishnan [7]. The modal participation factor for the ith mode can be written as 2p/u. The ratio of MPF and MPF is b and b (d where d is the permissible limit. The G
modal significance factor for the ith mode can be written as 1 uN 2p ; where b " G u u ((1!b)#2f b ) G G G G G and u is the exciting frequency and the ratio of MSF and MSF is b (d.
4. Automated meshing algorithm A general formulation developed by Sinha and Ramakrishnan [9], termed ‘growth meshing’ which is characterized by rule based local selection of elements resulting in an inward propagation of the mesh and includes several existing methods of unstructured gridding, is utilized for adaptive refinement of the domain. The 2-D quadrilateral meshes are generated mainly following the strategy adopted by Zhu et al. [13], where two triangles which share a common side are combined. The main feature of the quadrilateral generation is that first two triangles which have a common side are generated and then the triangles are combined to form a quadrilateral. The process
262
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
continues until the domain is fully covered by non-overlapping quadrilaterals. Furthermore, it is rule driven and offers flexibility through the incorporated constraints along with the possibility of adaptive modification. It is very effective in achieving a rapid gradation effect over irregular and complicated domains. However, one drawback of this type of quadrilateral generation by combining triangles is that the quality of the raw quadrilateral mesh produced is usually not good and mesh quality enhancement procedures are adopted following Zhu et al. [13] and Lee and Lo [14]. In order to capture the rapid variation of element size which commonly occurs in adaptive refinement meshes, the mesh generator employs the previous quadrilateral mesh as the background mesh. Nodal spacing values of the previous mesh obtained through error estimation is utilized for calculating node spacing at a point of the new mesh to be generated. The nodal spacing at a point is calculated as [15] L (1/r)d (19) d" G G G, L 1/r G G where r is the distance between the node to be generated and the nodal point of the previous mesh, G d is the nodal spacing at the nodal point of the previous mesh and ‘n’ is the total number of nodes. G The nodal points of the previous mesh which are neighbouring close to the point under consideration influence in arriving at the value of the nodal spacing at that point. For the generation of an entire quadrilateral element mesh over the domain, the boundary of the domain is transformed into a polygon with an even number of sides. A set of uniformly distributed sample points is first placed along the boundary line and node spacing on each sample point is calculated by interpolating its value obtained from the background mesh. The boundary nodes are then generated following a simple formula given by Lo [16] which allows the position of the new nodes to be calculated one by one along the boundary. The interior nodes are generated following the strategy given by Peraire et al. [17]. Starting from an active front chosen from the boundary segments, the nodal spacing values at the background mesh are utilized such that the elements with acceptable aspect ratio and minimum distortion are generated.
5. Computation of optimal mesh using modal superposition An overall strategy for accurate computation of dynamic response is given by Dutta and Ramakrishnan [7] for plane stress/plane strain problems. Following the same strategy, the number of basis modes are first obtained for a coarse mesh to begin with, which will satisfy the cut-off criteria based on b Od or b Od adopted. Discretization errors at element level m (i"1, number G of elements) are determined for those number of modes. Average nodal spacing values are calculated using m values and mesh is then adaptively refined based on those nodal spacing values. G The number of modes required for the refined mesh are determined satisfying the modal cut-off criteria. The overall domain discretization error g and the discretization error at element level m are G determined for chosen number of modes. Iteration is carried out until the mesh is an optimal one satisfying the prescribed domain discretization error limit gN corresponding to the number of modes chosen on b or b limit.
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
263
5.1. Numerical example A number of plane stress and plate bending problems have been solved to demonstrate the adaptive strategy and use of the automatic mesh generator in obtaining the optimal meshes. Example 1 (Plane stress plate with a hole). A square plane stress plate with a central hole subjected to harmonic loading of exciting frequency uN "5 rad/s is considered. Because of symmetry, only a quadrant could be considered for analysis and Fig. 2a is assumed as the starting mesh. The limits of MSF ratios and MPF ratios are chosen as 5e-2. The limit of domain discretization error gN is considered as 1%. The geometric and material properties of the plate are as below. Plate: 100 in;100 in;1 in with a central hole of radius of 5, E"0.3e5 lbf/in, l"0.3; o"1.0 lbm/in. For the problem with uN "5 rad/s, 15 eigenmodes are selected to satisfy the modal cut-off criteria corresponding to the initial mesh (Fig. 2a). The final mesh is obtained after one iteration. The automatic mesh generator described earlier is used for adaptive refinement. The exciting frequency is very low for this example and lies between first and second natural frequencies. The results are presented in Table 2. The final mesh is shown in Fig. 2b. NIGEN is the number of eigenmodes. A similar exercise is repeated with exciting frequency uN "25 rad/s and Fig. 3a is assumed as the starting mesh. The problem with higher exciting frequency uN "25 rad/s requires 19 basis modes corresponding to the initial mesh (Fig. 3a) in order to satisfy the modal cut-off requirement of b and b O5.e-2.
Fig. 2. Plane stress plate subjected to harmonic loading (u"5): (a) Mesh 1, 12 elements, g"5.38%; (b) Mesh 2, 127 elements, g"1.12%.
264
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
Table 2 Optimal mesh for the problem of a plate with a hole subjected to harmonic loading of uN "5 rad/s Mesh
Limit of b &b
NIGEN
g (%)
Fig. 2a Fig. 2a
5.e-2 5.e-2
15 17
5.383 1.122
Fig. 3. Plane stress plate subjected to harmonic loading (u"25): (a) Mesh 1, 12 elements, g"17.84%; (b) Mesh 2, 348 elements, g"2.04%; (c) Mesh 3, 393 elements, g"1.07%.
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
265
Table 3 Optimal mesh for the problem of a plate with a hole subjected to harmonic loading of uN "25 rad/s Mesh
Limit of b &b
NIGEN
g (%)
Fig. 3a Fig. 3b Fig. 3c
5.e-2 5.e-2 5.e-2
19 23 25
17.842 2.037 1.075
The final mesh is obtained after two iterations and results are shown in Table 3. The exciting frequency lies between the seventh and eighth natural frequency. The final mesh (Fig. 3c) corresponding to this problem is quite different from the one (Fig. 3b) corresponding to uN "5 rad/s because of the integration of different numbers of basis modes in the analysis. Example 2 (Plate problems). A number of plate examples are chosen with different boundary conditions, loading pattern and geometry of the plates so that the error estimation and adaptive analysis can reflect the variation in mesh gradation. A plate clampled on all the four edges and subjected to suddenly applied pressure load is considered. A quadrant is considered because of symmetry. The initial mesh (Fig. 4a) has nine elements. Modal truncation error is controlled by keeping the limit of b and b as 1.e-3 and the overall domain discretization error tolerance gN is kept as 4%. The geometric and material properties of the plate are as below. Plate: 300 in;300 in;3 in; E"3.e7 lbf/in; l"0.316; o"7.324 e-4 lbm/in. The optimal mesh (Fig. 4c) for the clampled plate subjected to suddenly applied pressure load is obtained iteratively following the strategy described in Section 3.8 and the results are shown in Table 4. Two iterations are required to arrive at the mesh for which g is less than the prescribed tolerance limit for overall domain discretization errors and the number of modes selected for the analysis satisfy the modal cut-off criteria. Intermediate mesh and the number of basis modes required at all the stages are mentioned in Table 4. The total number of elements corresponding to the final mesh (Fig. 4c) are 234 and the mesh clearly shows the refinement requirement near the clamped boundaries where the stress gradient is high. Similarly, the clamped plate is loaded with a suddenly applied point load at the centre of the plate. The discretization pattern for such a problem under static loading is quite clear; but a study of this problem under transient dynamic loading is made to understand as to how the mesh gradation vary. The optimal mesh is determined using the prescribed limit as mentioned above. Initial mesh used has nine elements as shown in Fig. 5a. The optimal mesh is determined using the prescribed limit as mentioned above. The optimal mesh (Fig. 5c) for the clamped plate subjected to a suddenly applied point load at the centre of the plate is also obtained iteratively and the results are shown in Table 5. Table 5 Two iterations are required to arrive at the optimal mesh satisfying prescribed error limits. The intermediate mesh (Fig. 5b) and the final mesh (Fig. 5c) show that the refinement requirement is
266
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
Fig. 4. Symmetric quadrant of clamped square plate (a/t"100) with uniformly distributed suddenly applied load: (a) Mesh 1, 9 elements, g"31.62%; (b) Mesh 2, 89 elements, g"6.65%; (c) Mesh 3, 234 elements, g"3.45%.
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
267
Table 4 Optimal mesh for a clamped plate subjected to a suddenly applied pressure load Mesh
Limit of b
Fig. 4a Fig. 4b Fig. 4c
1.e-3 1.e-3 1.e-3
&b
NIGEN
No. of elements
g (%)
7 9 10
9 89 234
31.62 6.65 3.45
more near the point of application of the load which is quite obvious. The effective stresses are calculated along the centroidal line parallel to the ½-axis for all the three meshes and these stresses are plotted as shown in Fig. 6 at some specified time steps. The stress plots clearly show that the stress distribution is affected by the finite element mesh refinement and stresses stabilize to the correct value when domain discretization errors g corresponding to a finite element mesh is low. Table 5 shows that the Mesh 1 has quite a high value of g when compared to Meshes 2 and 3. Thus, the differences in stresses corresponding to Meshes 2 and 3 are much lower when compared to that of Mesh 1 except near the centre of the plate. The refinement requirement near the centre of the plate will be more stringent as the concentrated load is acting at that point and hence the stress picture near that zone can be further improved accordingly. The study shows that the adaptive mesh refinement can yield a better stress distribution over the entire problem domain and the process is fully automatic. A simply supported skew plate subjected to a suddenly applied pressure load is also considered. Material properties are the same as mentioned above. The overall domain discretization error tolerance gN is kept as 9% and b O1.e-3. The analysis is carried out for the full plate of size 300;300;3 in 3. The initial mesh has nine elements and the mesh is shown in Fig. 7a. This problem is considered since it is expected that stress concentrations will occur near the corners where the edges meet at an obtuse angle and hence such a problem is very interesting from the point of view of adaptive analysis. The optimal mesh (Fig. 7c) for the simply supported skew plate subjected to a suddenly applied pressure load is obtained following the same strategy adopted above and the results are shown in Table 6. For proper representation of the boundary conditions along the skew edges, the following considerations are made. Simply supported boundary conditions along an edge are maintained by allowing only the rotation about that edge and restraining all other degrees of freedom. At a node ‘i’ (Fig. 8) along the skew edge, the rotations about normal to that skew edge can be written as b cos h!b sin h. Using the penalty approach, this amounts to adding a term to the potential energy which can be written as C(b cos h!b sin h), where C is a large number.
268
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
Fig. 5. Symmetric quadrant of clamped square plate (a/t"100) with suddenly applied point load at the centre of the plate: (a) Mesh 1, 9 elements, g"32.37%; (b) Mesh 2, 167 elements, g"5.01%; (c) Mesh 3, 313 elements, g"3.21%.
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
269
Table 5 Optimal mesh for a clamped plate subjected to a suddenly applied point load at the centre of the plate Mesh
Limit of b
Fig. 5a Fig. 5b Fig. 5c
1.e-3 1.e-3 1.e-3
&b
NIGEN
No. of elements
g (%)
10 12 13
9 167 313
32.37 5.01 3.21
Fig. 6. Effective stress (lbf/in) vs. distance along the centroidal line: (a) At 1/8th cycle; (b) at 1/4th cycle; (c) at 3/8th cycle; (d) at 1/2th cycle.
270
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
Fig. 7. Simply supported skew plate (a/t"100) with uniformly distributed suddenly applied load: (a) Mesh 1, 9 elements, g"39.57%; (b) Mesh 2, 175 elements, g"13.69%; (c) Mesh 3, 317 elements, g"8.82%.
Table 6 Optimal mesh for a simply supported skewed plate subjected to a suddenly applied pressure load Mesh
Limit of b
Fig. 7a Fig. 7b Fig. 7c
1.e-3 1.e-3 1.e-3
&b
NIGEN
No. of elements
g (%)
8 12 15
9 175 317
39.59 13.69 8.82
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
271
Fig. 8. Node ‘i’ along skew edge.
This can be rearranged as
C cos h 1 [b b ] 2 !C sin h cos h
!C sin h cos h
b
C sin h
b
. The terms C cos h, !C sin h cos h and c sin h get added to the stiffness matrix, for every node on the incline and the new stiffness matrix is solved for the unknown variables. Similar considerations are made for translational degrees of freedom. Two iterations are required to arrive at the final mesh. The intermediate mesh (Fig. 7b) and the final mesh (Fig. 7c) show that the refinement requirement is more near the corners where the two adjoining edges meet at an obtuse angle.
6. Conclusion The error estimation scheme for plane stress/plane strain and general plate and shell problems subjected to transient dynamic loads is computationally efficient and easily implementable. Local error indicators combined with adaptive meshing is very useful, leading to solutions of specified accuracy. The overall strategy is very effective in achieving high level of accuracy introducing only the required additional degrees of freedom during mesh refinement.
References [1] O.C. Zienkiewicz, Y.M. Xie, A simple error estimator and adaptive time stepping procedure for dynamic analysis, Earthquake Eng. Struct. Dyn. 20 (1991) 871—887. [2] L.F. Zeng, N.-E. Wiberg, X.D. Li, Y.M. Xie, A posteriori local error estimation and adaptive time-stepping for Newmark Integration in dynamic analysis, Earthquake Eng. Struct. Dyn. 21 (1992) 555—571. [3] L.F. Zeng, N.-E. Wiberg, Spatial mesh adaptation in semidiscrete finite element analysis of linear elastodynamic problems, Comput. Mech. 9 (1992) 315—332. [4] N.E. Wiberg, X.D. Li, A post-processed error estimate and an adaptive procedure for the semidiscrete finite element method in dynamic analysis, Int. J. Numer. Meth. Eng. 37 (1994) 3585—3603. [5] K.J. Joo, E.L. Wilson, An adaptive finite element technique for structural dynamic analysis, Comput. Struct. 30 (6) (1988) 1319—1339.
272
A. Dutta et al. / Finite Elements in Analysis and Design 31 (1999) 255—272
[6] R.D. Cook, J. Avrashi, Error estimation and adaptive meshing for vibration problems, Comput. Struct. 44 (3) (1992) 619—626. [7] A. Dutta, C.V. Ramakrishnan, Error estimation in finite element transient dynamic analysis using modal superposition method, Eng. Comput. 14 (1) 1997. [8] O.C. Zienkiewicz, J.Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, Int. J. Numer. Meth. Eng. 24 (1987) 337—357. [9] A. Sinha, C.V. Ramakrishnan, Growth methods and binary space partitioning for efficient finite element mesh generation, Internal Report, Applied Mechanics Dept., IIT, Delhi (1998). [10] S. Ahmad, B.M. Irons, O.C. Zienkiewicz, Analysis of thick and thin shell structures by curved finite elements, Int. J. Num. Meth. Eng. 2 (1970) 419—451. [11] E. Hinton, J.S. Campbell, Local and global smoothing of discontinuous finite element functions using a least squares method, Int. J. Num. Meth. Eng. 8 (1974) 461—480. [12] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique, Int. J. Num. Meth. Eng. 33 (1992) 1331—1364. [13] J.Z. Zhu, O.C. Zienkiewicz, E. Hinton, J. Wu, A new approach to the development of automatic quadrilateral mesh generation, Int. J. Num. Meth. Eng. 32 (1991) 849—866. [14] C.K. Lee, S.H. Lo, A new scheme for the generation of a graded quadrilateral mesh, Comput. Struct. 52 (5) (1994) 847—857. [15] R.J. Cass, S.E. Benzley, R.J. Meyers, T.D. Blacker, Generalized 3-D paving: An automated quadrilateral surface mesh generation algorithm, Int. J. Num. Meth. Eng. 39 (1996) 1475—1489. [16] S.H. Lo, Automatic mesh generation and adaptation by using contours, Int. J. Num. Meth. Eng. 31 (1991) 689—707. [17] J. Peraire, M. Vahdati, K. Morgan, O.C. Zienkiewicz, Adaptive remeshing for compressible flow computations, J. Comp. Phys. 72 (2) (1987) 449—466.