2.2 Stress analysis of a plate with a hole

This tutorial describes how to pre-process, run and post-process a case involving linear-elastic, steady-state stress analysis on a square plate with a circular hole at its centre. The plate dimensions are: side length 4  and radius  0.5 . It is loaded with a uniform traction of  10   over its left and right faces as shown in Figure 2.21. Two symmetry planes can be identified for this geometry and therefore the solution domain need only cover a quarter of the geometry, shown by the shaded area in Figure 2.21.

The problem can be approximated as 2-dimensional since the load is applied in the plane of the plate. In a Cartesian coordinate system there are two possible assumptions to take in regard to the behaviour of the structure in the third dimension: (1) the plane stress condition, in which the stress components acting out of the 2D plane are assumed to be negligible; (2) the plane strain condition, in which the strain components out of the 2D plane are assumed negligible. The plane stress condition is appropriate for solids whose third dimension is thin as in this case; the plane strain condition is applicable for solids where the third dimension is thick.

An analytical solution exists for loading of an infinitely large, thin plate with a circular hole. The solution for the stress normal to the vertical plane of symmetry is

(2.14)

Results from the simulation will be compared with this solution. At the end of the tutorial, the user can: investigate the sensitivity of the solution to mesh resolution and mesh grading; and, increase the size of the plate in comparison to the hole to try to estimate the error in comparing the analytical solution for an infinite plate to the solution of this problem of a finite plate.

2.2.1 Mesh generation

The domain consists of four blocks, some of which have arc-shaped edges. The block structure for the part of the mesh in the plane is shown in Figure 2.22. As already mentioned in subsubsection 2.1.1.1, all geometries are generated in 3 dimensions in OpenFOAM even if the case is to be as a 2 dimensional problem. Therefore a dimension of the block in the direction has to be chosen; here, 0.5  is selected. It does not affect the solution since the traction boundary condition is specified as a stress rather than a force, thereby making the solution independent of the cross-sectional area.

The user should start up FoamX as normal and open the plateHole case in their own $OpenFOAM_RUN/tutorials/stressedFoam directory. Open the blockMesh utility as described previously in subsubsection 2.1.1.1 and press the Edit Dictionary button. Until now, we have only specified straight edges in the geometries of previous tutorials but here we need to specify curved edges. The user should select edges to view how they are specified as shown in Figure 2.23. All 8 curved edges in this example are listed in a table; each can be edited which opens a selection window offering different types of curve, including arc, simpleSpline, polyLine etc., described further in subsection 6.3.1. In this example, all the edges are circular and so can be specified by the arc keyword entry. Editing the arc selection opens a window in which the user specifies the labels of the start and end vertices of the arc and a point vector through which the circular arc passes.

The blocks in this blockMeshDict do not all have the same orientation. As can be seen in Figure 2.22 the direction of block 0 is equivalent to the direction for block 4. This means care must be taken when defining the number and distribution of cells in each block so that the cells match up at the block faces.

5 patches are defined, one for each side of the plate and one for the hole. The blockMeshDict dictionary entries are given below.

The mesh should be generated using blockMesh and can be viewed in dxFoam as described in subsection 2.1.2. It should appear as in Figure 2.24.

Figure 2.24: Mesh of the hole in a plate problem.

2.2.1.1 Boundary and initial conditions

Once the mesh generation is complete, load the mesh into FoamX: remember this is done by clicking the right mouse button on Mesh in the case directory tree and selecting the Read Mesh function. The names of the patches will appear in a dictionary named Patches which should be set as follows: left and down patches are both symmetryPlane; the undefinedFaces are the front and back planes of the 2D geometry and should therefore be declared empty; the other patches are traction boundary conditions, set by traction boundary type.

The Fields must be set as before; here, the displacement U and temperature T. The traction boundary conditions are specified by a linear combination of: (1) a boundary traction vector; (2) a pressure that produces a traction normal to the boundary surface that is defined as negative when pointing out of the surface. The up and hole patches are zero traction so the boundary traction and pressure are set to zero. For the right patch the traction should be   and the pressure should be 0  as shown in Figure 2.25. All the displacement initial conditions should be set to  .

2.2.1.2 Mechanical properties

The physical properties for the case are set in the mechanicalProperties dictionary. For this problem, we need to specify the mechanical properties of steel given in Table 2.1. In the mechanical properties dictionary, the user must also set planeStress to yes.

Property Units Keyword Value Density rho 7854 Young’s modulus E Poisson’s ratio -- nu 0.3 Table 2.1: Mechanical properties for steel

2.2.1.3 Thermal properties

The temperature field variable T is present in the stressedFoam solver since the user may opt to solve a thermal equation that is coupled with the momentum equation through the thermal stresses that are generated. The user specifies at run time whether OpenFOAM should solve the thermal equation by the thermalStress switch in the thermalProperties dictionary. This dictionary also sets the thermal properties for the case, e.g. for steel as listed in Table 2.2.

Property Units Keyword Value Specific heat capacity C 434 Thermal conductivity k 60.5 Thermal expansion coeff. alpha Table 2.2: Thermal properties for steel

In this case we do not want to solve for the thermal equation. Therefore we must set the thermalStress keyword entry to no in the thermalProperties dictionary.

2.2.1.4 Control

As before, the information relating to the control of the solution procedure are read in from the controlDict dictionary. For this case, the startTime is 0 . The time step is not important since this is a steady state case; in this situation it is best to set the time step deltaT to 1 so it simply acts as an iteration counter for the steady-state case. The endTime, set to 100, then acts as a limit on the number of iterations. The writeInterval can be set to .

The controlDict entries should be saved as follows:

2.2.1.5 Discretisation schemes and linear-solver control

Let us turn our attention to the fvSchemes dictionary. Firstly, the problem we are analysing is steady-state so the user should select SteadyState for the time derivatives in timeScheme. This essentially switches off the time derivative terms. Not all solvers work for both steady-state and transient problems but stressedFoam does work, since the base algorithm is the same for both types of simulation.

The momentum equation in linear-elastic stress analysis includes several explicit terms containing the gradient of displacement. The calculations benefit from accurate and smooth evaluation of the gradient. Normally, in the finite volume method the discretisation is based on Gauss’s theorem, as discussed in subsection 4.4.3 and 2.4.6 of the Programmer’s Guide. The Gauss method is sufficiently accurate for most purposes but, in this case, the least squares method will be used. The user should therefore open the fvSchemes dictionary and select leastSquares for the default gradient discretisation scheme.

The fvSolution dictionary controls the linear equation solvers and algorithms used in the solution. The user should first look at the solvers subdictionary and select the ICCG solver for U. The solver tolerance (the scalar following the solver name) should be set to for this problem. The solver relative tolerance (the following scalar) sets the required reduction in the residuals within each iteration. It is uneconomical to set a high relative tolerance within each iteration since a lot of terms are explicit and are updated as part of the segregated iterative procedure. Therefore a reasonable value for the relative tolerance is , or possibly even higher, say .

The fvSolution dictionary contains a subdictionary, stressedFoam that contains some control parameters specific to the application solver. Firstly there is nCorrectors which specifies the number of outer loops around the complete system of equations, including traction boundary conditions within each time step. Since this problem is steady-state, we are performing a set of iterations towards a converged solution with the ’time step’ acting as an iteration counter. We can therefore set nCorrectors to 1.

The U keyword specifies a convergence tolerance for the outer iteration loop, i.e. sets a level of initial residual below which solving will cease. It should be set to the desired solver tolerance specified earlier, for this problem.

2.2.2 Running the code

The user should run the code here in the background, either from FoamX or from the command line as specified below, so he/she can look at convergence information in the log file afterwards.

2.2.3 Post-processing

Post processing can be performed as in subsection 2.1.4. The stressedFoam solver decomposes the stress tensor field into individual scalar field components, , etc. for purposes of post-processing. The viewer can therefore select the stresses from the volScalarField menu of the Field 1 visualisation window in dxFoam. It is perhaps worth reiterating that variables in the OpenFOAM solvers are usually represented by the mathematical symbol by which they are represented, In the case of Greek symbols, the variable is named phonetically, e.g. is named sigmaxx.

Figure 2.26: stress field in the plate with hole.

We would like to compare the analytical solution of Equation 2.14 to our solution. We therefore must output a set of data of along the left edge symmetry plane of our domain. The user should therefore open the Graph window from Window Control. The graph line can be set as described in subsubsection 2.1.5.7, e.g. here we could set the Grid Point to be and the Vector to be . The user should check the Do Plot and Reveal Line boxes and, in order to see the graph line, may need to turn the opacity of all the mesh components to 0 in the Options window. The user should see the stress fall from approximately 3  the applied load to the applied load along the -axis.

In order to compare the results with the analytical solution, the user must write the graph data to file by going to the Graph window and clicking the Write Now box in the Data Writing panel. The data is written to a file which the user can locate on his/her machine. Note that using the default settings, the file will be written to the working directory from where dxFoam was launched on the command line. The data is written as a simple text file that the user can load into a graph drawing package and simply compare with the analytical solution of Equation 2.14.

The user may also generate the required graph data using the sample utility. In the FoamX case server, the user should select select sample from the postProcessing -> miscellaneous sub-menu of the from the Foam Utilities menu. The user should select Edit Dictionary which opens a sampleDict dictionary, pre-specified in the tutorial case. The dictionary contains the entries as follows and summarised in Table 7.7. The sample line specified in sampleSets is set between and . Clicking Start Run will execute the sample utility.

The writeFormat is raw 2 column format. In an application such as GnuPlot, one could type the following at the command prompt would be sufficient to plot both the numerical data and analytical solution:

Figure 2.27: Normal stress along the vertical symmetry

2.2.4 Exercises

The user may wish to experiment with stressedFoam by trying the following exercises:

2.2.4.1 Increasing mesh resolution

Increase the mesh resolution in each of the and directions. Use mapFields to map the final coarse mesh results from subsection 2.2.3 to the initial conditions for the fine mesh.

2.2.4.2 Introducing mesh grading

Grade the mesh so that the cells near the hole are finer than those away from the hole. Design the mesh so that the ratio of sizes between adjacent cells is no more than 1.1 and so that the ratio of cell sizes between blocks is similar to the ratios within blocks. Mesh grading is described in subsection 2.1.6. Again use mapFields to map the final coarse mesh results from subsection 2.2.3 to the initial conditions for the graded mesh. Compare the results with those from the analytical solution and previous calculations. Can this solution be improved upon using the same number of cells with a different solution?

2.2.4.3 Changing the plate size

The analytical solution is for an infinitely large plate with a finite sized hole in it. Therefore this solution is not completely accurate for a finite sized plate. To estimate the error, increase the plate size while maintaining the hole size at the same value.