CHAPTER V NUMERICAL EXAMPLES
5.4 Dynamic Problems
To demonstrate the accuracy and efficiency of the proposed boundary element method, five representative examples are presented in this section. The materials, geometry, loads, and defects considered in these examples include isotropic or anisotropic, rectangular or triangular, harmonic or Heaviside-type, and hole/crack or interface.
Examples 14 and 15 consider the plate made of isotropic materials; one is steel and the other is set for comparison with [53], whose Young’s modulus E, Poisson’s ratio , and mass density
are Examples 16 considers the plate made of a T800/Epoxy unidirectional composite with the fibers oriented at 30o degree to the x1-axis, and the mechanical properties of T800/Epoxy in the material principal direction are11 22 33
where the subscript 1 denotes the fiber direction and the subscripts 2 and 3 are the directions transverse to the fiber. Example 18 is a bi-material plate consisted of two isotropic materials whose materials properties are varied based upon (5.36).
Example 14: Free vibration of an isotropic plate
Consider an isotropic plate which is fixed on one side and free on the other sides. To
vibration analysis, one is a square plate, and the other is a triangular plate. The material properties of the plate are given in (5.36). The geometry, element meshes, and boundary conditions are shown in Figure 5.25(a) and (b). Since the accuracy of boundary element analysis will depend on the element discretization, convergence study is usually presented by showing the physical response with respect to the element number. Moreover, due to the use of particular solutions in the present BEM, the effects of the number and location of the interior source points should also be reflected in the related figures. Thus, in this example the results of natural frequencies are plotted against M (element number) with several different np (number of interior source point, the same number with different location will be distinguished by a lower-case letter). The location of the interior source points is shown in Figure 5.26. To have a proper comparison, the solutions obtained from ANSYS are also shown here.
Figure 5.27(a)-(d) show that for an isotropic square plate the first four natural frequencies obtained from the present BEM (M=84, np=9) and ANSYS (M=576) are in good agreement, whose difference is within 2%. The differences of the natural frequencies obtained from np=0, np=5a, np=5b, and np=9 of BEM are within 18%. Figure 5.28(a)-(d) show the results of the first four natural frequencies for an isotropic triangular plate. The difference between the present BEM (M=67, np=16) and ANSYS (M=432) is within 4%, which is acceptable.
Through this example it is clear that the convergency of the results depend on the way we select the appropriate element discretization and interior source points. Without further detailed explanation, the results of all the following examples will then be presented using the most appropriate element discretization and interior source points.
Example 15: An isotropic plate subjected to a Heaviside-type load
A rectangular plate, with side lengths 2 m and h4 m, is made of isotropic materials whose properties are given in (5.37). The plate is fixed in normal direction but free in tangential direction on three sides and subjected to a Heaviside-type load
( ) ( ) Pa
p t H t on the fourth side (Figure 5.29), where H(t) is the step function. Because a Heaviside-type load cannot be expressed by a single harmonic function, this problem cannot be solved by (3.33) and (3.34) for the analysis of steady-state forced vibration.
Instead, due to the inclusion of the complete spectrum of frequencies in a Heaviside-type load, it is a standard loading case that tests the performance of the numerical methods for transient analysis. As stated in Section 3.2.5, the numerical evaluation of the ordinary differential equations formed by BEM can be performed by Houbolt’s algorithm or modal superposition method. Therefore, the numerical results shown in this example include those of BEM by these two different techniques. Also, for the purpose of comparison, the results calculated by the dual reciprocity BEM using new Fourier radial basis functions [53]
are also presented.
For the mesh of the plate boundary with 12 elements, the selection of 3 interior points for the particular solution, the use of the time step t 0.001 sec. and the initial conditions u(0)0 and (0)u 0 for all nodes, the numerical results of vertical displacement at point A and normal traction at point B are shown in Figure 5.30. From this Figure we see that the results begin to diverge after few time steps by using Houbolt’s algorithm. Except the unstable initial region (0<t<0.01) of Houbolt’s algorithm, all the numerical values are in good agreement.
Example 16: An anisotropic elastic plate with an elliptical hole subjected to a harmonic load
Consider a square plate with an elliptical hole subjected to a harmonic load ( ) cos100 GPa
p t t on the upper edge of the plate. The plate is made of a unidirectional composite whose mechanical properties are given in (5.38). The geometry, element meshes, interior points for particular solutions, boundary and loading conditions are shown in Figure 5.31. Table 5.3 shows the variation of the first four natural frequencies with respect to the ratio of b/a. Figure 5.32 shows the mode shapes corresponding to the first four natural frequencies for the problem with circular hole, i.e., b/a=1. These results show that the natural frequencies calculated by the present BEM and ANSYS well agree with each other. The main difference is that only 48 elements are needed for the present BEM to attain an accurate solution, whereas fine meshes are needed for ANSYS. This improvement is due to the fact that no meshes are required along the boundary of the elliptical hole for the present BEM.
Figure 5.33 and Figure 5.34 show the time variation of the vertical displacement at point C and the stresses at point D, calculated by the present BEM using modal superposition method and ANSYS, with initial conditions (0)u 0 and (0)u 0 for all nodes and a/b=1 for circular holes. Since the applied load is a harmonic load, the steady-state solution (3.34) without considering the effects of initial condition is also shown in Figure 5.33 and Figure 5.34. In real application, due to the existence of damping the complementary solutions associated with natural frequencies will die out with time, and the solutions obtained from the transient analysis should approach to the steady-state solution as time goes on. Although the transient responses shown in Figure 5.33 and Figure 5.34 do not approach to the steady-state solutions, they are reasonable since the damping effect is not considered in the present study. Both of the transient and steady-state solutions show the periodicity of the applied load ( )p t cos100t, i.e., T 2 /100
0.0628sec.Whereas the local period T 9.96 10 3sec. is observed for the transient responses, whose associated frequency 2 / T 630.7 is close to the second lowest natural frequency 2 630.9 (see Table 5.3). The mode shape corresponding to 2 shown in Figure 5.32 reveals that this is a mode dominated by the load in x2-direction, which agrees with the load applied in this example.
Theoretically, the stress 11 at point D should be zero due to the requirement of the traction-free condition around the hole. However, only the steady-state solution obtained by the present BEM (see Figure 5.34 (a)) meets this requirement. The transient solution of
11 of the present BEM is very close to zero (see Figure 5.34 (a)), whereas both of the transient solutions calculated by ANSYS (see Figure 5.34 (b)) are not so close to zero as the present BEM. Nevertheless, they are all relatively less than 22. In other words, they are all well approximated to zero, but the approximation of the present BEM is better than ANSYS.
Example 17: Dynamic analysis of a plate with a crack
To show the influence of crack length (a), or crack orientation (), or crack location () on the natural frequencies of a plate, a square plate made of isotropic material, whose properties are defined as Eqn.5.36, contains a crack with different parameters such as a, and is considered here (Figure 5.35). By varying the crack length, or crack orientation, or crack location, the numerical results of natural frequencies of the first four modes are shown in Table 5.4-Table 5.6. Their associated vibration mode shapes are shown in Figure 5.36 for the case of 2a/L=1/3, L=6 m, =0o, =3 m. For the purpose of comparison, some results are also calculated by ANSYS. The differences between the present BEM and
ANSYS are all within 2%, which is acceptable. The effects of crack length, orientation and position can then be studied from the last two columns of Table 5.4-Table 5.6, which show the ratio of reduced frequencies due to the appearance of cracks. With the aids of the associated mode shapes shown in Figure 5.36, reasonable explanation of these values can be made for the relatively larger or smaller numbers shown in these Tables. For example:
(1) Table 5.4 shows the longer the crack length is, the smaller the natural frequency is, which is reasonable since the crack could soften the plate. (2) Table 5.4 shows that the length effect is much more apparent in the second mode. This is also reasonable since Figure 5.36 shows that the second mode is a tensile mode, which is more sensitive than the other modes if the crack is lying horizontally. (3) Table 5.5 shows the opposite trend on the ratios of reduced frequency for mode 2 and mode 4. The former decreases with respect to the angle of orientation, and the latter behaviors on the contrary. This is conceivable by knowing that the most precarious condition is the tensile mode. The tensile mode of horizontal crack is the mode 2, whereas that of the upright crack is the mode 4. (4) Table 5.6 shows that the closer between the crack and the fixed end, the lower the first three natural frequencies. Opposite trend occurs for mode 4. It seems odd for this opposite trend.
Whereas when we check their associated mode shapes we found that the crack in different location may result in different mode shape for higher vibration modes, and hence the comparison made for higher modes should be carefully.
Besides the parameter discussed above, an anisotropic elastic plate with a crack (a2.4 mm, w20 mm, =40 mm ) subjected to a Heaviside-type load p t( )H t( ) GPa in the plane stress condition is performed (Figure 5.37). This example had already been treated by Zhu et al. [56] who used a finite difference formulation and by Albuquerquea et al. [57] who used BEM with subregion technique. Due to the symmetry of this problem,
only one-quarter of the plate is analyzed by Zhu and Albuquerquea where 17 boundary elements and 12 internal points are used for the latter one, whereas a full model is constructed for present BEM where 24 boundary elements and 32 internal points are used.
The material properties are
4 4
The dynamic stress intensity factors are calculated from the near-tip stress formula shown in Eqn. (5.4) and the normalized dynamic stress intensity factor K K , where I / 0
0 0
K a, are plotted in Figure 5.38. The results show that the one of present BEM agrees with them obtained from the references. In this example, the near-tip stresses are accurate and the stress intensity factor KI from Eqn. (5.4) with different values r=0.01 mm and r=0.001 mm are convergent enough to this problem because of the special fundamental solutions used and there is no element around the crack surfaces.
Example 18: A bi-material plate subjected to a Heaviside-type load
Consider a bi-material plate made by two dissimilar isotropic materials subjected to a Heaviside-type load ( )p t H t( ) GPa at the upper ends (Figure 5.39). The elastic properties of the material above the interface are given in (5.36), whereas those of the material below the interface are varied by a ratio k for Young’s modulus, i.e., E2 kE1 and all the other properties remain the same. The element meshes, interior points for particular solutions, and boundary conditions are shown in Figure 5.39. Figure 5.40 shows the variation of the first four natural frequencies with respect to the ratio k. Figure 5.41 shows the mode shapes corresponding to the first four natural frequencies for the plate with k=5. The response of the vertical displacement at point F is then shown in Figure 5.42. All
these results show that the vibration analysis of a bi-material plate can be performed successfully by the present BEM with no meshes along the interface.
CHAPTER VI CONCLUSIONS
In this dissertation, the elastostatic fundamental solutions for hole, crack, inclusion and interface problems based on Stroh formulism are applied successfully in the fields of piezoelectric, viscoelastic materials and dynamic BEM. By comparing the mathematical forms of the basic equations in elasticity, piezoelectricity and viscoelasticity in Laplace domain, it can be found that there are no differences between the mathematical forms of those. The only different is the contents and dimension of the matrix and vectors. Hence, by the basic knowledge discussed in Chapter II and III, the fundamental solutions and BIEs of piezoelectric or viscoelastic problems can be performed analogous to elastostatic one by carrying out some transformations. Moreover, by the aid of dual reciprocity, the static fundamental solutions can be used again for dynamic problems. Because the fundamental solutions used satisfy the boundary conditions exactly, no elements are needed along the specific boundary, which is shown in the examples of Chapter V. From a series of numerical examples some conclusions are given in the follows.
For the static analysis in hole, crack and inclusion problems, only 8 elements were used in present BEM to obtain accurate results, but in the interface problems (or called interfacial corner problems) more elements are needed (compared with hole, crack and inclusion problems) to grasp the behaviors near the interface. However, comparing with those in the traditional BEM or FEM, present BEM is more efficient and accurate especially for the case of infinite domain containing a crack. For viscoelastic problems, although the numerical inversion may raise the concerns of inaccuracy and
merits which attribute to the exact satisfaction of boundary conditions overtake the possible drawbacks induced by the numerical inversion. For the dynamic analysis of the free vibration cases, fine mesh is needed to perform with the accurate results of higher mode. The reason is that the mode shapes of higher mode are usually complex, hence the higher degrees of freedom is necessary.
Through the all numerical examples in Chapter V, the supremacy in accuracy and efficiency of present BEM is proved by comparing with FEM or traditional BEM. As previously concluded, in the cases of static hole/crack/ inclusion problems only 8 boundary elements are needed to attain precise results. Maybe the special boundary elements, which adopt the specific fundamental solutions satisfying the boundary conditions, can be the fundamental elements of FEM to improve the accuracy of the analysis.
REFERENCES
[1] Crouch, S. L., 1976, “Solution of plane elasticity problems by the displacement discontinuity method. I. Infinite body solution,” International Journal for Numerical Methods in Engineering, 10, pp. 301-343.
[2] Blandford,, G. E.,†, Ingraffea, A. R., Liggett, J. A., 1981, “Two-dimensional stress intensity factor computations using the boundary element method,” International Journal for Numerical Methods in Engineering, 17, pp. 387-404.
[3] Hong, K. H., Chen, J. T., 1988, “Derivations of Integral Equations of Elasticity,”
Journal of Engineering Mechanics, 114, pp. 1028-1044.
[4] Cruse, T. A., Snyder, M. D., 1975, “Boundary-integral equation analysis of cracked anisotropic plates,” International Journal of Fracture, 11, pp. 315-328.
[5] Hwu, C. and Yen, W. J., 1991, “Green’s functions of two-dimensional anisotropic plates containing elliptic hole,” Int. J. Solids and Structures, 27, pp. 1705-1719.
[6] Hwu, C. and Yen, W. J., 1992, “Plane problems for anisotropic bodies with an elliptic hole subjected to arbitrary loadings,” The Chinese Journal of Mechanics, 8, pp.
123-129.
[7] Clements, D. L., Haselgrove, M. A., 1983, “A boundary integral equation for a class of crack problems in anisotropic elasticity,” Int. J. Comp. Math., 12, pp. 267-278.
[8] Snyder, M. D. and Couse, T. A., 1975, “Boundary-integral equation analysis of cracked anisotropic plates,” Int. J. Fracture, 11, pp. 315-328.
[9] Hwu, C. and Yen, W.J., 1993, “On the anisotropic elastic inclusions in plane elastostatics,” Journal of Applied Mechanics, 60, pp. 626-632.
[10] Hwu, C. and Liao, C.Y., 1994, “A Special Boundary Element for the Problems of Multi-Holes, Cracks and Inclusions,” Computers and Structures, 51, pp.23-31.
[11] Ting, T.C.T., 1996, Anisotropic Elasticity, Oxford University Press, Oxford.
[12] Hwu, C., 2010, Anisotropic Elastic Plates, Springer, New York.
[13] Kuo, C. M. and Barnett, D.M., 1991, “Stress singularities of interface cracks in bonded piezoelectric half-spaces,” Modern Theory of Anisotropic Elasticity and Applications, pp.33-50.
[14] Suo, Z., Kuo, C. M., Barnett, D. M. and Willis, J. R., 1992, “Fracture Mechanics for Piezoelectric Ceramics,” Journal of the Mechanics and Physics of Solids, 40, pp.
739-765.
[15] Liang, Y. C. and Hwu, C., 1996, “Electromechanical analysis of defects in piezoelectric materials,” Smart Materials & Structures, 5(3), pp. 314-320.
[16] Hwu, C., 2008, “Some explicit expressions of extended Stroh formalism for two-dimensional piezoelectric anisotropic elasticity,” International Journal of Solids and Structures, 45(16), pp. 4460-4473.
[17] Chen, Y.C. and Hwu, C., 2010, “Green’s Functions for Anisotropic/Piezoelectric Bimaterials and Their Applications to Boundary Element Analysis,” Computer Modeling in Engineering & Sciences, 57, 1, pp. 31-50.
[18] Carini,A. and Gioda,G., 1986, “A boundary integral equation technique for visco-elastic stress analysis,” International Journal for Numerical and Analytical Methods in Geomechanics,10, pp. 585-608.
[19] Sensale,B. and Greus,G.J., 1993, “Boundary elements analysis of viscoelastic fracture,” Boundary Element XV: Stress Analysis, pp. 291-301.
[20] Sim, W.J. and Kwak, B.M., 1988, “Linear viscoelastic analysis in time domain by boundary element method,” Computers and Structures, 29, pp. 531-539.
[21] Mesquita, A.D., Coda, H.B. and Venturini, W.S., 2001, “Alternative time marching process for BEM and FEM viscoelastic analysis,” International Journal for Numerical Methods in Engineering, 51, pp. 1157-1173.
[22] Mesquita, A.D. and Coda, H.B., 2001, “An alternative time integration procedure for Boltzmann viscoelasticity: A BEM approach,” Computers and Structures, 79, pp.
1487-1496.
[23] Schanz, M., 1999, “A boundary element formulation in time domain for viscoelastic solids,” Communications in Numerical Methods in Engineering, 15, pp. 799-809.
[24] Manolis, G.D. and Beskos, D.E., 1988, Boundary element methods in elastodynamics, Unwin-Hyman, London.
[25] Nardini D. and Brebbia C.A., 1982, “A new approach to free vibration analysis using boundary elements,” Proceedings 4th International Conference on Boundary Element Methods.
[26] Partridge, P.W., Brebbia, C.A. and Weobel, L.C., 1992, The dual reciprocity boundary element method, Computational Mechanics Publications.
[27] Partridge, P.W., 2000, “Towards criteria for selecting approximation functions in the dual reciprocity method,” Engineering Analysis with Boundary Elements, 24, pp.
519–529.
[28] Kögl, M. and Gaul, L., 2000, “A 3-D boundary element method for dynamic analysis of anisotropic elastic solids,” Computer Modeling in Engineering & Sciences, 1, pp.
27-43.
[29] Kögl, M. and Gaul, L., 2003, “Free vibration analysis of anisotropic solids with the boundary element method,” Engineering Analysis with Boundary Elements, 27, pp.
107–114.
[30] Bathe, K.J. and Wilson, E.L., 1973, “Stabilityand accuracy analysis of direct integration methods,” Earthq. Engng Struct. Dyn., 1, pp. 283-291.
[31] Hilber, H.M., Hughes, T.J.R. and Taylor, R.L., 1977, “Improved numerical
Engineering & Structural Dynamics, 5, pp. 283-292.
[32] Houbolt, J.C., 1950, “A recurrence matrix solution for the dynamic response of elastic aircraft,” J. Aeronaut. Sci., 17, pp. 540-550.
[33] Brebbia, C.A., Telles, J.C.F. and Wrobel, L.C., 1984, Boundary element techniques : theory and applications in engineering, Springer-Verlag, Berlin , New York.
[34] Agnantiaris, J.P. and Polyzos, D. and Beskos, D.E., 1998, “Three-dimensional structural vibration analysis by the Dual Reciprocity BEM,” Computational Mechanics, 21, pp. 372-381.
[35] Haddad, T.M., 1995, Viscoelasticity of Engineering Materials, Chapman & Hall.
[36] Hoop, A.T. de, 1966, “An elastodynamic reciprocity theorem for linear, viscoelastic media,” Appl. Sci. Res., 16, pp. 39-45.
[37] Rogacheva, N.N., 1994, The Theory of Piezoelectric Shells and Plates, CRC Press, London.
[38] Lekhnitskii, S.G., 1968, Anisotropic Plates, Gordon and Breach Science Publishers, New York.
[39] Reismann, H. and Pawlik, P.S., 1980, Elasticity: Theory and Applications, John Wiley & Sons, New York.
[40] Schapery, R.A., 1962, “Approximate methods of transform inversion for viscoelastic stress analysis.” Proc.-4th U.S. National Congress on Applied Mechanics, pp.
1075-1085.
[41] Gaul L., Kögl M. and Wagner M., 2003, Boundary element methods for engineers and scientists, Springer, New York.
[42] Wylie, C.R., 1975, Advanced Engineering Mathematics, McGraw-Hill, New York.
[43] Blandford, G. E., Ingraffea, A. R. and Liggett, J. A., 1981, “Two dimensional stress intensity factor computations using the boundary element method,” Int. J. Numer.
Meth. Engrg., 17, pp. 387-404.
[44] Pan, E. and Amadei, B., 1999, “Boundary element analysis of fracture mechanicsin anisotropic bimaterials,” Engineering Analysis with Boundary Elements, 23, pp.
683–691.
[45] Pan, E., 1999, “A BEM analysis of fracture mechanics in 2D anisotropic piezoelectric solids,” Engineering Analysis with Boundary Elements, 23, pp.
67–76.
[46] Narayanan, G.V. and Beskos, E., 1982, “Numerical Operation Methods for Time-Dependent Linear Problem,” International Journal for Numerical Methods in Engineering, 18, pp. 829-1854.
[47] Carini, A. and Gioda, G., 1986, “A boundary integral equation technique for
[47] Carini, A. and Gioda, G., 1986, “A boundary integral equation technique for