A Six-node Triangle Finite Volume Method for Solids with a 一六节点三角形有限体积法的固体和

A Six-node Triangle Finite Volume Method for Solids with a 一六节点三角形有限体积法的固体和 A Six-node Triangle Finite Volume Method for Solids with a Rotational Degree of Freedom * , Marcus WheelWenke Pan Department of Mechanical Engineering, University of Strathclyde, Glasgow, UK. Abstract: A novel vertex centred finite-volume (FV) method is presented for determining the displacement field in loaded two-dimensional structures. The solid is discretized into six-node finite elements (FE) and the standard six-node FE shape function is employed to describe the displacement variation in each element. However, the displacements of the midside nodes are replaced with the displacements and rotations of the two neighbouring edge nodes. Thus, only the control volumes centred around each vertex node of the triangular element are required. These control volumes are constructed by joining lines between the centre of the element and the midside nodes. The force equilibrium in two perpendicular directions and one in-plane moment equilibrium equation are derived for each control volume. By solving the system equations, the displacements and rotations of the vertex nodes can be obtained and then the displacements of the midside nodes can be calculated. Following this, the strain and stress can be determined through the relationship between strain and displacement, and stress and strain, respectively. The results of simulation show that this novel FV method is not only better than the previously derived three-node triangle FV method with rotational degrees of freedom, but also better than the three-node FE method with rotations, and can also compete with the six-node triangle FE method. In addition, no locking problem was encountered as the value of Poisson’s ratio of the elastic material approached 0.5. 1 Keywords: Vertex centred finite volume method; Rotational degree of freedom; Six- node triangle element; Moment equilibrium equations. 1. Introduction The finite element method has dominated the simulation of solids and structures for many years, while the FV method has been used successfully for fluid dynamic analysis. Over the last fifteen years, there has been some progress is using the FV method for solids and structures. Initially, the FV method was used for small strain elastic compressible and incompressible material for 2D plane strain, plane stress, axial-symmetry and 3D analysis [1-6], and also dynamic analysis [7]. Following this, it was used for plate analysis [8-9] and then for large deformation hyperelastic [10-11] and elastic-plastic materials [12]. The FV method has also been used for coupled analysis of interacting fluids and structures in which the same discretization method is used for both media [13-14]. This unified method may show some potential advantage for this type of analysis. Two common methods to construct the FV control volume are the cell vertex method and the cell centre method. For the cell-centre method, the displacement field of the shape function cannot be higher than second order, whilst for the cell-vertex FV method, the standard FE shape function can be used directly. Some special FE methods incorporating a rotational degree of freedom(DOF) for 2D analysis have been presented and their shape functions can also be used for FV analysis. 2 Allman [15-16] presented a compatible triangular finite-element including vertex rotations for plane elasticity analysis. When replacing the middle-node displacement with the displacements and rotations of the neighbouring vertex node, Cook et al [17] compared the results between 3-node, 6-node triangular and 4-node and 8-node quadrilateral finite elements and pointed out that there are some advantages for the inclusion of rotations in some cases. Onate et al [18] presented an FV method for beam analysis which includes a rotation, the discritization procedure involving the four types of equations of bending moment, shear force, rotation and deflection. Based on [15], Wenke et al. [19] included the Allman rotation for 3-node FV analysis for 2D elasticity problems, the results showing that there are some advantages over the 3-node FE method with or without rotations and the FV method without rotations for some cases. In this paper, a 6-node triangular finite-element shape function was used for 2D elastic analysis. To incorporate rotations in the present FV analysis, control volumes were constructed around each vertex node of the six-node triangle only; there are no volumes around the midside node unlike other reported methods [20]. Two force and a moment equilibrium equations were written for each control volume. By using the assumption in [17], the total system of equilibrium equations can be transformed to including only the vertex node displacements and rotations and thus the DOF of each triangular element is then reduced from 12(6*2, 6 nodes with displacement components u and v) to 9 (3*3, 3 vertex nodes with u, v and rotation ,). Because of the special structure of the CV, the total number of unknown variables is equal to the total number of equilibrium equations. By solving these equations, all unknown edge displacements and rotations can be obtained, and then the midside node displacement can be determined. The second 3 section of this paper gives an overview of the triangular element and the third section discusses the construction of the CV, gives an overview of the formulation and details the displacement and stress boundary conditions. The fourth section provides some examples for testing the method and also a comparison between the different methods. Finally, conclusions are given. 2. Displacement, strain and stress fields for the 6 node triangular finite element 2.1 Displacement fields Fig. 1 shows a six-node 2-D triangular element. A conventional six-node triangular polynomial shape function is used and the displacement components, u and v, within the triangles are assumed to be 22,,，，，，，uaaxayaxaxyay123456 (1) ,22v,b，bx，by，bx，bxy，by,123456 where x and y are the Cartesian coordinates of any point in the interior of the element. Applying equations (1) at every node of the triangle, gives: ,1,a,Au (2) ,,1b,Av, TTT,,,,,,a,aa...ab,bb...bu,uu...uwhere , and 126126126 T,,v,vv...v, while matrix A is given in the Appendix. 126 2.2 Strain fields The expressions of the strain within the triangle can be deduced from the strain and displacement relationship to provide the following: 4 ,,u,,,Ba,xxx,x,,v,, (3) ,,Bb,yyy,y,,v,u,,,，,Bb，Baxyxy,,x,y, where , and are the usual direct and shear strain components and matrices ,,,yyxyxx and are given in the Appendix. BByx 2.3 Stress fields Elastic isotropic or anisotropic material can be considered, using Hook’s law, for the stress and strain relationship of the material: (4) σ,Dε TT,,σ,,,,,,ε,,,,where the stress and strain vectors are and xxyyxyxxyyxyand T，，DDDD1,,，，111213,,T,, (5) DDDD2,,,,D212223,,,,T,,,,,,DDDD3313233，，，， The vector , and are the first, second and third row of the matrix D ,,,,,,D1D2D3 respectively. Thus: T,,,,,D1ε,,,,xx,,,,T, (6) D2σ,,Dε,ε,,,,,,,,,,,,,,yyT,,,,,,,,,D3εxy,,,, 3. Control volume and the equilibrium equations 3.1 Control volume 5 The conventional cell vertex based FV method is incompatible with the six-node triangle element unless the construction of the CV’s is modified. An obvious modification is shown in Fig. 2(a). Another approach is described in [20]. There is a control volume surrounding every node including the midside ones. The control volumes are multi-faced polygons and are constructed by connecting the centre of the triangles to points midway between the midside and vertex nodes. For each vertex node, there are two CV edges located within each of the elements connected to it. The multi-faced polygons are the control volumes associated with a vertex nodes (edge node A and midside node B) located inside the solid and a vertex node (Edge node C) on the boundary. For a 2D problem, each control volume has two force equilibrium equations along perpendicular directions and each node has two unknown displacement variables, so that the total number of equations is equal to the total number of unknown variables. Solving these equations, the unknown nodal displacements can be obtained and the strains and stresses determined subsequently. In this paper, rotational degrees of freedom are incorporated in the analysis (see Fig.3), by replacing the midside node displacements with the displacements and rotations of the vertex nodes. In so doing, the final unknown variables will consist of only the displacements and rotations of the vertex nodes. The modified method for the construction of the control volumes in a mesh of 6 node triangles cannot be used here, however conventional vertex based CVs can be constructed. There are several ways to build these. The easiest way to construct them for this purpose is shown in Fig. 2(b). The control volume is constructed by the connection of the centre point of the neighbouring triangular element to the neighbouring midside node. The control volume for an internal node, A, and a boundary edge node, C, are shown in Fig. 2(b). Fig. 4 (a) 6 details the control volume around a vertex node. Each control volume has three equations that include two forces and a moment equilibrium equation. Thus, the total number of discrete equilibrium equations will be the same as the total number of unknown variables that comprise two displacements and a rotation for each node. 3.2 Equilibrium equations and stiffness matrix From Fig. 4(b), for any point P(x,y) along the line AB: ,,,xxSsin,A (7) ,y,y，Scos,A, where is the angle between the normal of line AB and the x-axis and S is the distance , from point A to P. The resultant force components and acting in x and y directions are: FFyx SAB,,,,,,，F(cossin)ds,,Xxxxy0 (8) ,SABF,(,cos,，,sin,)ds,,Yxyyy0, Also, the resultant moment about the origin of the coordinate system is: SSABAB (9) M,,,(cos,，,sin,)yds，(,cos,，,sin,)xdsoxxxyxyyy,,00 Also, the resultant forces and moments can be integrated analytically, to provide the following: ，，B0ds,x,1uA0，，,,,,TT,，，,,,,，0BFD1cos,D3sin,ds (10.1) ,,,Xy,,,1,,v0A,,，，,,BBdsds,,yx，， ，，B0ds,x,1uA0，，,,,,TT,，，,,,,，0BFD2sin,D3cos,ds (10.2) ,,,Yy,,,1,,v0A,,，，,,BBdsds,,yx，， 7 and: M,o ,,，，Byds0,x,,,,TT,,，，,,,,,,,D1cos，D3sin0Byds，,y,,,,,1,,BydsByds (10.3) u,,A0，，yx,,，，,,,,,,,1,,v0A，，B0xds,,，，,x,,,,TT,,，，,,,,D2sin,，D3cos,0Bxds,y,,,,,,,,BBxdsxds,,yx，，,, where the derivation of the matrices, , and are given BxdsBydsBxdsByds,,,,yyxx in the Appendix. For each control volume, the forces acting in the x and y directions and the moment about the origin of the coordinate system must be equivalent to those acting externally, where these can be expressed as follows: j,n,jEF，F,0,XX,j1,j,n,,jE (11) F，F,0,,YYj1,,j,njE,MM0，,,OO,j1,, EEEwhere , and M are the sum of any external forces and moments applied to the FFOXY control volume, and n is the total number of CV internal edges around each boundary edge node. The moment equilibrium equation for each volume is necessary to enable the rotation to be introduced and this is different from the standard finite-volume method for solids. For the conventional 2D FV method for solids, there is no need to consider the moment equilibrium equation, but for the FV method presented in this paper, the impending introduction of the extra unknown variable (the rotation), requires that the moment equilibrium equation must be considered. 8 By assembling all equations of every control volume the total equilibrium equations of the problem considered can be written as follows: uvuvEKKF,，，,,XuXvXu,,,,,,uvuvE (12) KKF,,,,,,YuYvY,,v,,uvuvE,,,,MMM,mumvO，，,, The unknown variables in the above equations are the displacements of all element nodes and the stiffness matrix is therefore not square. The displacements of the midside nodes of each element can be replaced by the displacements and rotations of the two neighbouring vertex nodes. Assuming that the transformation matrix is T [17], the following equation shows the relationship between the unknown variables of the midside nodes and that of the vertex nodes: TT,,,,u...uv...v,T*uuuvvv,,, (13) 1616123123123 For each finite element the linear transformation matrix T is a 12*9 matrix and can be found in [17]. Substituting (13) into (12) finally gives the system equations as follows: EKKKuF,,,，，,,,XuXvXX,,,,,,EKKKvF,, (14) ,,,,,YuYvYY,,E,,,,,,KKKωM,,MuMvMo，，,,,, where: TT,,,,uvω,u...uv...v,..., 111npnpnp TTEEEEEEEEE,,,,FFM,F...FF...FM...M 111XYOxnpynpnp and ‘np’ is the total number of vertex nodes. By solving Eqn. (14), the vertex displacements and rotations can be obtained, after which, using Eqn. (13), the midside node displacements can be calculated. From Eqns. (3) and (6), the strain and stress field can be determined subsequently. 9 4. Boundary conditions Displacement boundary conditions can be one or all of the three types: 000u,uv,v,,w, and pppppp They can be applied directly to the system equations (14) as in the conventional FE method. The force boundary can be a concentrated point force or a distributed stress. For a 000concentrated force at the edge node, , and . They can be F,FF,FM,MXpXpYpYppp put into Eqn. (14) directly in the usual way. For a distributed stress, for example the uniformly distributed stress shown in Fig. 5: , and F,,Scos,F,,Ssin,XA0ACYA0AC M,,,Scos,*(x，x)/2，,Ssin,*(y，y)/2OA0ACCA0ACCA which in effect enable the distributed loads to be replaced by equivalent nodal forces and moments. 5. Numerical results A number of benchmark problems, including simple tension and simple shear are considered. Then convergence and performance near incompressibility are tested for a square plate under a pure bending loading. Then the case of a rectangular plate with different ratios of length-to-width subjected to pure bending is investigated. Finally, a comparison of the present FV method and previous FV methods and FE methods for a cantilever beam subjected to varied load is performed. The values of Young’s modulus, E, and Poisson’s ratio, ν, were chosen as 1.0GPa and 0.3 respectively for test one to test four (but not for the near incompressibility tests). For the last test, the material properties are the same as that of the reference [15] for simplicity of the comparison. 10 5.1 Test 1: simple tension of a square plate The first test problem is the simple tension of a square thin plate with the length and width both being 2.0m. A quarter model is simulated, shown in Fig. 6, for which the =0, V=0, U=0,V=0 and also ,=,=,=0. A constant boundary conditions are UAABCABC horizontal tensile stress of 1MPa is applied on the right side of the plate by imposing equivalent forces to the right edge nodes and moments to the corner nodes as required. The whole model is discritized as four 6-node elements and hence there are five corresponding control volumes. Table 1 shows the comparison of the calculated and theoretical results of displacement, strain and stress at point D(0.75m, 0.2m). As expected the results are exact. 5.2 Test 2: simple shearing of a square plate The simple shear model is shown in Fig. 7. Constant shear stresses are applied along all the plate(size 1m*1m) sides, the value of the shear stress being 1.0MPa. The boundary conditions are as follows: U=0, V=0, U=0, ,=0. The same discritization is used as AABA for Test 1, and the investigation point is point D(0.75m, 0.2m). Table 2 shows the comparison of the simulated results and the exact results. From Table 1 and Table 2, the simulated results are noted to be the same as the exact results within the limits of computational accuracy. Actually, by investigating the entire set of results for both models, a linear displacement field and a constant strain and stress field are observed. Thus, both of the patch tests have been passed. Although the midside node displacements are replaced by the neighbouring edge node displacements and rotations, due to the assumption of quadratic displacement field (see equation (1)), the displacement field within the 6-node triangle element is still in quadratic form, the simulation results can predict the exact results for constant stress cases. 11 5.3 Test 3: pure bending of a square plate 5.3.1 Convergence The material and geometry of the third benchmark problem are the same as those of Test 2, but this test is not that of a constant-stress state. For the convergence test, a series of mesh refinements were performed. The whole model is discritized into 6-node triangular elements with different mesh sizes. A typical mesh is shown in Fig 8. The meshes were all constrained by suppressing the edge normal displacements and rotations of all of the left edge nodes as well as the tangential displacement of the node midway along the edge. A linear normal stress distribution, varying from 1.0MPa at the top corner to –1.0MPa at the bottom corner, was then applied to the right edge of the plate by applying suitable edge normal forces and moments to the right edge nodes. Table 3. shows the displacement component values predicted at the lower right corner of the plate with mesh refinement by the present method and that with previous two-field and three-field FV formulations. The results of U and V displacements with the present method and also the exact results are shown in the fourth and the seventh columns respectively. From these results, it can be seen that the simulated values converge to the exact values. For the coarsest mesh 4*4, by considering the U displacement, the relative error of the current method is less than 0.7%, while for the two-field formulation, it is 11.5% and for the three-field formulation, it is 1.4%. For the V displacement, the relative error of the present method and two-field and three-field are 0.6%, 10.9% and 1.3% respectively. Comparing all the cases in this table, the results show that at any given level of mesh refinement, the present method is more accurate than its three-field counterpart when 12 the number of degrees of freedom used in each method is the same. The number of degrees of freedom of the present method is about 50% greater than that of the two-field formulation with same mesh, but the results for the U and V displacements of the present method with the 8*8 mesh (135 degrees of freedom) are better than that of two-field 16*16 mesh(degrees of freedom 578). From Table.3, it can be seen that present method is better than both the two-field and previous three-field finite-volume methods. 5.3.2 Poisson’s ratio The same model as 5.3.1 is used for this test except Poisson’s ratio is now varied. When the value of Poisson’s ratio for the material approaches 0.5, the 3-node triangular finite element will exhibit locking. To investigate if this behaviour occurs for the present method, the value of Poisson’s ratio varying from 0.49 to 0.49….9 (15 digits of ‘9’) is performed. Fig. 9 shows the top right corner node displacements versus the logarithm function of Poisson’s ratio. The theoretical solution for completely incompressible material (Poisson’s ratio equals 0.5) is also shown in this figure. From both Fig. 9 (a) and (b), it can be seen that as the mesh is refined the computed results tend towards the exact results at all values of ν investigated. Clearly, the results do not appear to diverge from the exact solution as Poisson’s ratio approached 0.5. 5.4 Test 4: Pure bending of a slim cantilever with a high aspect ratio The fourth verification test involved predicting the deformation of a rectangular plate with an aspect ratio of 8 (L/H=8.0) under a state of pure bending. Two meshes with nx*ny of 8*2 and 16*4 were used for the comparison of the three FV methods. In this case the number of elements in the horizontal direction was four times the number in the 13 vertical direction. Thus, the aspect ratio or length of the horizontally aligned side relative to the vertically aligned side of each element was 2.0. The same boundary conditions as those of the previous case were used. The right edge had suitable forces and moments applied so as to impose a linear normal stress distribution, varying from 1.0MPa at the top right corner to –1.0MPa at the bottom right corner. Table 4 shows the comparison of the results from the three FV methods. From this table, for 81 degrees of freedom, the U and V displacements of the present method are much better than that of two-field FV method with 170 degrees of freedom. For the 8*2 mesh, the values for the present FV method are the same as those for three-field FV method, for only one layer of element above and under the neutral axis of the plate, while for the 16*4 mesh, the values for the present method are worse than those for the previous three-field method. This may be caused by the replacement of the middle node displacements by the displacements and rotations of the two neighbouring edge nodes and the greater aspect ratio of the horizontally aligned side being twice that of the vertically aligned side of the element. Although the quadratic displacement fields have been used for this study, the linear stress fields should be predicted exactly for pure bending cases, the simulation results only converge with mesh refinement. The reason is that the rotation boundary conditions used are not exactly the same as the true ones. 5.5 Test 5: bending of a cantilever subjected to an end load The last test problem is a cantilevered rectangular plate subjected to a transverse end load. This problem was previously used to assess the performance of FE procedures incorporating rotations [15, 17]. The geometry of the plate and the associated material 14 properties are shown in Fig. 10. For direct comparison with previous FE and FV results, the imperial units were not converted into metric units. The two meshes employed previously [15, 19] were again used to represent the geometry of the plate. The 3 lbs, was administered by applying a parabolic shear transverse end load, W = 40×10 Wy32stress distribution () across the right end of the cantilever. ,,[1,4()]xyx,LHtH2 The equivalent forces and moments were applied to the nodes on the right edge of the plate. The left edge of the plate was constrained in the manner employed in the preceding tests. Table 5 shows a comparison the results of the present FV method and those of the other FE and FV methods, which include the constant strain triangle (CST), the Allman triangle, the linear strain triangle(LST), a two-field FV and the three-field FV [15, 19]. For the present method with 81 degrees of freedom, the prediction of the vertical deflection of the lower right corner of the plate is 0.3483mm, which is better than that of the CST with 170 DOF and of the Allman Triangle and the two-field and three-field FV with the same degrees of freedom. With 255 DOF, the present method is also better than the aforementioned methods. For the comparison between the LST and the present method, it is seen that they are both converging quickly to exact results, with the LST performing slightly better than the present method. 6. Conclusions A new FV procedure for predicting deformations in two-dimensional plane problems has been presented. Control volumes centred on the vertices of a triangular FE mesh have been used for the construction of the force and moment equilibrium equations. In contrast to previous work the six-node standard polynomial shape function is used for 15 the displacement field. The displacements of the midside nodes are replaced by the displacements and rotations of the neighbouring vertex nodes. By using the side CVs and the replacement of the unknown displacements of the midside nodes, the total number of system unknown variables is the same as the total number of equilibrium equations. The simulation results show that the present method passes the patch tests with constant strain states and converges to the exact results for linear strain states, under mesh refinement. No locking phenomenon is observed when the value of Poisson’s ratio approaches 0.5. The method is better than two-field finite volume methods and also the constant strain triangle and the Allman triangle based finite element methods. It is also better than three-field FV method for regular mesh and can also compete with the linear strain triangle finite element method. Appendix The formulation of matrix A in Eqn. (2) is as follows: 221xyxxyy，，111111,,221xyxxyy222222,,A, (A-1) ,,..................,,221xyxxyy666666，， Band in Eqn. (2) are given as follows: Byx T,,B0102xy0,,,x (A-2) ,T,,B0010x2y,,y, BThe closedform integration of and in Eqns. (10.1 to 10.3) are given as follows: Byx ,, (A-3) Bds,0s02xy0,x 16 Bds,,,00s0x2y (A-4) y, with: 12 (A-5) x,xS,sin,SPABAB2 12 (A-6) y,yS，cosS,PABAB2 Where is the distance between point A and B. SAB References [1] I. Demirdzic, S. Muzaferija. Finite volume method for stress analysis in complex domains. International Journal for Numerical Methods in Engineering 1994; 37: 3751- 3766. [2] C. Bailey, M. Cross. A finite volume procedure to solve elastic solid mechanics problems in three dimensions on an unstructured mesh. International Journal for Numerical Methods in Engineering 1995; 38: 1757-1776. [3] M.A. Wheel. A mixed finite volume formulation for determining the small strain deformation of incompressible materials. International Journal for Numerical Methods in Engineering 1999; 44: 1843-1861. [4] M.A. Wheel. A finite volume approach to the stress analysis of pressurized axisymmetric structures. International Journal of Pressure Vessels and Piping 1996; 68: 311-317. [5] K. Zarrabi, A. Basu. A finite volume element formulation for solution of elastic axisymmetric pressurized components. International Journal of Pressure Vessels and Piping, 2000; 77: 479-484. 17 [6] J. Fainberg, H.-J. Leister. Finite volume miltigrid solver for thermo-elastic stress analysis in anisotropic materials. Computer Methods in Applied Mechanics and Engineering 1996; 137: 167-174. [7] A.K. Slone, C. Bailey, M. Cross. Dynamic solid mechanics using finite volume methods. Applied Mathematical Modelling 2003; 27: 69-87. [8] M.A. Wheel. A finite volume method for analysing the bending deformation of thick and thin plates. Computer Methods in Applied Mechanics and Engineering 1997; Issues 1-2, 147: 199-208. [9] N. Fallah. A cell vertex and cell centred finite volume method for plate bending analysis. Computer Methods in Applied Mechanics and Engineering 2004; 193: 3457- 3470. [10] K. Maneeratana, A. Ivankovic. Finite volume method for geometrically nonlinear th Annual ACME Conference, School of stress analysis applications. Proceeding of the 7 Engineering, University of Durham, 29-30 March, 1999. [11] W. Pan, M. A. Wheel. A finite volume method for predicting finite strain deformations in incompressible materials. European Conference on Computational Mechanics Solids, Structures and Coupled Problems in Engineering, Munich, Germany. August 31– September 3, 1999. [12] G Taylor, C. Bailey, M. Cross. A vertex based finite volume method applied to non-linear material problems in computational solid mechanics. International Journal for Numerical Methods in Engineering 2003; 56: 507-529. [13] A.B. Oldroyd, M.A. Wheel, T.J. Scanlon. An integrated finite volume based approach for analysing flow induced vibrations. Proceeding of European Conference on Computational Mechanics, Munich, Germany, August 31– September 3, 1999. 18 [14] Slone, A.K., Bailey, C., Pericleous, K.A., Cross, M. Dynamic fluid structure interaction using finite volume unstructured mesh procedures. Computers and Structures, 2002; 80: 371-390. [15] D. J. Allman. A compatible triangular element include vertex rotations for plane elasticity analysis. Computer & Structures 1984; 19: 1-8. [16] D. J. Allman. Evaluation of the constant strain triangle with drilling rotations. . International Journal for Numerical Methods in Engineering 1988; 26: 2645-2655. [17] R.D. Cook. On the Allman triangular and a related quadrilateral element. Computers & Structures 1986; 22: 1065-1067. [18] E. Onate, M. Cervera, O.C. Zienkiewicz. A finite volume format for structural mechanics. International Journal for Numerical Methods in Engineering 1994; 37: 181-201. [19] W. Pan, M. A. Wheel. A finite volume method for solid mechanics incorporating rotational degrees of freedom. Computers & Structures, 2003; Issue 5, 81: 321-329. [20] G.I. Tsamasphyros, C.D. Vrettos. Higher order finite volume formulation for the solution of classical 1D and 2D elasticity. Submitted to Computational Mechanics, 2008. 19 Table 1 Comparison of the displacement, strain and stress between the calculated and exact results under the simple tension condition. UV , , ,,,,D D xxDyyDxyDxxDyyDxyD-4-4-4-4-4(*10 m) (*10 m) (*10 ) (*10 ) (*10 ) (MPa) (Mpa) (MPa) FV 0.0750 -0.0060 0.1000 -0.0300 0.0000 1.000 0.000 0.000 Exact 0.075 -0.006 0.1 -0.3 0.0 1.0 0.0 0.0 Table 2 Comparison of the displacement, strain and stress between the calculated and the exact results under shearing conditions. UV , , ,,,,D D xxDyyDxyDxxDyyDxyD-4-4-4-4-4(*10 m) (*10 m) (*10 ) (*10 ) (*10 ) (MPa) (Mpa) (MPa) FV 0.0000 0.1950 0.0000 0.0000 0.2600 0.000 0.000 1.000 Exact 0.0 0.195 0.0 0.0 0.26 0.0 0.0 1.0 20 Table 3 The convergence test of a square plate subject to pure bending loading by the present method and also a comparison between the present method and two-field and three-field methods. The displacement shown in this table is at point C. Two-field Three- Two-field Three- Mesh (u, v) field Current (u, v) field Current --(*10(*10size Umethod Vmethod (u, v, ,) (u, v, ,) CC33----m) (*10m) (*10(nx*ny) UVU(*10V(*10CCCC3333from[19] m) from[19] m) m) m) from[19] from[19] 4*4 -0.88556 -0.986084 0.993694 -0.95784 -1.06055 -1.06811 8*8 -0.96310 -0.995865 0.998107 -1.03924 -1.07088 -1.07306 12*12 -0.98160 -0.997991 0.999076 -1.05772 -1.07305 -1.07408 16*16 -0.98886 -0.998802 0.999447 -1.06475 -1.07385 -1.07446 24*24 --- --- 0.999733 --- --- -1.07475 Exact 1.0 1.0 1.0 -1.075 -1.075 -1.075 21 Table 4 Comparison of the results between the three FV methods: the present method, and the two-field and three-field methods [19]. U U V V V Mesh UCCCCCC (u, v) (u, v) size (u, v, ,) (u, v, ,) (u, v, ,) (u, v, ,) 3-node 6-node 3-node 6-node triangle triangle triangle triangle 8×2 -0.373781 -0.765342 -0.765342 -2.986773 -6.12700 -6.12700 (54) (81) (81) (54) (81) (81) 16×4 -0.608709 -0.790950 -0.816499 -4.886667 -6.33638 -6.54391 (170) (255) (255) (255) (255) (255) Exact -0.8000 -0.8000 -0.8000 -6.4075 -6.4075 -6.4075 22 Table 5 Comparison of exact solution for and end loaded rectangular cantilever with predictions provided by various FE and FV procedures using three different mesh configurations (number of degrees of freedom in parenthesis) [15,19] Mesh CST Allman LST FV FV Current (Ref. 15,19) Triangle (u, v) method (u, v, ,) (u, v, ,) I 0.0909 0.2696 0.0903 0.3226 (20) (30) ------- (20) (30) ------- II 0.1988 0.3261 0.3550 0.1985 0.3431 0.3483 (54) (81) (170) (54) (81) (81) III 0.3115 0.3471 0.3556 0.3056 0.3520 0.3561 (170) (255) (594) (170) (255) (255) Exact 0.3558 0.3558 0.3558 0.3558 0.3558 0.3558 23 v 3 y 3 v u 5 3 v 6 5 u 5 v 2 6 u6 v 4 6 2 u v 2 1 4 u 4 1 u 1 x O Fig. 1. A two-dimensional 6-node triangular element with a degree of translation only. 24 Edge node of the Centre of the Midside node triangle element triangle element B A C C A (b) (a) Fig. 2. Discritization of a sample model with six-node triangle elements and CVs associated with typical nodes. (a) Conventional FV discritization, ‘A’ edge node(inside the model), ‘B’ midside node (inside the model), ‘C’ edge node(on the boundary of the model). (b) New FV discritization, ‘A’ edge node (inside the model), ‘C’ edge node (on the boundary of the model). 25 y v 3 3 ,u 3 3 v 2 2 u v 2 1 , 2 1 ,u 11 x Fig. 3. A two-dimensional 3-node triangular element with degrees of translation and rotation. 26 Edge node ,yy N Centre of the element B ,yx , P ,xx ,xy S O A A O (b) (a) Fig. 4. A typical FV control volume. (a) control volume of an edge node, (b) the stresses on the one side of the profile. 27 y N , 0B S AC , C (x,y)CC (x,y)AAA x O Fig. 5. Uniformly distributed stresses on the boundary. 28 ,,1.0MPa B (4) (1) (3) y D (2) x A C Fig. 6. Test 1: a quarter model of square plate subject to uniform tensile stress. 29 ,,1.0MPa B (4) (1) (3) y D (2) x A Fig. 7. Test 2: The simple shear of a square plate. 30 A C ,,1.0MPa . B Fig. 8. Test 3: A square plate under pure bending. 31 1.01 05101520 1.00 Exact -1.07 0.99Mesh 16*16m) m)-3Mesh 8*8-3 0.98-1.09Exact Mesh 4*4 U ( *10Mesh 16*16 V ( *100.97 Mesh 8*8-1.11 Mesh 4*40.96 0.95-1.13 05101520-log10(0.5-v) -log(0.5-v)10 (a) (b) Fig. 9. Top right corner node displacements against logarithmic () ,log(0.5,,)10variation in Poisson’s ratio: (a) horizontal displacement; (b) vertical displacement. 32 6E = 3010psi , , ,,,, 3W = 4010lbs thickness = 1.0 inch y 6E = 30×10psi x H = 12 inch , , ,,,, y L = 48 inch x Fig. 10. End loaded cantilever test problem. 33