- Email: [email protected]
Adaptive finite element analysis of structures subjected to transient dynamic loads using time marching scheme
Computers and Structures 80 (2002) 2313–2319 www.elsevier.com/locate/compstruc
Adaptive finite element analysis of structures subjected to transient dynamic loads using time marching scheme Anjan Dutta Department of Civil Engineering, Indian Institute Technology, Guwahati 781039, India Received 14 November 2000; accepted 4 July 2002
Abstract Presents a method of obtaining accurate response using the finite element analysis for linear elastodynamic problems under transient dynamic loading. The accuracy obtained would depend on the discretization of both space and time. An integrated method is presented where both space and time discretization errors are evaluated and iteratively converges to a solution of desired accuracy through adaptive refinement. In this paper, it is proposed to solve the response using Newmark-b method. A numbers of examples are solved to demonstrate the effectiveness of the procedure. Ó 2002 Civil-Comp Ltd. and Elsevier Science Ltd. All rights reserved. Keywords: Finite element; Transient dynamic; Discretization; Error; Refinement adaptivity
1. Introduction Estimation of errors and adaptive analysis of practical finite element problems are subject of great importance if confidence in results is needed for engineering application. 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. Babuska and Rheinboldt [1,2] have estimated the errors in finite element solution using the residuals in the equilibrium equations created by replacing the actual solution by finite element solution. These error estimates, although very accurate, are difficult to implement computationally. Ladeveze and Leguillon [3] and Ladeveze et al. [4] have proposed an error estimate based on the concept of error in the constitutive relation. Zienkiewicz and Zhu [5] have estimated the error based
E-mail address: [email protected] (A. Dutta).
on a smoothing of the finite element stresses. This is very popular method as this can be easily implemented in the existing finite element code. Using finite elements, linear dynamic transient response is usually solved either by mode superposition or by direct integration schemes. Procedures describing adaptive time stepping are available in Refs. [6,7] and also spatial mesh adaptation [8]. Wiberg and Li [9] have proposed a post-processed type of a posteriori estimates in space and also in time when direct integration is used for dynamic 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. Combe et al. [10] have proposed an error estimate based on the concept of error in the constitutive relation for transient dynamic loading. Schleupen and Ramm [11] have proposed local and global error estimations for linear structural dynamic problems. Global method is based on post-processing techniques starting from a semi-discrete formulation. Temporal error due to finite difference discretization and spatial error due to finite element discretization are calculated independently. The temporal error estimators are applied within one time step and the spatial error estimators at a time
0045-7949/02/$ - see front matter Ó 2002 Civil-Comp Ltd. and Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 2 ) 0 0 2 6 6 - 3
2314
A. Dutta / Computers and Structures 80 (2002) 2313–2319
point. The local error estimators are designed to evaluate the error of local variables in a certain region by applying duality techniques. 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. The work reported so far on discretization error estimation under dynamic loading and using modal superposition scheme are very few. Joo and Wilson [12] have arrived at the final mesh using Ritz vector as basis of transformation. In their investigation, the authors have made use of modal participation and amplification factor and obtained error estimates based on Babuska’s criterion using amplified Ritz modes. Cook and Avrashi [13] have discussed the procedure for estimating the discretization error of the finite element modeling as applied to the calculation of natural frequency of vibrations. Meshes are obtained corresponding to each mode. Dutta and Ramakrishnan [14] proposed a measure for discretization, which is a logical extension of Zienkiewicz and Zhu [5] error criterion by involving time integration to consider the variation of response with time. In this paper, we have adopted the adaptive strategy given by Dutta and Ramakrishnan [14] for arriving at an optimal mesh for a prescribed domain discretization error limit and integrated with an adaptive time stepping procedure proposed by Zienkiewicz and Xie [6] with some modification. The adaptive strategy is extended to estimation of errors and adaptive analysis for plates and shells, where the well-known degenerated shell element by Ahmad et al. [15] is used for evaluation of response under dynamic loads.
Such a ‘patch’ represents a union of elements containing this vertex node. This polynomial expansion will be used for each component of re and we get re ¼ Pa
ð3Þ
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; x2 ; . . .
ð4Þ
and a ¼ ½a1 ; a2 ; a3 ; a4 ; . . .T
ð5Þ
The determination of the unknown parameters a of the expression given in Eq. (3) is best made by ensuring a least square fit of this to the set of superconvergent points existing in the patch considered. To do this we minimize F ðaÞ ¼
n X
re ðxi ; yi Þ re ðxi ; yi Þ
i¼1
¼
n X
ðre ðxi ; yi Þ P ðxi ; yi ÞaÞ2
where ðxi ; yi Þ are the co-ordinates of a group of sampling points, n ¼ mk is the total number sampling points and k is the number of the sampling points on each element mj ðmj ¼ 1; 2; . . . ; mÞ of the element patch. The minimization condition of F ðaÞ implies that a satisfies n X
P T ðxi ; yi ÞP ðxi ; yi Þa ¼
n X
P T ðxi ; yi Þre ðxi ; yi Þ
ð7Þ
i¼1
This can be solved in matrix form as a ¼ A1 b
Finite element stresses are calculated as r ¼ DBz
ð6Þ
i¼1
i¼1
2. Error measure for discretization
2
ð1Þ
The approximate solution r containing discretization errors differs from the accurate solution r and the difference is the error er
where A ¼
ð8Þ n X
P T ðxi ; yi ÞP ðxi ; yi Þ and
i¼1
b¼
n X
P T ðxi yi Þre ðxi ; yi Þ
ð9Þ
i¼1
thus; er ¼ r r
ð2Þ
A recent and elegant technique for the evaluation of accurate stress is done using superconvergent patch recovery method by Zienkiewicz and Zhu [16]. In the recovery process, it is assumed that the accurate nodal values r belong to a polynomial expansion re of the same complete order p as that present in the basis function N and which is valid over an element patch surrounding the particular assembly node considered.
Once the parameter a are determined the recovered nodal vales of r are simply calculated by insertion of appropriate co-ordinate into the expression for re . The stresses in the nodes inside the patch are evaluated using a from Eq. (3). The stresses at boundary nodes can be determined using interior patches. Zienkiewicz and Boroomand [17] demonstrated that the accurate stress recovery is possible even without implementing the superconvergent stresses for some particular examples.
A. Dutta / Computers and Structures 80 (2002) 2313–2319
The mathematical basis and broad details of the algorithm are described in Ref. [14]. The complete procedure of discretization error measure is given below in an algorithmic form. 2.1. Algorithm for domain discretization error (1) Solve M€z þ Kz ¼ qðtÞ (1.1) Using Newmark-b method calculate dynamic response for te½0; T 1 2 zðt þ DtÞ ¼ zðtÞ þ z_ ðtÞDt þ Dt b €zðtÞ 2 þ b€zðt þ DtÞ z_ ðt þ DtÞ ¼ z_ ðtÞ þ Dtfð1 cÞ€zðtÞ þ c€zðt þ DtÞg Rearranging the first equation for displacement at ðt þ DtÞ, we can get 1 1 z_ ðtÞ ðzðt þ DtÞ zðtÞÞ €zðt þ DtÞ ¼ bDt2 bDt 1 1 €zðtÞ 2b where b ¼ 1=4 and c ¼ 1=2 (1.2) for e ¼ 1, n elements { Retrieve zeðtÞ from zðtÞ for j ¼ 1, number of Gaussian point ðjÞ { Calculate reðtÞ ¼ DBðjÞ zeðtÞ where D and B are usual elasticity and strain displacement matrices respectively. ðjÞ Calculate smooth stress reðtÞ ðjÞ
reðtÞ ¼ DBðjÞ zeðtÞ Energy norm for error e at element level Z 1=2
T ðjÞ ðjÞ ereðtÞ D1 e dX kekeðtÞ ¼ e reðtÞ X
Energy norm for FEM solution at element level Z 1=2
T ðjÞ ðjÞ kukeðtÞ ¼ reðtÞ D1 r dX e eðtÞ
2315
(2) Overall domain discretization error RT kekðtÞ dt kekav ¼ 0 T RT 0
kukav ¼
kukðtÞ dt T
h i1=2 uk2ðtÞ þ kek2ðtÞ kukðtÞ ¼ k Domain discretization error g¼
kekav 100% kukav
(3) Discretization error at element level For i ¼ 1, n elements {Error at element level, ni ¼
kekiðtÞ em
where " em ¼ g
k uk2ðtÞ þ kek2ðtÞ
#1=2
n
and g is the permissible domain discretization error. Time averaged error at element level, RT ni dt ni0 ¼ 0 T If g 6 g and If ni0 < 1 for i ¼ 1, n STOP Else New element size, hinew ¼
hiold ðni0 Þ/
where / ¼ 1=p for no singularity and / ¼ 1=k for element close to singularity (p is the order of defining polynomial and k is the strength of singularity) Endif } Go to step 1
X
} for the whole structure hX i1=2 kekðtÞ ¼ kek2eðtÞ
kukðtÞ ¼ } }
X 1=2 2 u eðtÞ
3. Adaptive time stepping In order to control the time discretization error in the time marching scheme, some means of estimating the errors should be introduced so that the steps can be suitably adjusted. The adaptive time stepping is aimed at maintaining largest possible step size while keeping the accuracy within the prescribed limit. In this work, the error measure proposed by Zienkiewicz and Xie [6] is adopted with some modification. The total time domain
2316
A. Dutta / Computers and Structures 80 (2002) 2313–2319
T is divided into a finite number ðnT Þ of sub-domains. A uniform value of Dt is maintained in a sub-domain and the time discretization errors are calculated. Time integration of errors is carried out over individual subdomain interval to arrive at a new value of step size for that sub-domain. The local error in time is estimated as 1 e ¼ Dt2 b ð10Þ ð€znþ1 €zn Þ 6
every time sub-domain. Mesh is then adaptively refined based on the space discretization error values at element level and new values of Dt are obtained for every time sub-domain. Iteration is carried out until the mesh is an optimal one and time discretization errors in all the subdomains are within the prescribed limit.
where €znþ1 is acceleration vector at ðn þ 1Þth time step. €zn is acceleration vector at nth time step and b ¼ 1=4. The L2 norm of the local error is 1 k€znþ1 €zn k ð11Þ gt ¼ kek ¼ Dt2 b 6
A number of example problems subjected to transient dynamic loads have been considered for demonstration of the integrated adaptive procedure. Since access to exact solution is not available in general, computation of real error is not always possible. However, a simple case is considered first, whose analytical solution has been compared with the solution obtained using adaptive finite element module. This example will demonstrate the convergence of the finite element result obtained using integrated space and time adaptivity towards the analytical solution, which can be treated as near to exact solution. A simply supported beam 10 m span with uniform cross-section 200 mm 500 mm and subjected to a concentrated constant force of 1000 N suddenly applied at mid-span is considered. The properties are E ¼ 2:1 1011 N/m2 , m ¼ 0:3 and q ¼ 7850 kg/m3 . The analytical solution for response of the beam at mid-span with time is given as [18] 2P X 1 np np yðL=2; tÞ ¼ sin tÞ sin ð1 cos x n mL n x2n 2 2
Usually, it is inconvenient to specify an absolute error tolerance. To measure the relative error, it may be defined as gR ¼
kek ðkukÞmax
ð12Þ
where ðkukÞmax is the maximum value of the corresponding norm of the displacement solution recorded during the computation. Time integration of error over sub-domain leads to gRi ¼
kekavei ðkukÞmax
ð13Þ
where R Ti kekavei ¼
0
kekdt Ti
for i ¼ 1; nT
Ti is the time interval of sub-domain ‘i’. In this section, the error gRi in a sub-domain i is used to calculate new step size for the next iteration so that the error is roughly equal to the prescribed tolerance. With this tolerance given as gt , an upper limit ðc1 gt Þ and lower limit c2 gt is also specified, where the parameter c1 and c2 are in the ranges of 0 6 c1 6 1 and c2 P 1. When the error gRi exceeds the upper limit the step size needs to be reduced and the new step size in a subdomain i may be predicted as !1=3 gt Dtnew ¼ Dtold ð14Þ gRi Similarly, if the error gRi is smaller than lower limit, then the step size may be increased using the above expression. Hence the overall strategy can be written as below. Starting with a coarse mesh and a value of Dt for every time sub-domain, overall domain discretization error and discretization errors at element level are calculated and time discretization errors are calculated for
4. Numerical examples
where m is the mass per unit length, P is the magnitude of suddenly applied load, L is the span of the beam, n is the number of mode shapes and ffithe natural frequency pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi xn is given as xn ¼ n2 p2 EI=mL4 . The analytical solution at time t ¼ 0:0426 (i.e. T1 =2) considering 10 number of modes is 9:522 105 . Finite element analysis is carried out by considering four elements and dividing the time domain into five equal sub-domains. Each subdomain is divided into five uniform steps. The time step in each sub-domain is taken as 3:408 103 . The response at mid-span for t ¼ 0:0426 is 8:865 105 . The percentage of error in the finite element solution compared to the analytical solution is 6.9%. Assuming g and gt as 1% and 1% respectively, space and time discretization errors are determined. Refined mesh contains eight elements and updated time step in each subdomain is shown in Table 1. The response at mid-span for t ¼ 0:0426 is 9:3827 105 and the percentage of error in the finite element solution compared to the analytical solution is 1.46%. This simple example nicely illustrates the performance of the proposed algorithm. A few more examples are considered for further demonstration.
A. Dutta / Computers and Structures 80 (2002) 2313–2319 Table 1 Simply supported beam subjected to suddenly applied load
2317
Table 2 Cantilever beam subjected to suddenly applied load
Mesh
No. of time sub-domain, nT
g
Dt
Mesh
No. of time sub-domain, nT
g
Dt
Four elements Zero elements
5 5
5.92 1.09
3:408 103 1:86 103 2:21 103 2:22 103 2:2 103 1:91 103
Fig. 1b Fig. 1c
5 5
16.78 9.89
4:640 104 2:851 104 3:706 104 3:706 104 3:706 104 2:316 104
citing frequency 25 rad/s. Because of symmetry, only a quadrant could be considered for analysis. Initial mesh is shown in Fig. 2a and time domain is divided into five equal sub-domains and each sub-domain is divided into 10 uniform steps. The time step in each sub-domain is taken as 0.0461 for the mesh in Fig. 2a. Assuming g and gt as 1% and 1% respectively, space and time discretization errors are determined. Final refined mesh
Fig. 1. (a) Cantilever beam, (b) initial mesh and (c) final mesh. Cantilever beam subjected to suddenly applied load at quarter span.
A cantilever beam 20 m long and 2.5 m deep (Fig. 1a) subjected to suddenly applied load of 6 105 N at quarter span is considered with E ¼ 2:1 1011 N/m2 , m ¼ 0:3 and q ¼ 7850 kg/m3 . Initial mesh (Fig. 1b) is having only eight elements. Time domain is divided into five equal sub-domains and each sub-domain is divided into 10 uniform steps. The time step in each sub-domain is taken as 4:64 104 for the mesh in Fig. 1b. Assuming g and gt as 10% and 1% respectively, space and time discretization errors are determined. Refined mesh is shown in Fig. 1c and the updated time step in each subdomain is shown in Table 2. A plane stress plate 100 in: 100 in: 1 in. with a central hole of radius 5 in. is considered. The material properties are E ¼ 0:3 105 lbf/in.2 , m ¼ 0:3 and q ¼ 1:0 lbm/in.3 The plate is subjected to harmonic loading ex-
Fig. 2. (a) Initial mesh and (b) final mesh. Plane stress plate subjected to harmonic loading.
2318
A. Dutta / Computers and Structures 80 (2002) 2313–2319
is obtained after two iterations and is shown in Fig. 2b and the updated time step in each sub-domain is shown in Table 3. Table 3 Plane stress plate with a central hole subjected to harmonic loading Mesh
No. of time sub-domain, nT
g
Dt
Fig. 2a Fig. 2b
5 5
17.84 1.07
0.0461 0.0352 0.0426 0.0415 0.0425 0.0332
Table 4 Clamped plate with a suddenly applied point load at the center of the plate Mesh
No. of time sub-domain, nT
g
Dt
Fig. 3a Fig. 3b
5 5
32.34 3.11
3:726 103 3:242 103 3:221 103 3:247 103 3:199 103 3:319 103
A clamped plate 300 in: 300 in: 3 in. is loaded with a suddenly applied point load at the center of the plate. The material properties are E ¼ 3 107 lbf/in.2 , m ¼ 0:316 and q ¼ 7:324 104 lbm/in.3 A quadrant is considered because of symmetry. The initial mesh is shown in Fig. 3a and time domain is divided again into five equal sub-domains and each sub-domain is divided into ten uniform steps. The time step in each subdomain is taken as 3:726 103 for the mesh in Fig. 3a. Assuming g and gt as 3% and 1% respectively, space and time discretization errors are determined. Final refined mesh is obtained after two iterations and is shown in Fig. 3b and the updated time step in each sub-domain is shown in Table 4. 5. Conclusion The error estimation scheme for plane stress/plane strain and plate problems subjected to transient dynamic loads is computationally efficient and easily implementable. Local error measures (both in space and time) combined with adaptive meshing is very useful, converging iteratively to a solution of desired accuracy. The overall strategy is very useful in attaining high level of accuracy introducing only the required additional degrees of freedom during mesh refinement. However, the proposed time averaged error estimation will not perform well in the case of wave propagation problem as the peak stress changes its location with time as the wave propagates and hence the algorithm is particularly valid for a class of transient dynamic load, where the load location does not change. References
Fig. 3. (a) Initial mesh and (b) final mesh. Symmetric quadrant of clamped square plate with suddenly applied load at the center of the plate.
[1] Babuska I, Rheinboldt WC. Error estimators for adaptive finite element computation. SIAM J Numer Anal 1978; 15(4):736–54. [2] Babuska I, Rheinboldt WC. Adaptive approaches and reliability estimations in finite element analysis. Comput Meth Appl Mech Engng 1979;17/18:519–40. [3] Ladeveze P, Leguillon D. Error estimate procedure in the finite element method and application. SIAM J Numer Anal 1983;20(3):485–509.
A. Dutta / Computers and Structures 80 (2002) 2313–2319 [4] Ladeveze P, Pelle JP, Rougeot P. Error estimation and mesh optimization for classical finite elements. Engng Comput 1991;8:69–80. [5] Zienkiewicz OC, Zhu JZ. A simple error estimator and adaptive procedure for practical engineering analysis. Int J Numer Meth Engng 1987;24:337–57. [6] Zienkiewicz OC, Xie YM. A simple error estimator and adaptive time stepping procedure for dynamic analysis. Earthquake Engng Struct Dyn 1991;20:871–87. [7] Zeng LF, Wiberg NE, Li XD, Xie YM. A posteriori local error estimation and adaptive time-stepping for Newmark integration in dynamic analysis. Earthquake Engng Struct Dyn 1992;21:555–71. [8] Zeng LF, Wiberg NE. Spatial mesh adaptation in semidiscrete finite element analysis of linear ealstodynamic problems. Comput Mech 1992;9:315–32. [9] Wiberg NE, Li XD. A post-processed error estimate and an adaptive procedure for the semidiscrete finite element method in dynamic analysis. Int J Numer Meth Engng 1994;37:3585–603. [10] Combe JP, Ladeveze P, Pelle JP. Constitutive relation error estimation for transient finite element analysis. Comput Meth Appl Mech Engng 1999;176:165–85.
2319
[11] Schleupen A, Ramm E. Local and global error estimations in linear structural dynamics. Comput Struct 2000;76:741– 56. [12] Joo KJ, Wilson EL. An adaptive finite element technique for structural dynamic analysis. Comput Struct 1988; 30:1319–39. [13] Cook RD, Avrashi J. Error estimation and adaptive meshing for vibration problems. Comput Struct 1992; 44:619–26. [14] Dutta A, Ramakrishnan CV. Error estimation in finite element transient dynamic analysis using modal superposition method. Engng Comput 1997;14(1):135–58. [15] Ahmad S, Irons BM, Zienkiewicz OC. Analysis of thick and thin shell structures by curved finite elements. Int J Numer Meth Engng 1970;2:419–51. [16] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique. Int J Numer Meth Engng 1992; 33:1331–64. [17] Zienkiewicz OC, Boroomand B. Recovery by equilibrium in patches. Int J Numer Meth Engng 1997;40:137–64. [18] Paz M. Structural Dynamics. New Delhi: CBS Publishers and Distributors; 1987.
Adaptive finite element analysis of structures subjected to transient dynamic loads using time marching scheme Anjan Dutta Department of Civil Engineering, Indian Institute Technology, Guwahati 781039, India Received 14 November 2000; accepted 4 July 2002
Abstract Presents a method of obtaining accurate response using the finite element analysis for linear elastodynamic problems under transient dynamic loading. The accuracy obtained would depend on the discretization of both space and time. An integrated method is presented where both space and time discretization errors are evaluated and iteratively converges to a solution of desired accuracy through adaptive refinement. In this paper, it is proposed to solve the response using Newmark-b method. A numbers of examples are solved to demonstrate the effectiveness of the procedure. Ó 2002 Civil-Comp Ltd. and Elsevier Science Ltd. All rights reserved. Keywords: Finite element; Transient dynamic; Discretization; Error; Refinement adaptivity
1. Introduction Estimation of errors and adaptive analysis of practical finite element problems are subject of great importance if confidence in results is needed for engineering application. 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. Babuska and Rheinboldt [1,2] have estimated the errors in finite element solution using the residuals in the equilibrium equations created by replacing the actual solution by finite element solution. These error estimates, although very accurate, are difficult to implement computationally. Ladeveze and Leguillon [3] and Ladeveze et al. [4] have proposed an error estimate based on the concept of error in the constitutive relation. Zienkiewicz and Zhu [5] have estimated the error based
E-mail address: [email protected] (A. Dutta).
on a smoothing of the finite element stresses. This is very popular method as this can be easily implemented in the existing finite element code. Using finite elements, linear dynamic transient response is usually solved either by mode superposition or by direct integration schemes. Procedures describing adaptive time stepping are available in Refs. [6,7] and also spatial mesh adaptation [8]. Wiberg and Li [9] have proposed a post-processed type of a posteriori estimates in space and also in time when direct integration is used for dynamic 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. Combe et al. [10] have proposed an error estimate based on the concept of error in the constitutive relation for transient dynamic loading. Schleupen and Ramm [11] have proposed local and global error estimations for linear structural dynamic problems. Global method is based on post-processing techniques starting from a semi-discrete formulation. Temporal error due to finite difference discretization and spatial error due to finite element discretization are calculated independently. The temporal error estimators are applied within one time step and the spatial error estimators at a time
0045-7949/02/$ - see front matter Ó 2002 Civil-Comp Ltd. and Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 2 ) 0 0 2 6 6 - 3
2314
A. Dutta / Computers and Structures 80 (2002) 2313–2319
point. The local error estimators are designed to evaluate the error of local variables in a certain region by applying duality techniques. 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. The work reported so far on discretization error estimation under dynamic loading and using modal superposition scheme are very few. Joo and Wilson [12] have arrived at the final mesh using Ritz vector as basis of transformation. In their investigation, the authors have made use of modal participation and amplification factor and obtained error estimates based on Babuska’s criterion using amplified Ritz modes. Cook and Avrashi [13] have discussed the procedure for estimating the discretization error of the finite element modeling as applied to the calculation of natural frequency of vibrations. Meshes are obtained corresponding to each mode. Dutta and Ramakrishnan [14] proposed a measure for discretization, which is a logical extension of Zienkiewicz and Zhu [5] error criterion by involving time integration to consider the variation of response with time. In this paper, we have adopted the adaptive strategy given by Dutta and Ramakrishnan [14] for arriving at an optimal mesh for a prescribed domain discretization error limit and integrated with an adaptive time stepping procedure proposed by Zienkiewicz and Xie [6] with some modification. The adaptive strategy is extended to estimation of errors and adaptive analysis for plates and shells, where the well-known degenerated shell element by Ahmad et al. [15] is used for evaluation of response under dynamic loads.
Such a ‘patch’ represents a union of elements containing this vertex node. This polynomial expansion will be used for each component of re and we get re ¼ Pa
ð3Þ
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; x2 ; . . .
ð4Þ
and a ¼ ½a1 ; a2 ; a3 ; a4 ; . . .T
ð5Þ
The determination of the unknown parameters a of the expression given in Eq. (3) is best made by ensuring a least square fit of this to the set of superconvergent points existing in the patch considered. To do this we minimize F ðaÞ ¼
n X
re ðxi ; yi Þ re ðxi ; yi Þ
i¼1
¼
n X
ðre ðxi ; yi Þ P ðxi ; yi ÞaÞ2
where ðxi ; yi Þ are the co-ordinates of a group of sampling points, n ¼ mk is the total number sampling points and k is the number of the sampling points on each element mj ðmj ¼ 1; 2; . . . ; mÞ of the element patch. The minimization condition of F ðaÞ implies that a satisfies n X
P T ðxi ; yi ÞP ðxi ; yi Þa ¼
n X
P T ðxi ; yi Þre ðxi ; yi Þ
ð7Þ
i¼1
This can be solved in matrix form as a ¼ A1 b
Finite element stresses are calculated as r ¼ DBz
ð6Þ
i¼1
i¼1
2. Error measure for discretization
2
ð1Þ
The approximate solution r containing discretization errors differs from the accurate solution r and the difference is the error er
where A ¼
ð8Þ n X
P T ðxi ; yi ÞP ðxi ; yi Þ and
i¼1
b¼
n X
P T ðxi yi Þre ðxi ; yi Þ
ð9Þ
i¼1
thus; er ¼ r r
ð2Þ
A recent and elegant technique for the evaluation of accurate stress is done using superconvergent patch recovery method by Zienkiewicz and Zhu [16]. In the recovery process, it is assumed that the accurate nodal values r belong to a polynomial expansion re of the same complete order p as that present in the basis function N and which is valid over an element patch surrounding the particular assembly node considered.
Once the parameter a are determined the recovered nodal vales of r are simply calculated by insertion of appropriate co-ordinate into the expression for re . The stresses in the nodes inside the patch are evaluated using a from Eq. (3). The stresses at boundary nodes can be determined using interior patches. Zienkiewicz and Boroomand [17] demonstrated that the accurate stress recovery is possible even without implementing the superconvergent stresses for some particular examples.
A. Dutta / Computers and Structures 80 (2002) 2313–2319
The mathematical basis and broad details of the algorithm are described in Ref. [14]. The complete procedure of discretization error measure is given below in an algorithmic form. 2.1. Algorithm for domain discretization error (1) Solve M€z þ Kz ¼ qðtÞ (1.1) Using Newmark-b method calculate dynamic response for te½0; T 1 2 zðt þ DtÞ ¼ zðtÞ þ z_ ðtÞDt þ Dt b €zðtÞ 2 þ b€zðt þ DtÞ z_ ðt þ DtÞ ¼ z_ ðtÞ þ Dtfð1 cÞ€zðtÞ þ c€zðt þ DtÞg Rearranging the first equation for displacement at ðt þ DtÞ, we can get 1 1 z_ ðtÞ ðzðt þ DtÞ zðtÞÞ €zðt þ DtÞ ¼ bDt2 bDt 1 1 €zðtÞ 2b where b ¼ 1=4 and c ¼ 1=2 (1.2) for e ¼ 1, n elements { Retrieve zeðtÞ from zðtÞ for j ¼ 1, number of Gaussian point ðjÞ { Calculate reðtÞ ¼ DBðjÞ zeðtÞ where D and B are usual elasticity and strain displacement matrices respectively. ðjÞ Calculate smooth stress reðtÞ ðjÞ
reðtÞ ¼ DBðjÞ zeðtÞ Energy norm for error e at element level Z 1=2
T ðjÞ ðjÞ ereðtÞ D1 e dX kekeðtÞ ¼ e reðtÞ X
Energy norm for FEM solution at element level Z 1=2
T ðjÞ ðjÞ kukeðtÞ ¼ reðtÞ D1 r dX e eðtÞ
2315
(2) Overall domain discretization error RT kekðtÞ dt kekav ¼ 0 T RT 0
kukav ¼
kukðtÞ dt T
h i1=2 uk2ðtÞ þ kek2ðtÞ kukðtÞ ¼ k Domain discretization error g¼
kekav 100% kukav
(3) Discretization error at element level For i ¼ 1, n elements {Error at element level, ni ¼
kekiðtÞ em
where " em ¼ g
k uk2ðtÞ þ kek2ðtÞ
#1=2
n
and g is the permissible domain discretization error. Time averaged error at element level, RT ni dt ni0 ¼ 0 T If g 6 g and If ni0 < 1 for i ¼ 1, n STOP Else New element size, hinew ¼
hiold ðni0 Þ/
where / ¼ 1=p for no singularity and / ¼ 1=k for element close to singularity (p is the order of defining polynomial and k is the strength of singularity) Endif } Go to step 1
X
} for the whole structure hX i1=2 kekðtÞ ¼ kek2eðtÞ
kukðtÞ ¼ } }
X 1=2 2 u eðtÞ
3. Adaptive time stepping In order to control the time discretization error in the time marching scheme, some means of estimating the errors should be introduced so that the steps can be suitably adjusted. The adaptive time stepping is aimed at maintaining largest possible step size while keeping the accuracy within the prescribed limit. In this work, the error measure proposed by Zienkiewicz and Xie [6] is adopted with some modification. The total time domain
2316
A. Dutta / Computers and Structures 80 (2002) 2313–2319
T is divided into a finite number ðnT Þ of sub-domains. A uniform value of Dt is maintained in a sub-domain and the time discretization errors are calculated. Time integration of errors is carried out over individual subdomain interval to arrive at a new value of step size for that sub-domain. The local error in time is estimated as 1 e ¼ Dt2 b ð10Þ ð€znþ1 €zn Þ 6
every time sub-domain. Mesh is then adaptively refined based on the space discretization error values at element level and new values of Dt are obtained for every time sub-domain. Iteration is carried out until the mesh is an optimal one and time discretization errors in all the subdomains are within the prescribed limit.
where €znþ1 is acceleration vector at ðn þ 1Þth time step. €zn is acceleration vector at nth time step and b ¼ 1=4. The L2 norm of the local error is 1 k€znþ1 €zn k ð11Þ gt ¼ kek ¼ Dt2 b 6
A number of example problems subjected to transient dynamic loads have been considered for demonstration of the integrated adaptive procedure. Since access to exact solution is not available in general, computation of real error is not always possible. However, a simple case is considered first, whose analytical solution has been compared with the solution obtained using adaptive finite element module. This example will demonstrate the convergence of the finite element result obtained using integrated space and time adaptivity towards the analytical solution, which can be treated as near to exact solution. A simply supported beam 10 m span with uniform cross-section 200 mm 500 mm and subjected to a concentrated constant force of 1000 N suddenly applied at mid-span is considered. The properties are E ¼ 2:1 1011 N/m2 , m ¼ 0:3 and q ¼ 7850 kg/m3 . The analytical solution for response of the beam at mid-span with time is given as [18] 2P X 1 np np yðL=2; tÞ ¼ sin tÞ sin ð1 cos x n mL n x2n 2 2
Usually, it is inconvenient to specify an absolute error tolerance. To measure the relative error, it may be defined as gR ¼
kek ðkukÞmax
ð12Þ
where ðkukÞmax is the maximum value of the corresponding norm of the displacement solution recorded during the computation. Time integration of error over sub-domain leads to gRi ¼
kekavei ðkukÞmax
ð13Þ
where R Ti kekavei ¼
0
kekdt Ti
for i ¼ 1; nT
Ti is the time interval of sub-domain ‘i’. In this section, the error gRi in a sub-domain i is used to calculate new step size for the next iteration so that the error is roughly equal to the prescribed tolerance. With this tolerance given as gt , an upper limit ðc1 gt Þ and lower limit c2 gt is also specified, where the parameter c1 and c2 are in the ranges of 0 6 c1 6 1 and c2 P 1. When the error gRi exceeds the upper limit the step size needs to be reduced and the new step size in a subdomain i may be predicted as !1=3 gt Dtnew ¼ Dtold ð14Þ gRi Similarly, if the error gRi is smaller than lower limit, then the step size may be increased using the above expression. Hence the overall strategy can be written as below. Starting with a coarse mesh and a value of Dt for every time sub-domain, overall domain discretization error and discretization errors at element level are calculated and time discretization errors are calculated for
4. Numerical examples
where m is the mass per unit length, P is the magnitude of suddenly applied load, L is the span of the beam, n is the number of mode shapes and ffithe natural frequency pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi xn is given as xn ¼ n2 p2 EI=mL4 . The analytical solution at time t ¼ 0:0426 (i.e. T1 =2) considering 10 number of modes is 9:522 105 . Finite element analysis is carried out by considering four elements and dividing the time domain into five equal sub-domains. Each subdomain is divided into five uniform steps. The time step in each sub-domain is taken as 3:408 103 . The response at mid-span for t ¼ 0:0426 is 8:865 105 . The percentage of error in the finite element solution compared to the analytical solution is 6.9%. Assuming g and gt as 1% and 1% respectively, space and time discretization errors are determined. Refined mesh contains eight elements and updated time step in each subdomain is shown in Table 1. The response at mid-span for t ¼ 0:0426 is 9:3827 105 and the percentage of error in the finite element solution compared to the analytical solution is 1.46%. This simple example nicely illustrates the performance of the proposed algorithm. A few more examples are considered for further demonstration.
A. Dutta / Computers and Structures 80 (2002) 2313–2319 Table 1 Simply supported beam subjected to suddenly applied load
2317
Table 2 Cantilever beam subjected to suddenly applied load
Mesh
No. of time sub-domain, nT
g
Dt
Mesh
No. of time sub-domain, nT
g
Dt
Four elements Zero elements
5 5
5.92 1.09
3:408 103 1:86 103 2:21 103 2:22 103 2:2 103 1:91 103
Fig. 1b Fig. 1c
5 5
16.78 9.89
4:640 104 2:851 104 3:706 104 3:706 104 3:706 104 2:316 104
citing frequency 25 rad/s. Because of symmetry, only a quadrant could be considered for analysis. Initial mesh is shown in Fig. 2a and time domain is divided into five equal sub-domains and each sub-domain is divided into 10 uniform steps. The time step in each sub-domain is taken as 0.0461 for the mesh in Fig. 2a. Assuming g and gt as 1% and 1% respectively, space and time discretization errors are determined. Final refined mesh
Fig. 1. (a) Cantilever beam, (b) initial mesh and (c) final mesh. Cantilever beam subjected to suddenly applied load at quarter span.
A cantilever beam 20 m long and 2.5 m deep (Fig. 1a) subjected to suddenly applied load of 6 105 N at quarter span is considered with E ¼ 2:1 1011 N/m2 , m ¼ 0:3 and q ¼ 7850 kg/m3 . Initial mesh (Fig. 1b) is having only eight elements. Time domain is divided into five equal sub-domains and each sub-domain is divided into 10 uniform steps. The time step in each sub-domain is taken as 4:64 104 for the mesh in Fig. 1b. Assuming g and gt as 10% and 1% respectively, space and time discretization errors are determined. Refined mesh is shown in Fig. 1c and the updated time step in each subdomain is shown in Table 2. A plane stress plate 100 in: 100 in: 1 in. with a central hole of radius 5 in. is considered. The material properties are E ¼ 0:3 105 lbf/in.2 , m ¼ 0:3 and q ¼ 1:0 lbm/in.3 The plate is subjected to harmonic loading ex-
Fig. 2. (a) Initial mesh and (b) final mesh. Plane stress plate subjected to harmonic loading.
2318
A. Dutta / Computers and Structures 80 (2002) 2313–2319
is obtained after two iterations and is shown in Fig. 2b and the updated time step in each sub-domain is shown in Table 3. Table 3 Plane stress plate with a central hole subjected to harmonic loading Mesh
No. of time sub-domain, nT
g
Dt
Fig. 2a Fig. 2b
5 5
17.84 1.07
0.0461 0.0352 0.0426 0.0415 0.0425 0.0332
Table 4 Clamped plate with a suddenly applied point load at the center of the plate Mesh
No. of time sub-domain, nT
g
Dt
Fig. 3a Fig. 3b
5 5
32.34 3.11
3:726 103 3:242 103 3:221 103 3:247 103 3:199 103 3:319 103
A clamped plate 300 in: 300 in: 3 in. is loaded with a suddenly applied point load at the center of the plate. The material properties are E ¼ 3 107 lbf/in.2 , m ¼ 0:316 and q ¼ 7:324 104 lbm/in.3 A quadrant is considered because of symmetry. The initial mesh is shown in Fig. 3a and time domain is divided again into five equal sub-domains and each sub-domain is divided into ten uniform steps. The time step in each subdomain is taken as 3:726 103 for the mesh in Fig. 3a. Assuming g and gt as 3% and 1% respectively, space and time discretization errors are determined. Final refined mesh is obtained after two iterations and is shown in Fig. 3b and the updated time step in each sub-domain is shown in Table 4. 5. Conclusion The error estimation scheme for plane stress/plane strain and plate problems subjected to transient dynamic loads is computationally efficient and easily implementable. Local error measures (both in space and time) combined with adaptive meshing is very useful, converging iteratively to a solution of desired accuracy. The overall strategy is very useful in attaining high level of accuracy introducing only the required additional degrees of freedom during mesh refinement. However, the proposed time averaged error estimation will not perform well in the case of wave propagation problem as the peak stress changes its location with time as the wave propagates and hence the algorithm is particularly valid for a class of transient dynamic load, where the load location does not change. References
Fig. 3. (a) Initial mesh and (b) final mesh. Symmetric quadrant of clamped square plate with suddenly applied load at the center of the plate.
[1] Babuska I, Rheinboldt WC. Error estimators for adaptive finite element computation. SIAM J Numer Anal 1978; 15(4):736–54. [2] Babuska I, Rheinboldt WC. Adaptive approaches and reliability estimations in finite element analysis. Comput Meth Appl Mech Engng 1979;17/18:519–40. [3] Ladeveze P, Leguillon D. Error estimate procedure in the finite element method and application. SIAM J Numer Anal 1983;20(3):485–509.
A. Dutta / Computers and Structures 80 (2002) 2313–2319 [4] Ladeveze P, Pelle JP, Rougeot P. Error estimation and mesh optimization for classical finite elements. Engng Comput 1991;8:69–80. [5] Zienkiewicz OC, Zhu JZ. A simple error estimator and adaptive procedure for practical engineering analysis. Int J Numer Meth Engng 1987;24:337–57. [6] Zienkiewicz OC, Xie YM. A simple error estimator and adaptive time stepping procedure for dynamic analysis. Earthquake Engng Struct Dyn 1991;20:871–87. [7] Zeng LF, Wiberg NE, Li XD, Xie YM. A posteriori local error estimation and adaptive time-stepping for Newmark integration in dynamic analysis. Earthquake Engng Struct Dyn 1992;21:555–71. [8] Zeng LF, Wiberg NE. Spatial mesh adaptation in semidiscrete finite element analysis of linear ealstodynamic problems. Comput Mech 1992;9:315–32. [9] Wiberg NE, Li XD. A post-processed error estimate and an adaptive procedure for the semidiscrete finite element method in dynamic analysis. Int J Numer Meth Engng 1994;37:3585–603. [10] Combe JP, Ladeveze P, Pelle JP. Constitutive relation error estimation for transient finite element analysis. Comput Meth Appl Mech Engng 1999;176:165–85.
2319
[11] Schleupen A, Ramm E. Local and global error estimations in linear structural dynamics. Comput Struct 2000;76:741– 56. [12] Joo KJ, Wilson EL. An adaptive finite element technique for structural dynamic analysis. Comput Struct 1988; 30:1319–39. [13] Cook RD, Avrashi J. Error estimation and adaptive meshing for vibration problems. Comput Struct 1992; 44:619–26. [14] Dutta A, Ramakrishnan CV. Error estimation in finite element transient dynamic analysis using modal superposition method. Engng Comput 1997;14(1):135–58. [15] Ahmad S, Irons BM, Zienkiewicz OC. Analysis of thick and thin shell structures by curved finite elements. Int J Numer Meth Engng 1970;2:419–51. [16] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique. Int J Numer Meth Engng 1992; 33:1331–64. [17] Zienkiewicz OC, Boroomand B. Recovery by equilibrium in patches. Int J Numer Meth Engng 1997;40:137–64. [18] Paz M. Structural Dynamics. New Delhi: CBS Publishers and Distributors; 1987.