A Finite Volume Method for Solids with a Rotational Degree of Freedom Based on the 6-Node Triangle
W. Pan*¹, M.A. Wheel¹ and Y. Qin2
1. Department of Mechanical Engineering, University of Strathclyde, Glasgow, UK.
2. Department of Design, Manufacture Engineering and Management, University of Strathclyde, Glasgow, UK.
* corresponding author, email:
Abstract:A finite-volume (FV) cell vertex method is presented for determining the displacement field for solids exhibiting with incompressibility. The solid is discretized into six-node finite elements and the standard six-node finite-element shape function is employed for each element. Only control volumes around vertex node of the triangular element are considered. For considering the material incompressibility, a constant hydrostatic pressure as an extra unknown variable within each element is assumed. The force equilibrium in two perpendicular directions and one in-plane moment equilibrium equation are derived for each control volume. The volume conservation is satisfied by setting the integration of volumetric strain as zero within each element. By solving the system control equations, the displacements and rotations of the vertex nodes and the hydrostatic pressure for each element can be obtained and then the displacements of the midside nodes can be calculated. The simulation results show that this FV method passes the patch tests and converges to theoretical results under mesh refinement for material behaviour incompressibility.
Keywords: Finite Volume Method, Control volume, Vertex centred method, Rotational degree, Incompressibility.
1. Introduction
Finite element (FE) method has dominated solid and structure analysis for many years, while FV method has been successfully used for fluid dynamic analysis. Over the last fifteen years, some progresseshave been made for employing FV method for solids and structures. Apart from the application of FV method to elastic[1-6], elastic-plastic[7], and plate analysis [8,9] for compressible materials. FV method has also been applied to large deformation compressible and near-incompressible hyperelastic material analysis[10,11].For incompressible material deformation analysis, a three-node and four-node element has been used for analysis 2D elastic media by using the cell central FV method [12], a bi-linear displacement field was assumed within the quadrilateral constructed by the connection of cell central and adjoining edge nodes, and also the constant hydrostatic pressure within each element. Mixed form FV method was also presented for 2D and 3D and also fracture problems with incompressible linear elastic materials[13], the non-physical stresses oscillations near crack tip was eliminated.FV method that included the Allman degree of rotation for 3-node[14] and 6-node FV [15] analysis for 2D elastic problems were presented, some advantages have been demonstrated for some cases. In this paper, a mixed form FV cell vertex method is presented for determining the displacement field for solids exhibiting incompressibility. The similar discretization method with linear strain distribution six-node triangle element as [15] is used. To include the rotation degree of freedom, the midside displacements are replaced by its neighbouring vertex nodes[16], and the constant hydro-static pressure within each element is assumed although there is no difficulty for the assumption of a linear distribution of hydrostatic pressure within each element. The force equilibrium in two perpendicular directions and one in-plane moment equilibrium equation are derived for each control volume. The volume conservation is satisfied by setting the integration of volumetric strain as zero within each element for small strain elastic problem. The total number of the force and moment equilibrium equations plus the number of volume conservation equationsis the same as the number of the total unknown variables due to the special construction of the CV. By solving these equations, the displacements and rotations of the vertex nodes and the hydrostatic pressure for each element can be obtained and then the displacements of the middle-side nodes can be calculated. In section 2, formulation for displacement, strain and stress filed will be given, and also the compensation equation related with material incompressibility.The construction of special CV and the system equations are given in section 3, and the test problems will be followed in section 4, and the final section is the conclusions.
2. Formulation for displacement, strain and stress fields
The six-node triangle element is considered for the discretization of the 2D incompressible material problem. The same formulation for displacement and strainfield within each element as [15] are used in this paper and are summarized as following:
2.1 Displacement fields
The displacement components, u and v, within a six-node triangular element are:
(1)
Where and are coefficients. By applying equations (1) at every node of the triangle, gives:
(2)
Where , andui,vi, are displacement components of vertex nodes (for i=1,2,3) and midside nodes (for i=4,5,6) along x and y direction respectively.
For the introduction of rotation degree of freedom of each vertex nodes, the midside node displacement will be replaced by the neighbouring vertex node displacements and rotations. The relationship can be repressed by the transformation matrix T,
(3)
Where i, (i=1,2,3) are rotations of the vertex nodes of 6-node triangle, andmatrix T can be found in [16].
With the use of relationship between strain and displacement, the direct strain components , and shear strain are
(4)
Where A is a constant matrix, the components of vectorand are constant or linear function of coordinators x or y, and A, and can be found in[15].
By using the constitutive relationship between stresses and strain and the hydrostatic pressure , the stress components are
(5)
where G is the shear modulus, p is hydrostatic pressure. Put (4) into (5), gives
(6)
3. Control volume and the equilibrium equations
3.1 Control volume for cell vertex FV procedure
One difference between FE and FV method is that equilibrium equations must be built on each individual CV for FV, under the cell centre FV procedure, the CV is the same as FE element, while for cell vertex FV procedure, the CV is completely different from FE element. There are some ways to build CV for cell vertex FV procedure.
One method to construct CV for cell vertex FV procedure is to connect the middle point between neighbouring nodes and the centre of the element to form a multi-side complex shape CV. Other method was presented in [17], some position of the points needed to be connected to construct CV was obtained by trial and error method. In this paper, due to the introduction of rotational degree of freedom, special CV as used in [15] was adopted. Figure 1(a) shows the typical CVs(constructed by dotted lines) around a vertex nodeA and B inside and on the boundary of the model respectively.
Assuming the total number of vertex node and element are ‘nd’ and ‘ne’ respectively, thus the total numer of the system unknown variables will be ‘3*nd+ne’. There are two force equilibrium equations along perpendicular directions and one moment equilibrium equation for each CV, so the total number of equation related with forces and moment are ‘3*nd’, and it is less than the total number of unknown variables. Therefore, the compensation equation has to be built. The incompressibility condition for each element can be used for this purpose.
Figure 1Typical CVsand inner stresses. Solid circle, open circle and small triangle represent vertex nodes, midside nodes and geometrical centre of FE element. (a) CVs around the inside vertex node A and boundary vertex node B (b) inner stresses on one of the segments of the profile of CV.
3.2 Force and moment equilibrium equations
To obtain the force and moment equilibrium equations for each CV, the resultant forces and moment have to be integrated along the profile of each CV. From formula (5), it can be seen that stress components within each CV is linearly distributed, so the resultant forces and moment can be integrated analytically along the profile of each CV. For one typical segment of the CV shown in Figure 1(b), the coordinators of any point along the segment can be expressed as
(7)
where the angle between the normal of line CD and the x-axis is and the distance from point C to P is S.
By integrating explicitly the stresses components along the profile of each CV, the resultant force components and acting in x and y directions are.
(8.1)
(8.2)
The resultant moment about the origin of the coordinate system is:
(8.3)
By summarizing all the resultants force and moment along the profile of each CV, the equilibrium equation for each CV can be written as following:
(9)
Where , and are external forces and moment acting on the CV.
3.3 Compensation equation for incompressibility
For linear elastic material, the incompressibility of each element can be expressed as
(10)
For simplicity, the strain components at the centre of the triangular element were used for above formula (10), therefore,
(11)
Each element has one compensation equation and one unknown hydrostatic pressure variable.
3.4 system equations
Combining formulation (9) and (11), the system equation can be assembled as following:
(12)
where
3.5. Boundary conditions
For concentrated point force or moment acting on each CV, they can be directly put into the right-hand side of equations (12). For distributed stresses on boundaries, resultant forces alone x and y direction and moment to the origin of coordinator can be easily integrated along the profile of CV, then the same procedure as concentrated forces and moments can be used. The stiffness matrix of the whole structure shown in equations (12) is a singular matrix, displacement and rotation boundary conditions have to be applied to the equations, the standard procedure as the treatment of displacement boundary conditions for FE can be used. That is, let the diagonal element of the stiffness matrix corresponding to the ithconstraint as one, while the others in the same row as zero, also let the ith element of the right-hand side vector as the exact value of the constraint. By solving the modified Eqn. (12), the vertex displacements, rotations and hydrostatic pressure can be obtained, after which, using Eqn. (3), the midside node displacements can be calculated. From Eqns. (4) and (6), the strain and stress field can be determined subsequently.
4. Numerical test results
Four 2D plane strain benchmark problemswith incompressible material behaviour were considered. The first and the second examples are isotropic plate under simple tension and simple shear, the third one is case of asquare plate subject to pure bending. The fourth case is an infinite long tube subjecte to internal pressure. The value of material’s shear modulus, G, is chosen as 1.0GPa. For each case, the calculated results were compared with that of theoretical results.
4.1 Test 1and 2: simple tension and shear of a square plate
An infinitely long square plate with the length and width of 2m*2m subjected to simple tension stress 1.0MPa is shown by the solid lines in Figure 2. A quarter model is analysed and the model is discretized into four different size 6-node triangle elements. The dotted lines in Figure 2 show the expected deformed mesh. The calculated results show that the constant strain field is exactly predicted. Figure 3 shows the initial and the deformed mesh of an infinitely long square plate with length and width of 1m*1m subject to simple shear. As the simple tension case, the solid and dotted lines represent initial and deformed mesh. From this figure, the exact deformation shape is obtained.
4.2 Test problem 3: Square plate subjected to pure bending
The third test example, shown in Figure 4 (a half model), is an isotropic square plate with length 2.0m subject to pure bending loading and the linearly varying edge normal stress is given by
(13)
Theoretical results of the displacement components are as following
(14)
where Young’s modulus E=2G(1+) is and =0.5.
Figure 4 A thick square plate subject to pure bending
Figure 5 Top right corner vertical displacement changes with the number of division along x direction
Figure 6 Top right corner horizontal displacement changes with the number of division along x direction
The half plate model is tested using a consistent refinement of a regular mesh with the same number of division along x and y directions. The boundary conditions are: The horizontal displacements and in-plane rotations along left edge are zeros; the vertical displacement at the middle point of the left edge is zero. Figure 5 showsthe comparison between the calculated and theoretical results of vertical displacement of top right corner under mesh refinement. From this figure, it can be seen that the calculated results converge to the theoretical results. Figure 6 shows the results of the horizontal displacement at the top right corner, similar trend as Figure 5 can be observed although the horizontal displacement converges slower than that of vertical displacement.
4.3 Test problem 4: Inflation of an infinitely long tube under internal pressure
Figure 7(a) is a schematic drawing of an infinitely long tube(a quarter model) subjected to internal pressure, the tube material is linear incomprensible material with shear modulus 100.0GPa. The analytical solution of radial displacement is
(15)
Figure 7Infinitely long tube subjected to internal pressure. (a) a quarter model with loading and boundary shown, (b) the undeformed(solid line) and deformed mesh(dashed line), solid and empty circle represent nodes attached to undeformed and deformed mesh respectively, the number of division along the radial and hoop directions are all 4, displacement amplification factor is6e4.
Figure 8 The relative error varies with the number of division along radial direction
Figure 9 The comparison between numerical and analytical results
The tube inner and outside radius are 0.1m and 0.4m respectively, the internal pressure is 1.0MPa. The symmetry boundary conditions were applied for the line AB and CD, and the rotations along these two lines are all constrained.Figure 7(b) shows a the undeformed and deformed mesh of a special case with the number of division along radial and hoop direction are all 4. From this figure, the axi-symmetrical deformation is observed. Figure 8 shows the relative error of the inner born radial displacement against the number of divisions along the radial direction. With the mesh refinement, the relative error becomes smaller and smaller. In Figure 9, with the division along radial and hoop direction of 32, the comparison between the numerical and analytical results of radial displacement profile along boundary AB is presented. It can be clearly seen that the presented method can give good results.
5. Conclusions
Six-node triangle element with high order displacement field was used for cell vertex FV procedure to tackle the 2D problem with material behaviour incompressibility. Constant hydrostatic pressure within each element was assumed. With special construction of CV and the replacement of midside node displacement by neighbouring edge node displacement and rotations, the total number of unknown variables equals to the total number of force and moment equilibrium equations. The benchmark test examples shows that the presented FV method can predict the constant strain field exactly, and can also predict the linear and more complex strain field with the mesh refinement
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 finite volume approach to the stress analysis of pressurized axisymmetric structures. International Journal of Pressure Vessels and Piping 1996; 68: 311-317.
[4] 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.
[5] 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.
[6] 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.
[7] 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.
[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 stress analysis applications. Proceeding of the 7th Annual ACME Conference, School of 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] 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.
[13] I. Bijelonja, I. Demirdzic, S. Muzaferija. A finite volume method for incompressible linear elasticity. Computer methods in applied mechanics and engineering, 195(2006), 6378-6390.
[14] 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.
[15] W. Pan, M. A. Wheel. A Six-node Triangle Finite Volume Method for Solids with a Rotational Degree of Freedom. Submitted to Communications in Numerical Methods in Engineering, 2008.
[16]R.D. Cook. On the Allman triangular and a related quadrilateral element. Computers & Structures 1986; 22: 1065-1067.
[17] 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.
Figure captions
Table captions
1