.. _3dsolidtutorial: 3D Solid Tutorial with Coreform Cubit\ |reg|\ ============================================== Introduction ---------------- This tutorial gives an introduction to the usage of |FOURC| for simulating plastic material behaviour. It is assumed that |FOURC| has been built on your machine according to the instructions and has passed the tests without error messages. In addition, the pre-processing of the finite element model is built with Coreform Cubit, but the input file can also be generated without this software. The analyzed structure is a simple round tensile bar, which is commonly used to identify the material's stress-strain curve by comparing the simulation with the respective experiment. The geometry, which has a cross sectional diameter of 10 mm and a length of constant diameter of 50 mm is shown in the following figure: .. figure:: /_assets/tut_solid-geometry.jpg :alt: Problem definition and geometrical setup :width: 500px :align: center Problem definition and geometrical setup The material under consideration is an Aluminum alloy with a yield strength of 330 MPa and a mild isotropic hardening. During the test, a force *F* acts on the top of the rod such that it elongates. For elasto-plastic simulations |FOURC| provides mainly the following seven material models: :: MAT_Struct_PlasticLinElast MAT_Struct_PlasticNlnLogNeoHooke MAT_Struct_ThermoPlasticHyperElast MAT_Struct_ThermoPlasticLinElast MAT_Struct_Damage MAT_Struct_PlasticGTN MAT_Struct_DruckerPrager The following flow chart depicts for which cases these seven models are applicable: .. figure:: /_assets/tut_solid_flow_chart.jpg :alt: plasticity models flow chart :width: 800px :align: center Capabilities of plasticity models available in 4C Preprocessing ---------------- Due to the symmetry of this model, we may restrict our analysis to one half in longitudinal direction. In addition, one can see that the structure has a rotational symmetry. Since axisymmetric elements are not available in |FOURC|, a 3D structure is modelled with double symmetry in both lateral directions, so that only one eighth of the structure is left. The final FE model is shown here. .. figure:: /_assets/tut_solid-mesh.jpg :alt: tutorial_solid mesh :width: 500px :align: center Mesh of one eighth of the tensile bar Symmetry boundary conditions are applied to the three symmetry planes, denoted with ``XSYMM, YSYMM, ZSYMM`` in the above figure. The top surface is subject to a Dirichlet boundary condition with linearly increasing displacement in y-direction. To this end, we model 200 time steps spanning from :math:`t_0 = 0` s to :math:`t_1 = 1` s with a time step of :math:`\Delta t = 0.005` s. Assuming nonlinear kinematics (large strain theory) and plasticity with isotropic hardening following an exponential curve, namely :math:`\sigma_Y = \sigma_{Y,0} + (\sigma_{Y,\infty} - \sigma_{Y_0}) \left[ 1 - \exp \left( - k \, \varepsilon_p \right) \right]` with :math:`\sigma_{Y,0}=330, \sigma_{Y,\infty} = 1000, k=5` . We choose the model ``MAT_Struct_PlasticNlnLogNeoHooke`` to accomplish this task, which is a plasticity model that uses von Mises yield criterion .. math:: \Phi = \tilde{\mathbf{\tau}} - \sqrt{\frac{2}{3}} \sigma_Y with :math:`\tilde{\mathbf{\tau}}` being the deviatoric Kirchhoff stress. Additionally, a compressible Neo-Hooke elasticity model expressed by the free energy potential .. math:: \rho_0 \Psi(\mathbf{B}^e) = \frac{K}{2} \left[ \frac{1}{2} (J^2 - 1) - \ln J \right] + \frac{1}{2} \mu \left[ \text{tr} \mathbf{\tilde{B}}^e - 3 \right] Here, :math:`J` is the determinant of the deformation gradient :math:`\mathbf{F}`, :math:`\mathbf{B}^e = \mathbf{F}^e {\mathbf{F}^e}^T` is the elastic part of the left Cauchy Green tensor, and :math:`K, \mu` are elastic constants. The geometry, finite element mesh and node sets for the boundary conditions are created using Coreform Cubit\ |reg|\ . The mesh is rather fine and 27-node brick elements have been used herein. The corresponding journal file can be run to reproduce this mesh and output a binary EXODUS mesh information file ``tutorial_solid.e``. The geometry dimensions (actually mainly the radius, all other dimensions depend on it) and element size can be varied, since they are parametrized in the journal file. .. literalinclude:: /tutorial_solid.jou :linenos: The general definition of the geometry in this cubit journal file, which can be run from the directory ``/tests/tutorials/`` includes the complete mesh and four surface node sets. If you don't have Coreform Cubit, you may use the resulting EXODUS file directly, which is also located in this directory. After the file containing the geometry, which is binary, has been created, another input file is needed, which contains - simulation parameters, - boundary conditions, and - further information on the geometry, e.g., the material assignment for the elements. This file should be written in yaml format, here it is called ``tutorial_solid.4C.yaml``. First of all, it contains the problem type, which is a 3D solid mechanics problem here: .. code-block:: yaml PROBLEM TYPE: PROBLEMTYPE: Structure PROBLEM SIZE: DIM: 3 Then we need the information about the geometry to be read in; this is done here: .. code-block:: yaml STRUCTURE GEOMETRY: FILE: tutorial_solid_geo.e ELEMENT_BLOCKS: - ID: 1 ELEMENT_NAME: SOLID ELEMENT_DATA: MAT 1 KINEM nonlinear The element type is a quadratic (27 node) 3D solid element(Solid HEX27) with nonlinear kinematics, that is, including a large strain formulation. Since the block number and material number can be different, the material, kinematics (small or large deformation) and other 4C-specific data are to be entered after the keyword ``ELEMENT_DATA``. The conditions section defines all conditions that are applied to the structure. In the current structure, three symmetry conditions (based on the surfaces defined in Cubit) and one loading condition have to be defined. They are thus Dirichlet type conditions. Note that only surface conditions are given; the lines that belong to two surfaces automatically get the combined boundary conditions. The loading condition has a function number (1) included for the y-direction, and also gets a tag: ``TAG monitor_reaction``. The function is defined below the condition with the identifier (``FUNCT1``), to which the Dirichlet boundary condition refers. The tag is for an extra output of the total forces, which are printed in a csv file, if the output is requested in the header file, see below. .. code-block:: yaml DESIGN SURF DIRICH CONDITIONS: - E: 1 NUMDOF: 3 ENTITY_TYPE: node_set_id ONOFF: - 0 - 1 - 0 VAL: - 0 - 0 - 0 FUNCT: - 0 - 0 - 0 - E: 2 NUMDOF: 3 ENTITY_TYPE: node_set_id ONOFF: - 0 - 0 - 1 VAL: - 0 - 0 - 0 FUNCT: - 0 - 0 - 0 - E: 3 NUMDOF: 3 ENTITY_TYPE: node_set_id ONOFF: - 1 - 0 - 0 VAL: - 0 - 0 - 0 FUNCT: - 0 - 0 - 0 - E: 4 NUMDOF: 3 ENTITY_TYPE: node_set_id ONOFF: - 0 - 1 - 0 VAL: - 0 - 1 - 0 FUNCT: - 0 - 1 - 0 TAG: monitor_reaction FUNCT1: - COMPONENT: 0 SYMBOLIC_FUNCTION_OF_SPACE_TIME: pull - VARIABLE: 0 NAME: pull TYPE: linearinterpolation NUMPOINTS: 2 TIMES: - 0 - 20 VALUES: - 0 - 20 After the boundary conditions, additional parameters are defined for the material, solvers, and output information. Particularly, the output of the forces is requested by a particular section, ``IO/MONITOR STRUCTURE DBC``, which controls accuracy and frequency of the output. This and all other output information is given in section ``IO`` and its subsections. .. code-block:: yaml IO: STRUCT_STRESS: Cauchy STRUCT_STRAIN: GL IO/RUNTIME VTK OUTPUT: INTERVAL_STEPS: 1 IO/RUNTIME VTK OUTPUT/STRUCTURE: OUTPUT_STRUCTURE: true DISPLACEMENT: true STRESS_STRAIN: true GAUSS_POINT_DATA_OUTPUT_TYPE: nodes IO/MONITOR STRUCTURE DBC: INTERVAL_STEPS: 1 You'll find two solver sections in this file as well, which are selected in the section ``STRUCTURAL DYNAMIC`` by the parameter ``LINEAR_SOLVER``. The value of this parameter refers to the section ``SOLVER x``, in which ``x`` is the number of the selected solver. In the section ``SOLVER 1``, the direct solver is defined, it does not have any parameters (except an optional name). The iterative solver in section ``SOLVER 2`` has some parameters, namely an xml file, which contains more parameters. The given xml file as well as some other examples for different multiphysics problems are given in the folder ``/tests/input-files/xml/*/*.xml``. .. code-block:: yaml STRUCTURAL DYNAMIC: INT_STRATEGY: Standard DYNAMICTYPE: Statics RESTARTEVERY: 1000 NUMSTEP: 10 MAXTIME: 10 TOLDISP: 0.0001 TOLRES: 0.0001 NORM_RESF: Rel PREDICT: TangDis LINEAR_SOLVER: 1 SOLVER 1: SOLVER: Superlu NAME: Direct_Solver SOLVER 2: SOLVER: Belos AZPREC: MueLu AZTOL: 1e-05 AZOUTPUT: 1000 AZSUB: 100 MUELU_XML_FILE: elasticity_template.xml NAME: Iterative_Solver If you choose ``LINEAR_SOLVER: 2``, you'll get the iterative solver. Note that you have to copy the xml file to your folder before executing |FOURC| with the iterative solver setting. The ``STRUCTURAL DYNAMIC`` section also defines the time integration, particularly the time stepping. Time step size is given by ``TIMESTEP``, the maximum time is given by ``MAXTIME`` and the maximum number of increments is ``NUMSTEP``. The end of the simulation is defined either when the maximum time or when the number of increments is reached, whatever comes earlier. Note that the number of time steps is currently set to 10, in order to keep the simulation time short for the CI/CD tests. If you want to run the whole simulation including necking of the specimen, you have to change the parameter to at least ``NUMSTEP: 200``. Besides output and solver details, the materials are defined in the input file. For the single material used here, we choose a large strain plasticity model, which is named ``MAT_Struct_PlasticNlnLogNeoHooke`` with parameters, which might represent a high strength Aluminum. .. code-block:: yaml MATERIALS: - MAT: 1 MAT_Struct_PlasticNlnLogNeoHooke: YOUNG: 70000 NUE: 0.33 DENS: 1 YIELD: 330 SATHARDENING: 1000 HARDEXPO: 5 VISC: 0.01 RATE_DEPENDENCY: 1 Last but not least, a section called ``RESULT DESCRIPTION`` is given in the input file, which is only used for the deployment tests. Some particular values are checked in the automated testing procedure. For the user, this section has no relevance. The yaml file and the geometry file are already contained in the ``/tests/tutorials`` directory. Simulation ----------- After the file ``tutorial_solid.4C.yaml`` has been created, the simulation is run on a single processor by the command :: /4C tutorial_solid.4C.yaml tutorial_solid Since the solver output is quite lengthy, one might want to redirect it into a file. For running the simulation with several processes, use ``mpirun -np /4C [further parameters]``. Post processing ---------------- The VTK results are collected in the directory ``tutorial_solid-vtk-files/``. If they shall be opened by Paraview, the data file ``tutorial_solid-structure.pvd`` contains the complete collection of files, so only this file is to be opened, all other files are then read automatically. Using the filter ``Warp by vector`` with Coloring defined by the scalar ``accumulated_plastic_strain``, one may obtain a contour plot on the deformed geometry: .. figure:: /_assets/tut_solid-deformed.jpg :alt: Paraview post processing output :width: 600px :align: center Deformed mesh with contours of accumulated plastic strain The output of the forces needed to apply the Dirichlet boundary condition are given in the directory ``tutorial_solid_monitor_dbc``. Each boundary condition using the tag ``monitor_reaction`` creates a single ASCII file with suffix ``.csv``, which for the present case provides the time, the force in y-direction and the initial and current area to which the force is applied. Since the displacement itself is not given in this file, the displacement is to be calculated from the time. A force-displacement graph may be created using gnuplot by the lines :: set terminal png size 1200,900 font Arial 14 set output 'tutorial_solid_graph.png' set key off set xlabel "displacement"; set ylabel "force"; set title "Force-displacement curve" plot "tutorial_solid_monitor_dbc/tutorial_solid_10004_monitor_dbc.csv" using ($2*12):(abs($6)) with lines title "force-displacement" .. figure:: /_assets/tutorial_solid_graph.png :alt: result curves for solid tutorial :width: 700px :align: center Force displacement curve for the tensile bar