Plate with a hole under uniaxial tension - GCMLab/GCMLab-FEM GitHub Wiki

This tutorial describes how to pre-process, run, and post-process the simulation of a square plate with a circular hole under uniaxial stress.

• Pre-requisites
• Problem Statement
• Mesh Generation
  • Open GMSH
  • Create Geometry
  • Create Mesh
• Configuration File Setup
  • Mesh Properties
  • Material Properties
  • Boundary Conditions
  • Computation Control
• Running the Simulation
• Visualizing the Results
  • Loading the Results
  • Contours
  • Deformed Configuration
  • Verification

• Installed Matlab with an active license
• Installed GMSH (Download here)
• Installed Paraview (Download here)
• Downloaded GCMLab-FEM on the local machine

^ Back To Top

The solid deformation of a square plate is simulated in this tutorial. The plate has side lengths of 4 m, and contains a circular hole in its center of radius R = 0.5 m. This problem has two symmetry planes along the horizontal and vertical axis. As such, only a quarter of the domain is used in the simulation, shown by the shaded area in the figure below.

A two-dimensional approximation is used for this problem, which is valid since the load is applied along the plane of the plate. Plane stress conditions are assumed since the thickness of the plate is assumed to be much smaller than the other plate dimensions.

Kirsch's solution for stresses around a hole for an infinite plate under uniaxial tension can be used to verify the solution. The full set of equations can be found here. Along the vertical axis of symmetry, the analytical solution for the horizontal stress is given as

The analytical solution for the vertical stress along the vertical axis of symmetry is

The simulation results will be compared with this analytical solution.

^ Back To Top

The mesh is generated using the GMSH software. The following steps are outlined below:

• Open GMSH
• Create Geometry
• Create Mesh

• Under the top navigation, go to File > New
• Save the new geometry file as 'Config Files/PlateWithHole.geo'
• Select the OpenCascade geometry kernel

Tips:

• Hold the left mouse button and drag to change the orientation of the axes
• Hold the right mouse button and drag to pan

^ Back To Top

The navigation menu on the left contains commands for editing the geometry, mesh, and solver. We will first work with the geometry commands to generate the model geometry.

Define the geometry points

• Open up the menu items under Geometry > Elementary entities > Add

• Select Point

• Using the popup menu, set the X, Y, and Z axes as 0, 0, and 0, respectively. Set the prescribed mesh size at point to 0.5. Click Add to add the first point.

• Repeat the process for the following list of points:

  X	Y	Z
  0.5	0	0
  2	0	0
  2	2	0
  0	2	0
  0	0.5	0

• Click q on your keypad to finish adding points.

  

Tip: You can manually edit the point coordinates by selecting Geometry > Edit script and modifying the text file that contains the geometry information. Once the edits have been made, save the text file, close it, and select Geometry > Reload script

Define the geometry boundary

• Select Geometry > Elementary entities > Add > Circle arc to create the boundary around the center of the domain

• Select the start point as Point 2 (0.5, 0, 0) with your mouse

• Select the center point as Point 1 (0, 0, 0) with your mouse

• Select the end point as Point 6 (0, 0.5, 0) with your mouse

• Click q on your keyboard to finish

  Tip: The ESC key turns the mouse selection ON and OFF. Ensure it is ON to be able to select the points with your mouse.

• Select Geometry > Elementary entities > Add > Line to create the remaining boundaries around the domain

• Select the start point as Point 2 (0.5, 0, 0) with your mouse

• Select the end point as Point 3 (2, 0, 0) with your mouse

• Repeat the two steps above with the following start and end points:

  Start	End
  Point 3 (2, 0, 0)	Point 4 (2, 2, 0)
  Point 4 (2, 2, 0)	Point 5 (0, 2, 0)
  Point 5 (0, 2, 0)	Point 6 (0, 0.5, 0)

• Click q to finish adding lines

Delete construction point

• Select Geometry > Elementary entities > Delete

• Select Point 1 (0, 0, 0)

• Click e to end selection and click q to finish adding curves

  

Define the domain surface

• Select Geometry > Elementary entities > Add > Plane surface
• Select the boundary of the domain
• Click e to end selection and click q to finish adding curves

^ Back To Top

Define the number of elements along the boundary of the domain

• Select Mesh > Define > Transfinite > Curve
• Set the Number of Points to 20, and leave the default Type and Parameter
• Select all boundary edges of the domain with your mouse
• Click e to end selection and click q to finish adding curves

Generate mesh

• Select Mesh > 2D to create the mesh

  Note that the element type is triangular by default

  

Export mesh file

• In the top navigation, select File > Export
• Save the file as 'Config Files\PlateWithHole.msh'
• Select the type 'Mesh - Gmsh MSH (*.msh)'
• Click Save
• Select the 'Version 2 ASCII' format
• Leave all check boxes unticked
• Click Ok to export mesh

^ Back To Top

A new configuration file is created to define this problem.

• Make a copy of 'Config Files\MasterConfigFile.m' and save the file as 'Config Files\PlateWithHole.m'.
• Modify the function header so that it reads

function [Mesh, Material, BC, Control] = PlateWithHole(config_dir, progress_on)

The remainder of the Configuration file consists of:

• Mesh Properties
• Material Properties
• Boundary Conditions
• Computation Controls

^ Back To Top

The mesh is contained in the Mesh structure in the program. Define the mesh properties as specified below:

%% Mesh Properties
    MeshType = 'GMSH';        
    meshFileName = 'PlateWithHole.msh';
    nsd = 2;
    Mesh = BuildMesh_GMSH(meshFileName, nsd, config_dir, progress_on);            

This tells the program that the mesh is imported from the GMSH mesh file PlateWithHole.msh and that the simulation has two dimensions. The program thenruns the mesh generation algorithm BuildMesh_GMSH which stores the generated mesh in the Mesh structure.

^ Back To Top

The physical properties are stored in the Material structure. For this problem, we set the mechanical properties of steel given in the table below:

The plate is assumed to follow plane stress conditions and the plate has a unit thickness. The code used to define the material properties in the config file is:

%% Material Properties
    Material.E = @(x) 2e11;   
    Material.nu = @(x) 0.3;
    Material.Dtype = 'PlaneStress'; 
    Material.t = @(x) 1;

^ Back To Top

Boundary conditions are stored in the BC structure.

First, symmetry conditions are applied to the left and bottom edges by setting the horizontal displacements on the left edge and the vertical displacements on the bottom edge to zero. The BuildMesh_GMSH function has already defined some DOF sets that we can use to set up the boundary conditions. Mesh.left_dofx is a vector containing the horizontal degrees of freedom on the left edge of the domain. Similarly, Mesh.bottom_dofy is a vector containing the vertical degrees of freedom on the bottom edge of the domain. The code to specify the Dirichlet boundary conditions is:

%% Boundary Conditions
    % Dirichlet boundary conditions (essential)
        BC.fix_disp_dof = [Mesh.left_dofx; Mesh.bottom_dofy];
        BC.fix_disp_value = zeros(length(BC.fix_disp_dof),1);  

Next, tractions are applied to the right edge of the domain to represent the uniaxial tensile load of 10 kPa on the plate. Again, we take advantage of the pre-defined set Mesh.right_nodes to specify the nodes on the right edge along which the traction load is applied. The code for defining the forces on the right edge is given below:

    % Neumann BC
        % uniform tensile stress applied to right edge
        t = 10e3;                               
        % length of right edge
        H = max(Mesh.x(:,2));                   
        % number of tributary lengths along the edge
        num_h = length(Mesh.right_nodes) - 1;   
        % load applied to each node 
        Fnode = t*H/num_h;
        
        % find the nodes in the top right and bottom right corners
        toprightnode = find(Mesh.x(Mesh.right_nodes,2) == max(Mesh.x(:,2)));
        botrightnode = find(Mesh.x(Mesh.right_nodes,2) == min(Mesh.x(:,2)));
        cornernodes = [toprightnode, botrightnode];

        BC.traction_force_node = Mesh.right_nodes;
        BC.traction_force_value = [Fnode*ones(size(Mesh.right_nodes)),... 
                                    zeros(size(Mesh.right_nodes))];
        BC.traction_force_value(cornernodes,1) = Fnode/2;

^ Back To Top

Other controls used for the simulation are set in the Control structure. The controls defined include the number of quadratures used in the integration, the stress recovery method used, and the linear solver type used. The code for the controls is given below:

%% Computation controls
    Control.qo = 2;                         % number of quadrature points
    Control.stress_calc = 'L2projection';   % stress recovery method
    Control.LinearSolver = 'LinearSolver1'; 

^ Back To Top

To run the simulation, first modify RunSim.m. The DirFolder variable should be set to the directory in which the configuration file is stored. The FileList variable should be set to the name of the configuration file. Define the directory where the output will be stored in the VTKFolder variable and turn on output of VTK files and progress indicators by setting plot2vtk and progress_on variables to 1.

%% Input
    DirFolder = 'Config Files';
    FileList = 'PlateWithHole';
    VTKFolder ='C:\Users\USERNAME\Documents\Matlab Results\';    
    plot2vtk = 1;
    progress_on = 1;

Run the function RunSim.m to run the simulation.

^ Back To Top

Once the simulation is complete, a VTK file containing the results will be exported to the directory defined in the variable VTKFolder in RunSim.m. The tutorial will guide you through the following steps:

• Loading the Results
• Contours
• Deformed Configuration
• Verification

^ Back To Top

• Open Paraview
• In the top navigation, select File > Open
• Select the VTK file from the directory: C:\Users\USERNAME\Documents\Matlab Results\Config Files\PlateWithHole\PlateWithHole.vtk
• Select Ok to load the file

Paraview has various toolbars. The Pipeline Browser, Properties, and Information toolbars are displayed by default on the left-hand side of the application screen. The Pipeline Browser shows the VTK file that has been opened. In the Properties toolbar, select Apply to load the data into Paraview.

^ Back To Top

The different data fields can be selected from the Active Variables toolbar located near the top of the application. Change the variable to Sxx to view the horizontal stress contours on the domain.

The mesh can be viewed by changing the domain representation to Surface with Edges in the Representation toolbar. This toolbar is located next to the Active Variables toolbar by default. The results show the mesh on the initial undeformed configuration.

Explore different contour variables and representation views to get acquainted with these settings. Note that vector fields can display the different components or the magnitude separately.

^ Back To Top

• In the representation toolbar, set the view to Wireframe
• In the Active Variables toolbar, set the contour variable to Solid Color to remove the contour

• In the top navigation, select Filters > Common > Warp By Vector. This tool can also be found in the Common toolbar with the icon: . This tool will create a new object, or Pipeline, that is visible in the Pipeline browser. The new object is by default called WarpByVector1
• In the Properties toolbox, under the Properties header, select the Vector U and a Scale Factor of 1e6.
• Select Apply to load the object
• Change the contour variable to U in the Active Variables toolbar to view the displacement contours
• In the Pipeline Browser, you will notice that the WarpByVector1 object is visible (open eye) and the initial object PlateWithHole.vtk is not visible (closed eye). Click on the eye next to the PlateWithHole.vtk object to make it visible. Now both objects should be rendered, showing the deformed and underformed states like the image below:

^ Back To Top

Calculate the analytical stress at each mesh point

• In the top navigation, select Filters > Common > Calculator. This tool can also be found in the Common toolbar with the icon: . This tool will create a new object, called Calculator1, visible in the Pipeline Browser.
• In the Properties toolbox, under the Properties header, change the Results Array Name to Sxx_analytical
• Add the following calculation in the calculation box: 10e3*(1+0.5^2/2/coordsY^2 + 3*0.5^3/2/coordsY^4)
• Select Apply. The Calculator1 object will now have the same fields as the PlateWithHole.vtk object in addition to the field Sxx_analytical
• With the Calculator1 object highlighted, select the Calculator tool again to add another object (Calculator2). The Result Array Name for the new field is Syy_analytical and the calculation is: 10e3*(3*0.5^2/2/coordsY^2 + 3*0.5^3/2/coordsY^4)
• Select Apply

Plot the analytical and numerical solutions

• With the Calculator2 object highlighted in the Pipeline Browser, select Filters > Data Analysis > Plot Over Line from the top navigation. This will create a new object called PlotOverLine1 in the Pipeline Browser

• In the Properties toolbar, under the Properties header, select the following points for the line:

  Point	X	Y	Z
  Point1	0	0.5	0
  Point2	0	2	0

• Select Apply in the Properties toolbar to load the new object. Since the object is a line plot, it will automatically create a new View tab

• In the Properties toolbar, under the Display header, toggle the checkbox for Use Index For XAxis so that it is unchecked. Then select Points_Y from the drop-down menu for X Array Name.

• The Series Parameters header contains a table of all possible fields to display in the vertical axis of the plot. Deselect all variables to reset the plot and select only the following variables: Sxx and Sxx_analytical.

• Explore the options under the View header in the Properties toolbar to modify the plot labels Tip: Selecting the gear icon in the Properties toolbar toggles the advanced settings

• To export the plot, select File > Export Scene... from the top navigation. This tool will allow you to export the plot in a variety of formats, including PDF, EPF, and SVG. A sample plot is shown below:

  

• Repeat the same steps to plot Syy and Syy_analytical

^ Back To Top

1. Increase the mesh resolution
2. Introduce mesh grading
3. Change the plate size

https://www.openfoam.com/documentation/tutorial-guide/tutorialse9.php