Home

Akantu User's Guide - Computational Solid Mechanics Laboratory

1. 66 B 5 Standard linear solid 66 B 6 Elasto plastic linear isotropic hardening 66 B Z Damage Marigo nn a en LO x CIR RH on oe OR Dro e Dod hac de je 66 BA amare Mazars aeq xo bee Ee co EXPO Na d EN Qasaqa deb dg AA am 67 BY ICohesive Wheat is essaie kr bale aed an POP nee ew tuwa pE quc EX Vd eg ce 67 B 10 Rohe mnm 67 B 11 Cohesive exponentiall 67 C Package dependencies 69 73 5 Chapter 1 Introduction Akantu means little element in Kinyarwanda a Bantu language From now on itis also an open source object oriented Finite Element library with the ambition to be generic and efficient Akantu is developed within the LSMS Computational Solid Mechanics Laboratory at the Ecole Polytechnique Federale of Lausanne Switzerland The open source philosophy is important for any scientific software project evolution The collaboration permitted by shared codes enforces sanity when users and not only developers can scrutinize and possibly criticize the implementation details Akantu was born with the vision to associate genericity robustness and efficiency while ben efiting from the open source visibility Genericity is necessary to allow the easy exploration of mathematical formulations through algorithmic ideas Robustness and reliability is
2. loop over nodes in contact surface to create contact elements cout lt lt Adding slave master pairs lt lt endl for auto nit cs node_begin nit cs node_end nit point type p amp coords nit ignore slave node if it doesn t lie within the bounding box if bb p continue int el a search amp coords nit if el 1 cd addPair nit master type model triangle 3 el We then start the loop over displacement increments and at each step we call solveContactStep and save the displacement resulting force and maximum pressure in an array that will be used to print the results at the end of the simulation const size_t steps 30 Real data 3 steps store results for printing Real step 0 001 top displacement increment size_t k 0 for Real delta 0 delta lt step steps delta step apply displacement to the top surface of the half sphere model applyBC BC Dirichlet FixedValue delta BC y top surface solve contact step this method also dumps Paraview files cd solveContactStep amp a data 0 k delta data 1 k cd getForce data 2 k cd getMaxPressure k The last portion of the output of code above is as follows Disp Force Max pressure 0 0 0 0 001 1 29 6 068e 04 58 002 003 004 005 006 007 008 009 0 01 011 012 013 014 015 016 017 018 019 0 02 021 022 023 0
3. ElementTypeMapArray lt Ulnt gt connectivities mesh getConnectivities Then the specific array associated to a given element type can be obtained by 6 CHAPTER 2 GETTING STARTED Array lt UInt gt amp connectivity_triangle connectivities _triangle_3 where the first order 3 node triangular element was used in the presented piece of code 2 5 1 Vector amp Matrix The Array iterators as presented in the previous section can be shaped as Vector or Matrix This objects represent 1 and 2 order tensors As such they come with some functionalities that we will present a bit more into detail in this here Vector lt T gt 1 Accessors o v i givesthe j component of the vector v o v i gives the j component of the vector v o v size gives the number of component 2 Level 1 results are scalars o v norm returns the geometrical norm L2 o v norm lt N gt returns the Ly norm defined as Xi Iva IN e N can take any positive integer value There are also some particular values for the most commonly used norms L 1 for the Manhattan norm L 2 for the geometrical norm and L_inf for the norm infinity o v dot x return the dot product of v and x o v distance x return the geometrical norm of v x 3 Level 2 results are vectors o V S V S v S v s those are element wise operators that sum sub stract multiply or divide all the component of v by the scalar s O V X V x sums Or substracts
4. addDumpfField 49 addDumpFieldExternal addDumpFieldTensor addDumpFieldVector Material 26 create a new material 35 Elastic 27 Neohookean 29 registerParam B8 Small deformation Plasticity Packages AKANTU CPPARRAY 73 AKANTU IMPLICIT 70 AKANTU LAPACK 70 AKANTU MPI 70 AKANTU MUMPS 70 AKANTU NLOPT 70 AKANTU OPTIMIZATION 70 AKANTU SCOTCH 70 71 SolidMechanicsModel initFull setTimeStep solveStep 21 25 StableTimeStep StructuralMechanicsModel assembleStiffnessMatrix solve up dateResidual 43 References 1 C Geuzaine and J F Remacle Gmsh A 3 D finite element mesh generator with built in pre and post processing facilities International Journal for Numerical Methods in Engineering 79 11 1309 1331 2009 2 Simulia ABAQUS FEA 3 TNO DIANA 4 M Ortiz and A Pandolfi Finite deformation irreversible cohesive elements for three dimensional crack propagation analysis International Journal for Numerical Methods in Engineering IJNME 44 1267 1282 1999 5 A Curnier A Curnier and A Curnier M thodes num riques en m canique des solides EPFL 1992 6 T Hughes and T Belytschko A precis of developments in computational methods for transient analysis Journal of Applied Mechanics ASME 50 1033 1041 1983 7 B M Ted Belytschko Wing Kam Liu Nonlinear Finite Elements for Continua and Structures Wiley 2000 8 J Simo
5. o oi 2GAe Atr Ae I 4 53 2 Check the Yielding criteria f Go gir a e 4 54 3 Compute the Plastic multiplier t k k _ Off 3GAP Oy Ap Ap dAp 4 56 o oy hAp 4 57 4 Compute the plastic strain increment tr Ag zap 4 58 e 5 Compute the stress increment Ao 2G Ae amp Ae Atr Ae A amp I 4 59 6 Update the variables gP CA AE 4 60 o 0 AG 4 61 We use an implicit integration technique called the radial return method to obtain the plastic multiplier This method has the advantage of being unconditionally stable however the accuracy remains dependent on the step size The plastic parameters to indicate in the material file are Cy Yield stress and h Hardening modulus In addition the elastic parameters need to be defined as previously mentioned E Young s modulus nu Poisson s ratio 4 4 CONSTITUTIVE LAWS 33 4 4 5 Damage In the simplified case of a linear elastic and brittle material isotropic damage can be represented by a scalar variable d which varies from 0 to 1 for no damage to fully broken material respectively The stress strain relationship then becomes o 1 d C e where o are the Cauchy stress and strain tensors and C is the elastic stiffness tensor This formulation relies on the definition of an evolution law for the damage variable In Akantu many possibilities exist and they are listed below M
6. 12 In Akantu is defaults to 1 2 see section 4 3 The initialization of the simulation is similar to the static and implicit dynamic version The model is created from the SolidMechanicsModel class In the initialization the explicit scheme is selected using the explicit lumped mass constant SolidMechanicsModel model mesh model initFull SolidMechanicsModelOptions explicit lumped mass Note Writing model initFull or model initFull SolidMechanicsModelOptions is equivalent to use the explicit lumped mass keyword as this is the default case The explicit time integration scheme implemented in Akantu uses a lumped mass matrix M reducing the computational cost This matrix is assembled by distributing the mass of each element onto its nodes The resulting M is therefore a diagonal matrix stored in the mass vector of the model The explicit integration scheme is conditionally stable The time step has to be smaller than the stable time step which is obtained in Akantu as follows 4 3 DYNAMIC METHODS 25 time_step model getStableTimeStep The stable time step is defined as Atos Ax as n 4 21 where Ax is a characteristic length e g the inradius in the case of linear triangle element u and A are the first and second Lame s coefficients and p is the density The stable time step corresponds to the time the fastest wave the compressive wave needs to travel the characteristic length of the mesh Howe
7. 4 4 6 Parameters o Sigma c Real Critical stress 0 beta Real B parameter G_cI Real Mode I fracture energy G cII Real Mode II fracture energy kappa Real x parameter penalty Real penalty coefficient a O O O O O B 10 Cohesive bilinear Keyword cohesive_bilinear Description here 4 4 6 Parameters o Sigma c Real Critical stress 0 beta Real B parameter G cI Real Mode I fracture energy G cII Real Mode II fracture energy kappa Real x parameter penalty Real Penalty coefficient a delta 0 Real Elastic limit displacement o O O 00 0 0 B 11 Cohesive exponential Keyword cohesive_exponential Description here 4 4 6 Parameters o Sigma c Real Critical stress 0 o beta Real B parameter o delta c Real Critical displacement c 67 Appendix C Package dependencies During the configuration cmake offers several Akantu options which have dependencies with each other or with external packages and software Each of these are described in details now AKANTU_CORE This package is the core engine of Akantu It depends on o A C compiler GCC gt 4 or Intel o The cross platform open source CMake build system o The Boost C portable libraries o The zlib compression library Under Ubuntu 14 04 LTS the installation can be performed using the commands gt sudo apt get install build essential cmake curses gui libboost dev zlibig dev Under Mac OS X the ins
8. CETERIS IO TTE 512 Setting Boundary Conditions 52 111 iv CONTENTS 6 Heat Transfer Model 45 6 P heor oe eae E ene hue aney ase eceweraseevausees 45 6 2 Using the Heat Transfer Model 45 6 2 1 Explicit DYNAMIC se s essre ern 46 7 49 49 50 8 51 51 52 Distributing Mesh Partitions 52 Launching a Parallel Program 52 9 53 Implicit Contact SOIVEN cuu gea pega EA PROC OH rel 53 9 1 1 Implementationl esee cl 53 A 59 1D Shape Functions uy uoces A OR RE ae sn ant 59 A IL QT mail 2b ac SSS ee bad Pd ERSTE I DE 59 1 2 Seoment 3b eue ee een 59 A 2 2D Shape Functions 60 A ATL LS Ol NeMEMC 1 60 ADD Manli dead E ale ad es e oor de er ek wee eS 60 ADS Quadrangle AL dei ia RE tee ee een 61 Ag A Ouadranele lis d Rr Ra 61 A 3 BD Shape Functions 62 T Tetrahedron s a L ol tea gal Ret a a oan Rx E ee ea 62 Tetrahedron OU 62 A 3 3 Hexahedion 8 u ler ab ae dh idos 63 B Material parameters 65 BE ee eee ee 65 B 2 Linear elastic anisotropicl 65 B 3 Linear elastic orthotropic 65 B 4 Neohookean finite strains
9. Neumann FromStress surface stress Top If the boundary conditions need to be removed during the simulation a functor is called from the Neumann boundary condition to free those boundary conditions from the desired boundary model applyBC BC Neumann FreeBoundary Side User specified functors can also be implemented A full example for setting both initial and boundary conditions can be found in examples boundary conditions cc The problem solved in this example is shown in Fig It consists of a plate that is fixed with movable supports on the left and bottom side On the right side a traction which increases linearly with the number of time steps is applied The initial displacement and velocity in x direction at all free nodes is zero and two respectively a i BAA Figure 4 3 Plate on movable supports 4 1 3 Material Selector If the user wants to assign different materials to the different finite elements of choice in Akantu a material selector has to be used By default Akantu assigns the first valid material in the material file to all elements present in the model regular continuum materials are assigned to the regular elements and cohesive materials are assigned to cohesive elements or element facets 4 1 MODEL SETUP 19 To assign different materials to specific elements mesh data information such as tag infor mation or specified pyhsical names can be used MeshDataMaterialSelector class uses this infor
10. Set the options that you need gt make gt make install Or use the Makefile we added for your convenience to handle the cmake configuration gt cd akantu gt make config gt make gt make install All the Akantu options are documented in Appendix C 2 3 Writing a main Function First of all Akantu needs to be initialized The memory management included in the core library handles the correct allocation and de allocation of vectors structures and or objects Moreover in parallel computations the initialization procedure performs the communication setup This is achieved by a pair of functions initialize and finalize that are used as follows include aka_common hh include using namespace akantu int main int argc char argv initialize material dat argc argv your code finalize 4 CHAPTER 2 GETTING STARTED The initialize function takes the material file and the program parameters which can be interpreted by Akantu in due form Obviously it is necessary to include all files needed in main In this manual all provided code implies the usage Of akantu as namespace 2 4 Creating and Loading a Mesh In its current state Akantu supports three types of meshes Gmsh 1 Abaqus 2 and Diana 3 Once a Mesh object is created with a given spatial dimension it can be filled by reading a mesh input file The method read of the class Mesh infers the mesh type from the file extension If
11. Tresca and Gurson can be easily implemented in this class as well In the von Mises yield criterion the hydrostatic stresses have no effect on the plasticity and consequently the yielding occurs when a critical elastic shear energy is achieved 3 tr tr 2 p f Ceff Oy 59 3 90 o E 4 49 f 0 Elastic deformation f 0 Plastic deformation 4 50 where oy is the yield strength of the material which can be function of plastic strain in case of hardening type of materials and o is the deviatoric part of stress given by 1 o 2g str o 1 4 51 32 CHAPTER 4 SOLID MECHANICS MODEL After yielding f 0 thenormality hypothesis of plasticity determines the direction of plastic flow which is normal to the tangent to the yielding surface at the load point Then the tensorial form of the plastic constitutive equation using the von Mises yielding criterion see equation 4 34 may be written as 0 Y 3 2 do 2 Oe ff tr Ae Ap 4 52 In these expressions the direction of the plastic strain increment or equivalently plastic strain rate is given by while the magnitude is defined by the plastic multiplier Ap This can be obtained using the consistency condition which impose the requirement for the load point to remain on the yielding surface in the plastic regime Here we summarize the implementation procedures for the small deformation plasticity with linear isotropic hardening 1 Compute the trial stress
12. and T Hughes Computational Inelasticity Springer 1992 9 J J Marigo Formulation d une loi d endommagement d un materiau lastique C R Acad Sci Paris S r II 292 19 1309 1312 1981 10 J Lemaitre A Course on Damage Mechanics Springer Berlin Heidelberg 1996 11 J Mazars Application de la m canique de l endommagement au comportement non lin aire et la rupture du b ton de structure Ph D thesis Universit Paris 6 1984 12 L Snozzi and J F Molinari A cohesive element model for mixed mode loading with frictional contact capability International Journal for Numerical Methods in Engineering 93 5 510 526 2013 13 G Camacho and M Ortiz Computational modelling of impact damage in brittle materials Inter national Journal of Solids and Structures 33 20 22 2899 2938 1996 14 F Frey and J Jirousek M thodes des elements finis Analyse des structures et milieux continus Trait de g nie civil de l Ecole Polytechnique F d rale de Lausanne PPUR 2009 15 ParaView Open Source Scientific Visualization 16 VisIt Visualization Tool 17 The MayaVi Data Visualizer 18 SCOTCH Static Mapping Graph Mesh and Hypergraph Partitioning 19 T Laursen Computational Contact and Impact Mechanics Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis Engineering online library Springer 2002 75
13. naturally expected from any simulation software even more in the context of parallel computations In order to achieve these goals we made noticeable choices in the architecture of Akantu First we decided to use the object oriented paradigm through C Then in order to prevent extra cost associated to virtual function calls we designed the library as a hybrid architecture with objects at high level layers and vectorization at low level layers Thus Akantu benefits from inheritance and polymorphism mechanisms without the counterpart of having virtual calls within critical loops This coding philosophy which was demonstrated to be highly efficient is innovative in the field of Finite Element software This document is appropriate for researchers and engineers willing to use Akantu in order to perform a finite element calculation for solid mechanics structural mechanics contact mechanics or heat transfer The solid mechanics solver which is the most complete and functional part of Akantu is presented in details in the remainder of this document Chapter 2 Getting Started 2 1 Downloading the Code The Akantu source code can be requested using the form accessible at the URL http 1lsms epf1 ch akantu There you will be asked to accept the LGPL license terms 2 2 Compiling Akantu Akantu is a cmake project so to configure it you can either follow the usual way gt cd akantu gt mkdir build gt cd build gt ccmake
14. some of its properties are listed in Table 3 5 Note Beam elements are of mixed order the axial displacement is linearly interpolated while transverse displacements and rotations use cubic shape functions 2 y 1 2 Figure 3 6 Schematic depiction of a Bernoulli beam element applied to 2D and 3D in Akantu The node numbering as used in Akantu is indicated and also the quadrature points are highlighted gray circles Element type Dimension nodes quad points d o f _bernoulli_beam_2 2D 2 3 6 _bernoulli_beam_3 3D 2 3 12 Table 3 5 Some basic properties of the beam elements in Akantu Chapter 4 Solid Mechanics Model The solid mechanics model is a specific implementation ofthe Model interface dedicated to handle the equations of motion or equations of equilibrium The model is created for a given mesh It will create itsown FEEngine object to compute the interpolation gradient integration and assembly operations A SolidMechanicsModel object can simply be created like this SolidMechanicsModel model mesh where mesh is the mesh for which the equations are to be solved A second parameter called spatial_dimension can be added after mesh if the spatial dimension of the problem is different than that of the mesh This model contains at least the the following six Arrays blocked_dofs contains a Boolean value for each degree of freedom specifying whether that de gree is blocked or not A Dirichlet boundary conditio
15. that case the number of components per tuple must be specified at the Array creation For example if we want to create an Array to store the coordinates sequences of three values of ten nodes the appropriate code is the following UInt nb nodes 10 UInt spatial dimension 3 Array Real position nb nodes spatial dimension 2 5 USING ARRAYS 5 In this case the x position of the eighth node number will be given by position 7 0 inC numbering starts at 0 and not 1 If the number of components for the sequences is not specified the default value of 1 is used It is very common in Akantu to loop over arrays to perform a specific treatment This ranges from geometric calculation on nodal quantities to tensor algebra in constitutive laws for example The Array object has the possibility to request iterators in order to make the writing of loops easier and enhance readability For instance a loop over the nodal coordinates can be performed like accessing the nodal coordinates Array Array lt Real gt nodes mesh getNodes creating the iterators Array lt Real gt vector_iterator it nodes begin spatial_dimension Array lt Real gt vector_iterator end nodes end spatial dimension for it end it Vector lt Real gt amp coords it do what you need In that example each Vector lt Real gt is a geometrical array of size spatial_dimension and the iteration is conveniently performed by
16. the Array iterator The Array object is intensively used to store second order tensor values In that case it should be specified that the returned object type is a matrix when constructing the iterator This is done when calling the begin function For instance assuming that we have a Array storing stresses we can loop over the stored tensors by creating the iterators Array lt Real gt matrix_iterator it stresses begin spatial_dimension spatial_dimension Array lt Real gt matrix_iterator end stresses end spatial_dimension spatial_dimension for it end it Matrix lt Real gt amp stress it do what you need In that last example the Matrix objects are spatial_dimension x spatial_dimension matrices The light objects Matrix and Vector can be used and combined to do most common linear algebra In general a mesh consists of several kinds of elements Consequently the amount of data to be stored can differ for each element type The straightforward example is the connectiv ity array namely the sequences of nodes belonging to each element linear triangular elements have fewer nodes than say rectangular quadratic elements etc A particular data structure called ElementTypeMapArray is provided to easily manage this kind of data It consists of a group of Arrays each associated to an element type The following code can retrieve the ElementTypeMapArray which stores the connectivity arrays for a mesh
17. the elements that are used within that mesh The element types that can be used depend on the mesh but also on the dimensionality of the problem 1D 2D or 3D In Akantu several isoparametric Lagrangian element types are supported and one serendipity element Each of these types is discussed in some detail below starting with the 1D elements all the way to the 3D elements More detailed information shape function location of Gaussian quadrature points and so on can be found in Appendix A 3 1 Isoparametric Elements 1D In Akantu there are two types of isoparametric elements defined in 1D These element types are called _segment_2 and _segment_3 and are depicted schematically in Figure Some of the basic properties of these elements are listed in Table 1 1 2 1 3 2 a _segment_2 b _segment_3 Figure 3 1 Schematic overview of the two 1D element types in Akantu In each element the node numbering as used in Akantu is indicated and also the quadrature points are highlighted gray circles Element type Order nodes quad points _segment_2 linear 2 1 _segment_3 quadratic 3 2 Table 3 1 Some basic properties of the two 1D isoparametric elements in Akantu 2D In Akantu there are four types of isoparametric elements defined in 2D These element types are called _triangle_3 _triangle_6 _quadrangle_4 and _quadrangle_8 and all of them 9 10 CHAPTER 3 ELEMENTS are depicted in Figure 2 As with the 1D elements s
18. the vector x to from v o v mul A x alpha stores the result of Ax in v a is equal to 1 by default o v solve A b stores the result of the resolution of the system Ax b in v o v crossProduct v1 v2 computes the cross product of v1 and v2 and stores the result in v Matrix lt T gt 1 Accessors o AG j gives the component A of the matrix A o A i gives the i column of the matrix as a Vector o A k gives the k component of the matrix matrices are stored in a column major way which means that to access Ajj k i jM o A rows gives the number of rows of A M o A cols gives the number of columns of A N o A size gives the number of component in the matrix M x N 2 Level 1 results are scalars 2 6 MANIPULATING GROUP OF NODES AND OR ELEMENTS 7 1 N o A norm lt N gt returns the Ly norm defined as Xi AG IN N can take any positive integer value o A norm is equivalent to A norm lt L_2 gt o A trace return the trace of A o A det return the determinant of A A doubleDot B return the double dot product of A and B A B o 3 Level 3 results are matrices o A eye s Matrix T eye s fills creates a matrix with the sI with I the identity matrix o A inverse B stores Bina o A transpose returns A o A outerProduct vl v2 stores vv in A o C mul lt t_A t_B gt A B alpha stores the result of the product of A and codeB time the scalar alpha in C t_A a
19. 1 2 2n 20 4 49 40 3 4 4n 4C 3 2 1 0 0 28 1 4 1 0 0 3 0 1 0 n n 1 0 4n 1 0 4 0 0 1 CQC 1 0 0 4 1 5 1 0 0 4E 1 E 9 0 4 88 4g AC 4 4 6 33 0 4En 4n 4 0 7 0 4 0 4n l1 E n O 4n 4 4 8 AC 4n 8 0 0 1 4L 1 n C 47 4C 4 4 4n amp 9 1 0 1 4EC 4C 0 4 10 0 4 4 Anl 0 4C 40 Gaussian quadrature points Coord 9 9 SP SP SB Se oe oP Weight 1 24 1 24 Coord 9 0 22 552 89 959 058 Weight 1 24 1 24 A 3 3D SHAPE FUNCTIONS A 3 3 Hexahedron 8 Element properties Node i Coord 1 0 1 1 1 1 2 1 1 1 3 1 1 1 4 1 1 1 5 1 1 1 6 1 1 1 7 1 1 1 8 1 1 1 Gaussian quadrature points Shape function N Ad 63 Derivative 9N 0 ON dn ON 00 reed SO 9 3 1 90 md 4 0 580 90 0 amp 0 8 mQ 0 g 1 1 1 0 g 1 4 1 m 1 0 1a 50 9 10 2a 09 10 20 9 1a 5a 9 10 230 0 10 20 9 10 90 9 10 20 09 10 2049 1a 90 0 10 230 9 10 20 9 1a 90 0 10 530 0 10 20 93 1a 90 9 10 20 0 10 50 1a 90 9 10 50 0 10 20 9 1a 50 90 10 20 0 10 50 En decken Gene Gee Ce Weight 1 1 1 1 oa ager saan Mr Were ere Weight 1 1 1 1 Appendix B Material parameters B 1 Linear elastic isotropic Keyword elastic Description here 4 4 1 Parameters o rho Real Density o E Real Youn
20. 24 025 026 027 028 029 ooo eo e 02090960 ooo oo Oo o oe oo oo 09090960 6 702e 06 1 832e 07 3 349e 07 5 2e 07 7 211e 07 9 377e 07 1 183e 08 1 41e 08 1 7e 08 963e 08 263e 08 581e 08 907e 08 244e 08 593e 08 051e 08 317e 08 823e 08 237e 08 639e 08 058e 08 483e 08 917e 08 363e 08 681e 08 332e 08 577e 08 281e 08 O OO YN O 00 n UU d 4 4 UU YN NN N ER UJ UJ UJ UJ N ON N SD N ON NON N PRP FA M M FB RP RP RP H 6 78e 09 453e 09 091e 10 171e 10 298e 10 481e 10 624e 10 616e 10 688e 10 678e 10 758e 10 805e 10 821e 10 877e 1 954e 10 131e 10 115e 10 201e 10 315e 10 379e 10 405e 10 485e 10 536e 10 854e 10 099e 10 291e 10 399e 1 426e 10 CHAPTER 9 CONTACT The following figure shows the deformed state of the half sphere at the end of the simulation together with the contact pressure distribution Figure 9 1 in 3D State of pressure and deformation at the end of the simulation of the example of Hertz Appendix A Shape Functions Schmatic overview of all the element types defined in Akantu is described in Chapter B In this appendix more detailed information shape function location of Gaussian quadrature points and so on of each of these types is listed For each element type the coordinates of the nodes are given in the isoparametric frame of reference together with the shape functions and t
21. 7 the resulting tangential one Figure B 5 12 CHAPTER 3 ELEMENTS Figure 3 4 Cohesive element in 2D for quadratic triangular elements T6 gt x S Figure 3 5 Insertion of a cohesive element For the static analysis of the structures containing cohesive elements the stiffness of the cohesive elements should also be added to the total stiffness of the structure Considering a 2D quadratic cohesive element as that in Figure B 4 the opening displacement along the mid surface can be written as No s A s u N s ba os Ni s NKAU PU 3 2 U3 V0 4 9 U5 U2 N s The U A and N are as following U Uo Vo U Vy mn U3 a Uy V4 U U5 3 3 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 O 0 0 10 0 0 0 0 10 0 O fs 0 0 0 1 0 0 0 0 0 10 0 3 4 0000100 0 0 0 1 0 000001 0 0 0 0 O 1 No s 0 Nis 0 Nas 0 k _ No 1 Bom 0 No s 0 NO 0 N Se The consistent stiffness matrix for the element is obtained as oT k PP 3 6 s 96 d5o 3 6 where T is the cohesive traction and the opening displacement for more details check Section 4 3 3 STRUCTURAL ELEMENTS 13 3 3 Structural Elements Bernoulli Beam Elements These elements allow to compute the displacements and rotations of structures constituted by Bernoulli beams Akantu defines them for both 2D and 3D problems respectively in the element types _bernoulli_beam_2 and _bernoulli_beam_3 A schematic depiction of a beam element is shown in Figure B 6Jand
22. ANTU_BLAS This package provides access to a BLAS implementation Under Ubuntu 14 04 LTS the installation can be performed using the following command gt sudo apt get install libatlas base dev AKANTU_IOHELPER This package activates the IOHelper facilities withing Akantu This is mandatory if you want to be able to output Paraview files as well as any Dumper within Akantu AKANTU_SCOTCH This package enables the use the Scotch library in order to perform a graph partitioning leading to the domain decomposition used within Akantu Under Ubuntu 14 04 LTS the installation can be performed using the commands gt sudo apt get install libscotch dev If you activate the advanced option AKANTU_USE_THIRD_PARTY_SCOTCH the make sys tem of akantu can automatically compile Scotch If the automated download fails due to a SSL access not supported by your version of CMake please download the file scotch 5 1 12b esmumps tar gz and then place it in the direc tory lt akantu source gt third party Index Constitutive_laws Elements 9 59 1 8159 Segment 2 59 Segment 3 59 20 960 Quadrangle 4 61 Quadrangle 8 61 Triangle 3 60 Triangle 6 60 30 10 62 Hexahedron 8 63 Tetrahedron 10 62 Tetrahedron 4 62 Bernoulli Isoparametric 9 Structural Forward Euler 46 HeatTransferModel assembleCapacity Lumped 46 explicitCorr explicitPred solveExplicitLumped 47 StableTimeStep updateResidual 47 IO
23. Akantu User s Guide October 6 2014 Version 2 2 Contents 1 Introduction Getting Started 2 1 Downloading the Code 22 Compiling Akahti l cortar 2 3 Writing a main Function 2 4 Creating and Loading a Mesh 2 5 2 5 1 Vector amp Matrix ee ae ea Berne here Manipulating group of nodes and or elements 2 6 1 Ihe NodeGroup object 2 6 2 The ElementGroup object Isoparametric Elements Cohesive Elements Structural Elements 2 6 3 1 3 2 3 3 4 Solid Mechanics Model 4 1 Model Setup 4 1 1 Setting Initial Conditions aso x o ce En serre meer rene A 4 1 2 Setting Boundary Conditions 4 1 3 Material Selector exi 48 51 222242222228 a ee ea eed due 4 1 4 Insertion of Cohesive Elements 4 2 Static Analysis sus de rtr destin y SE done e e Sad here 43 Dynamic Methods 4 3 1 Implicit Time Integration ee s ee ee Neo Hookean B 4 Wisco Elasticity B D l 4 4 4 4 4 Small Deformation Plasticity B 6 T NETTEN ccc HMM n T 4 4 6 Cohesivelaws RR RR A ls 4 5 Adding a New Constitutive Lawl Structural Mechanics Model 5d Model Semple u a uy passa dies ity ee ded enfe lius 5 1 1 Initialization
24. ERIAL_XXX_HH__ define __AKANTU_MATERIAL_XXX_HH__ __BEGIN_AKANTU__ class MaterialXXX public Material declare here the interface of your material Ji In the header file the user also needs to declare all the members of the new material These include the parameters that a read from the material input file as well as any other material parameters that will be computed during the simulation and internal variables In the following the example of adding a new damage material will be presented In this case the parameters in the material will consist of the Young s modulus the Poisson coefficient the resistance to damage and the damage threshold The material will then from these values compute its Lam coefficients and its bulk modulus Furthermore the user has to add a new internal variable damage in order to store the amount of damage at each quadrature point in each step of the simulation For this specific material the member declaration inside the class will look as follows class LocalMaterialDamage public Material declare constructors destructors here declare methods and accessors here Class Members purae e erue Msn A apart e TERE A O O A m S A P TET ere A A E O ery A 4 5 ADDING A NEW CONSTITUTIVE LAW 37 AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST Damage damage Real private the young modulus Real E Poisson coefficient Real nu First Lame coefficient Real lambda Seco
25. K i B D c B 4 82 40 CHAPTER 4 SOLID MECHANICS MODEL Therefore in this method the tangent matrix D is computed for a given strain Note The tangent matrix is a 4 order tensor which is stored as a matrix in Voigt notation Figure 4 16 Tangent to the stress strain curve o getStableTimeStep The stable time step should be calculated based on for the conditionally stable methods This function depends on the longitudinal wave speed which changes for different materials and strains Therefore the value of this velocity should be defined in this method for each new material Note In case of need the calculation of the longitudinal and shear wave speeds can be done in separate functions getPushWaveSpeed and getShearWaveSpeed Therefore the first function can be called for calculation of the stable time step Once the declaration and implementation of the new material has been completed this material can be used in the user s example by including the header file include material_XXX hh For existing materials as mentioned in Section 4 4 by default the materials are initialized inside the method initFull Ifa local material should be used instead the initialization of the material has to be postponed until the local material is registered in the model Therefore the model is initialized with the boolean for skipping the material initialization equal to true model initialization model initFull SolidMechani
26. Plane_Stress bool Plane stress simplification only 2D problems O O O O B 5 Standard linear solid Keyword sls_deviatoric Description here Parameters rho Real Density E Real Young s modulus nu Real Poisson s ratio Plane_Stress bool Plane stress simplification only 2D problems Eta Real Viscosity Ev Real Stiffness of the viscous element OO O O O 0 0 B 6 Elasto plastic linear isotropic hardening Keyword plastic_linear_isotropic_hardening Description here Parameters o rho Real Density o E Real Young s modulus o nu Real Poisson s ratio o h Real Hardening modulus o sigma y Real Yielding stress B 7 Damage Marigo Keyword marigo Description here 4 4 5 Parameters rho Real Density E Real Young s modulus nu Real Poisson s ratio Plane Stress bool Plane stress simplification only 2D problems Yd Random Rupture criterion S Real Damage Energy o o o o o o B 8 DAMAGE MAZARS B 8 Damage Mazars Keyword mazars Description here Parameters o rho Real Density KO Real Damage threshold beta Real Shear parameter o E Real Young s modulus o nu Real Poisson s ratio o At Real Traction post peak asymptotic value o Bt Real Traction decay shape o Ac Real Compression post peak asymptotic value o Bc Real Compression decay shape o o B 9 Cohesive linear Keyword cohesive linear Description here
27. U_IMPLICITI AKANTU_SCOTCH AKANTU_MUMPS AKANTU_COHESIVE_ELEMENT This package activates the cohesive elements engine within Akantu It depends on o A fortran compiler o An implementation of BLAS LAPACK Dependencies AKANTU_LAPACK AKANTU_CONTACT This package enables the contact mechanics engine for Akantu Dependencies AKANTU CPPARRAY AKANTU IMPLICITIAKANTU SCOTCH AKANTU MUMPS AKANTU_OPTIMIZATION AKANTU_NLOPT AKANTU_DAMAGE_NON_LOCAL This package activates the non local damage feature of AKANTU Dependencies AKANTU_LAPACK AKANTU_OPTIMIZATION This activates the optimization routines of Akantu This is cur rently needed by the contact detection algorithms Dependencies AKANTU_NLOPT AKANTU_IMPLICIT This package activates the sparse solver necessary to solve implicitely static dynamic finite element problems It depends on o MUMPS a parallel sparse direct solver o a graph partitioner Dependencies JAKANTU_SCOTCH AKANTU_MUMPS AKANTU_PARALLEL This option activates the parallel features of AKANTU Dependencies JAKANTU_MPI AKANTU_SCOTCH AKANTU_SCOTCH AKANTU_CPPARRAY This package provides access to the open source project If internet is accessible when configuring the project during cmake call this package will be auto downloaded 71 AKANTU_MPI This is a meta package providing access to MPI Under Ubuntu 14 04 LTS the installation can be performed using the commands gt sudo apt get install libopenm
28. a non standard file extension is used the mesh type has to be specified UInt spatial dimension 2 Mesh mesh spatial dimension Reading Gmsh files mesh read my gmsh mesh msh mesh read my gmsh mesh _miot_gmsh Reading Abaqus files mesh read my abaqus mesh inp mesh read my abaqus mesh _miot_abaqus Reading Diana files mesh read my diana mesh dat mesh read my diana mesh miot diana The Gmsh reader adds the geometrical and physical tags as mesh data The physical values are stored as a UInt data called tag 6 if a string name is provided it is stored asa std string data named physical names The geometrical tag is stored as a UInt data named tag 1 The Abaqus reader stores the ELSET in ElementGroups and the NSET in NodeGroups The material assignment can be retrieved from the std string mesh data named abaqus material 2 5 Using Arrays Data in Akantu can be stored in data containers implemented by the Array class In its most basic usage the Array class implemented in Akantu is similar to the vector class of the Standard Template Library STL for C A simple Array containing a sequence of nb element values can be generated with Array type example array nb element where type usually is Real Int UInt or bool Each value is associated to an index so that data can be accessed by typing type amp val example array index Arrays canalso contain tuples of values for each index In
29. a reference only since for some problems it is actually overestimated and convergence cannot be obtained 9 1 1 Implementation In Akantu the object that handles the implicit contact can be found in implicit_contact_manager hh The object that handles the implicit contact resolution stage is the class template 53 54 CHAPTER 9 CONTACT template lt int dim class model_type gt struct ContactData This object takes the command line parameters during construction which can be used to set up the behavior during contact resolution The object can take the following parameters default values in brackets e auto Penalty parameter for augmented Lagrangian formulation alpha 1 Multiplier for values of the penalty parameter utol 0 001 Tolerance used for multipliers in the Uzawa method ntol 0 001 Tolerance used in the Newton Raphson inner convergence loop usteps 100 Maximum number of steps allowed in the Uzawa loop nsteps 100 Maximum number of steps allowed in the Newton Raphson loop Also the following flags can be given to the command line dump Dumping within Newton iterations v Verbose output The ContactData object stores the state of the contact mechanics simulation The state is contained within the following variables sm_ slave master map multipliers_ Lagrange multiplier map areas_ slave areas map penalty_ penalty parameter map gaps_ gap function map model_ reference to solid mechanics model multip
30. al JE O A ELE M et M a ad a PE Pr Ve Ren ses Constructors Destructors 27 NU E RENE CE MM M Rac M MC E RUN aa NO Y public LocalMaterialDamage SolidMechanicsModel amp model const ID id J The user can now define the implementation of the constructor in the material_XXX cc file include local_material_damage hh include solid_mechanics_model hh __BEGIN_AKANTU__ LocalMaterialDamage LocalMaterialDamage SolidMechanicsModel amp model const ID amp id Material model id damage damage this AKANTU_DEBUG_IN this gt registerParam E E 0 pat parsable Young s modulus this registerParam nu nu 0 5 pat parsable Poisson s ratio this registerParam lambda lambda pat readable First Lame coefficient DE this gt registerParam mu mu _pat_readable Second Lame coefficient this gt registerParam kapa kpa _pat_readable Bulk coefficient this gt registerParam Yd Yd 50 pat parsmod this gt registerParam Sd Sd 5000 pat parsmod damage initialize 1 AKANTU_DEBUG_OUT C During the intializer list the reference to the model and the material id are assigned and the constructor of the internal field is called Inside the scope of the constructor the internal values have to be initialized and the parameters that should be printed out are registered with the function registerParam void registerParam name of the parameter key in th
31. arigo This damage evolution law is energy based as defined by Marigo 910 It is an isotropic damage law 1 Y ze C e 4 62 F Y Y Sd 4 63 Y Ya m 3 zi fF gt 0 4 64 unchanged otherwise In this formulation Y is the strain energy release rate Y the rupture criterion and S the damage energy The non local version of this damage evolution law is constructed by averaging the energy Y Mazars B 8 This law introduced by Mazars is a behavioral model to represent damage evolution in concrete The governing variable in this damage law is the equivalent strain eq V lt gt 4 lt gt 4 with lt gt the positive part of the tensor The damage the is defined as D aD 1 a D 4 65 1 A D 1 Rr A exp Pa Ko 4 66 1 A D 1 Sen A exp Pea o 4 67 3 q1 SEF OLE i a Ze 4 68 Eeq With xy the damage threshold A and B the damage parameter in traction A and B the damage parameter in compression f is the shear parameter a is the coupling parameter between traction and compression the e are the eigenstrain and the 44 are the eigenvalues of the strain if the material were undamaged The coefficients A and B are the post peak asymptotic value and the decay shape parameters 4 4 6 Cohesive laws Snozzi Molinari Law B 10 Akantu includes the Snozzi Molinari linear irreversible cohesive law see Figure 4 14 It is an extension to the Camacho Ortiz 13 cohesive law in order
32. aw denoted as a spring displacements as follows 5 Vou ji Vou 4 28 where represents the infinitesimal strain tensor Vou the displacement gradient tensor according to the initial configuration The constitutive equation for isotropic homogeneous media can be expressed as o Atr e I Que 4 29 where o is the Cauchy stress tensor A and u are the the first and second Lame s coefficients In Voigt notation this correspond to 011 1 v v v 0 0 0 11 022 v 1 v v 0 0 0 E22 033 E v v 1 v 0 0 0 33 Se e 4 30 oz 1 v 1 2v 0 0 0 Lo o 2e93 0 013 0 0 0 0 12 O 2813 012 0 0 0 0 0 FF 2en 28 CHAPTER 4 SOLID MECHANICS MODEL Linear anisotropic This formulation is not sufficient to represent all elastic material behavior Some materials have characteristic orientation that have to be taken into account To represent this anisotropy a more general stress strain law has to be used For this we define the elastic modulus tensor as follow 011 C11 C12 C13 C14 C15 C16 11 022 C21 C22 C23 C24 C25 C26 22 033 C31 32 C33 C34 C35 C36 33 4 31 023 C41 C42 C43 C44 C45 C46 2 23 013 C51 C52 C53 C54 C55 C56 2 13 012 C61 C62 C63 Cog C65 Ceo 2812 e2 n ny ey Figure 4 10 Material basis To simplify the writing of input files the C tensor is expressed in the material basis And this basis as to be given too This basis Omar n n2 na is used to define the rotation R nj ej An
33. been discussed how to solve this equation in the static case where 0 Here the method to solve this equation in the general case will be presented For this purpose a time discretization has to be specified The most common discretization method in solid mechanics is the Newmark method which is also the default in Akantu For the Newmark f method becomes a system of three equations see 6 for more details M 1 Citys Ku fos nil 4 10 1 Una Un 1 a At n a Af G a AP iin 4 11 Un 1 Un 1 B At n BAt n 1 4 12 In these new equations n iz and u are the approximations of t t and u t Equa tion is the equation of motion discretized in space finite element discretization and equa tions and are discretized in both space and time Newmark discretization The a and B parameters determine the stability and the accuracy of the algorithm Classical values for a and f are usually 6 1 2 for no numerical damping and 0 lt a lt 1 2 a Method p 1 2 Type O central difference explicit 1 6 Fox Goodwin royal road implicit 1 3 Linear acceleration implicit 1 2 Average acceleration trapezoidal rule implicit The solution of this system of equations 4 10 4 12 is split into a predictor and a corrector system of equations Moreover in the case of a non linear equations an iterative algorithm such as the Newton Raphson method is applied The system of equatio
34. cit Contact Solver The contact formulation corresponds to the augmented Lagrangian method 19 which seeks to minimize the following energy functional between two bodies with contact interface ro 2 1 1 al i 4 i 2 2 II u An II u la E An ENS Na dl 9 1 where II corresponds to the energy functional of the i th body in contact which is a function of the displacement field u and the Lagrangian multiplier Ay In the equation above ey is a penalty parameter g the contact gap function and the Macaulay bracket The solution procedure seeks to make the equation above stationary with respect to both u and An 0 DTP w hen u w n An ENS g dr Vw eV re 9 2 1 0 DaI qn An eng An qu dl Van M EN r A solution is then found by using Uzawa s method for which we solve for u with a fixed GEE r k w f i AQ engu gdr 20 Voeg 9 3 E followed by an update of the multipliers on ro AD AY ewg u 9 4 It is worth noting that the results of the implicit contact resolution depend largely on the choice of the penalty parameter ey Depending on this parameter the computational time needed to obtain a converged solution can be increased dramatically or a convergence solution could not even be obtained at all The code provides a flag that allows the user to rely on an automatic value of en for each slave node Yet this value should be used as
35. csModelOptions _explicit_lumped_mass true Once the model has been initialized the local material needs to be registered in the model model registerNewCustomMaterials lt XXX gt C name_of_local_material Only at this point the material can be initialized model initMaterials A full example for adding anew damage law can be found in examples new_material Chapter 5 Structural Mechanics Model Static structural mechanics problems can be handled using the StructuralMechanicsModel So far Akantu provides 2D and 3D Bernoulli beam elements 14 Just as for the SolidMechanicsModel this model is created for a given Mesh The model will create its own FEEngine object to compute the interpolation gradient integration and assembly operations The StructuralMechanicsModel constructor is called in the following way StructuralMechanicsModel model mesh spatial dimension where mesh is a Mesh object defining the structure for which the equations of statics are to be solved and spatial dimension is the dimensionality of the problem If spatial dimension is omitted the problem is assumed to have the same dimensionality as the one specified by the mesh Note 1 Dynamic computations are not supported to date Note 2 Structural meshes are created and loaded as described in Section 2 4 with MeshIOMSHStruct instead of MeshIOMSH This model contains at least the following Arrays blocked_dofs contains a Boolean value for each degree of fr
36. d C can be rotated in the global basis Q e1 ez es as follow Co RiCo R2 4 32 R R R 2R 2 R sR 3z R 2aR z RuRi3 Ri R12 RaRa R22R22 R23R23 Ro Ros R R2 R21KR Ra R31R31 Rs2R32 Ra3R33 Ra2R33 Ra1R33 R31R32 4 33 RaR3a R22R32 R23R33 R22R33 R2R33 Ra1R32 RuR3a Ri2R32 R13R33 Ri12R33 R11R33 R11R32 RuRa Ri2R22 R13R23 Ri2Ro Ri1R23 RuR2 RuRu RaRa Ra Rz Ra Rz R Rz Ri Rai R aR 2 R22R22 R32R32 R22R32 Ri2R3s2 Ri2R2 Ro Ri3R13 R23R23 R33R33 R23R33 R13R33 R13R23 4 34 R 2R 3z R22R23 R32R33 R22R33 RiR33 R12R33 R R z Ra R23 R31R33 RaiR33 Ru1R33 Ry Raz R R 2 Ra R22 R31R32 Ra1R32 Ri1R32 Ri R22 4 35 Linear orthotropic A particular case of anisotropy is when the material basis is orthogonal in which case the elastic modulus tensor can be simplified and rewritten in terms of 9 independents material parameters 4 4 CONSTITUTIVE LAWS 29 011 jy ty qq 0 0 0 11 022 c2 c3 0 0 0 E22 033 C33 0 0 0 33 4 023 C44 0 0 2623 36 013 sym 55 0 2613 012 C66 2815 C11 E1 1 vosva2 T cag E2 1 viava c33 E3 1 v12v21 1 4 37 C12 E va v31v23 T E2 vi2 v32V13 T 4 38 c13 Fi va1 V21V32 T E2 v13 V21V23 0 4 39 C23 E2 v32 v12V31 1 E3 v23 v21V13 0 4 40 C44 H23 C55 H13 C66 U12 4 41 1 ps 4 42 1 V12V21 V13V31 V32V23 2V21V32V13 The poison ratios follow the rule v vj Ej E 44 2 Neo Hookean The hypere
37. e material file member variable default value optional parameter access permissions description The available access permissions are as follows _pat_internal Parameter can only be output when the material is printed _pat_writable User can write into the parameter The parameter is output when the material is printed _pat_readable User can read the parameter The parameter is output when the material is printed _pat_modifiable Parameter is writable and readable 4 5 ADDING A NEW CONSTITUTIVE LAW 39 o _pat_parsable Parameter can be parsed i e read from the input file o _pat_parsmod Parameter is modifiable and parsable In order to implement the new constitutive law the user needs to specify how the additional material parameters that are not defined in the input material file should be calculated Fur thermore it has to be defined how stresses and the stable time step should be computed for the new local material In the case of implicit simulations in addition the computation of the tangent stiffness needs to be defined Therefore the user needs to redefine the following functions of the parent material void initMaterial for explicit and implicit simulations void computeStress ElementType el_type GhostType ghost_type _not_ghost for implicit simulations void computeTangentStiffness const ElementType amp el_type Array lt Real gt amp tangent_matrix GhostType ghost_type _not_gh
38. eedom specifying whether that de gree is blocked or not A Dirichlet boundary condition can be prescribed by setting the blocked_dofs value of a degree of freedom to true The displacement is computed for all degrees of freedom for which the blocked dofs value is set to false For the remaining degrees of freedom the imposed values zero by default after initialization are kept displacement rotation contains the generalized displacements i e displacements and rotations of all degrees of freedom It can be either a computed displacement for free degrees of freedom or an imposed displacement in case of blocked ones u in the following force moment contains the generalized external forces forces and moments applied to the nodes fex in the following residual contains the difference between the generalized external and internal forces and mo ments On the blocked degrees of freedom residual contains the support reactions r in the following It should be mentioned that at equilibrium residual should be zero on the free degrees of freedom An example to help understand how to use this model will be presented in the next section 41 42 CHAPTER 5 STRUCTURAL MECHANICS MODEL 5 1 Model Setup 5 1 1 Initialization The easiest way to initialize the structural mechanics model is model initFull The method initFull computes the shape functions initializes the internal vectors mentioned above and allocates the memory for the
39. ements the virtual method search struct Assignator public ContactData lt 3 SolidMechanicsModel gt SearchBase typedef Point 3 point type 9 1 IMPLICIT CONTACT SOLVER 55 typedef SolidMechanicsModel model_type typedef ModelElement lt model_type gt master_type model_type model_ Assignator model type amp model model model virtual int search const Real ptr point type p ptr ElementGroup amp rs model getMesh getElementGroup rigid surface for ElementGroup type iterator tit rs firstType tit rs lastType QO tit for ElementGroup const element iterator it rs element begin tit it rs element end tit 1t master type m model triangle 3 it if point has projection to triangle p m point lt 3 gt 0 m point lt 3 gt 1 m point lt 3 gt 2 1 return m element return 1 br The first thing we do in the main file is to add some type definitions int main int argc char argv set dimension static const UInt dim 3 type definitions typedef Point dim point type typedef BoundingBox dim bbox type typedef SolidMechanicsModel model type typedef ModelElement model type master type typedef ContactData dim model type contact type typedef std chrono high resolution clock clock typedef std chrono seconds seconds We initialize the library create a mesh and set the solid mechanics model up initial
40. er the initialization blocked_dofs contains a Boolean value for each degree of freedom specifying whether the de gree is blocked or not A Dirichlet boundary condition can be prescribed by setting the blocked_dofs value of a degree of freedom to true The temperature and the tempera ture_rate are computed for all degrees of freedom where the blocked_dofs value is set to false For the remaining degrees of freedom the imposed values zero by default after initialization are kept residual contains the difference between external and internal heat generations The residual contains the supported heat reactions R Q K T on the nodes that the temperature imposed Only a single material can be specified on the domain A material text file e g material dat provides the material properties as follows heat name material capacity XXX density XXX conductivity XXX XXX where the capacity and density are scalars and the conductivity is specified as a3 x3 tensor 6 2 1 Explicit Dynamic The explicit time integration scheme in Akantu uses a lumped capacity matrix C reducing the computational cost see Chapter 4 This matrix is assembled by distributing the capacity of each element onto its nodes Therefore the resulting C is a diagonal matrix stored in the capacity Array of the model model assembleCapacityLumped Note Currently only the explicit time integration with lumped capacity matrix is impleme
41. ero plane strain otherwise plane stress eta dashpot viscosity and Ev stiffness of the viscous element t 4 46 4 4 CONSTITUTIVE LAWS 31 Note that the current standard linear solid model is applied only on the deviatoric part of the strain tensor The spheric part of the strain tensor affects the stress tensor like an linear elastic material 444 Small Deformation Plasticity The small deformation plasticity is a simple plasticity material formulation which accounts for the additive decomposition of strain into elastic and plastic strain components This formulation is applicable to infinitesimal deformation where the additive decomposition of the strain is a valid approximation In this formulation plastic strain is a shearing process where hydrostatic stress has no contribution to plasticity and consequently plasticity does not lead to volume change Figure 4 13 shows the linear strain hardening elasto plastic behavior according to the additive decomposition of strain into the elastic and plastic parts in infinitesimal deformation as e e e 4 47 o 2G amp Atr e I 4 48 Linear Isotropic Hardening Q a Perfect Plasticity re Figure 4 13 Stress strain curve for the small deformation plasticity with linear isotropic hardening In this class the von Mises yield criterion is used In the von Mises yield criterion the yield is independent of the hydrostatic stress Other yielding criteria such as
42. eters are defined as a base value to which we add a random number that follows the chosen distribution The Uniform distribution is gives a random values between in min max The Weibull distribution is characterized by the following cumulative distribution function F x 1 ed 4 27 4 4 CONSTITUTIVE LAWS 27 which depends on m and A which are the shape parameter and the scale parameter The following sections describe the constitutive models implemented in Akantu In Appendix B a summary of the parameters for all materials of Akantu is provided 4 4 1 Elasticity The elastic law is a commonly used constitutive relationship that can be used for a wide range of engineering materials e g metals concrete rock wood glass rubber etc provided that the strains remain small i e small deformation and stress lower than yield strength The elastic laws are often expressed as o C with where o is the Cauchy stress tensor represents the infinitesimal strain tensor and C is the elastic modulus tensor Linear isotropic The linear isotropic elastic behavior is described by Hooke s law which states that the stress is linearly PES The to the applied strain material behaves like an ideal spring as illustrated in Figure 4 9 The equation that relates the strains to the displacements is point from the VES a b Figure 4 9 a Stress strain curve for elastic material and b schematic representation of Hooke s l
43. f time the method solveStep which has already been used in the static example in Section 4 2 is called inside a time loop time loop Real time 0 for UInt s 1 time lt max_time s time time_step model solveStep lt _scm_newton_raphson_tangent_modified _scc_increment gt le 12 100 24 CHAPTER 4 SOLID MECHANICS MODEL An example of solid mechanics with an implicit time integration scheme is presented in examples implicit implicit_dynamic cc This example consists of a 3D beam of 10m x 1m x 1m blocked on one side and is on a roller on the other side A constant force of 5kN is applied in its middle Figure 4 presents the geometry of this case The material used is a fictitious linear elastic material with a density of 1000 kg m a Young s Modulus of 120 MPa and Poisson s ratio of 0 3 These values were chosen to simplify the analytical solution An approximation of the dynamic response of the middle point of the beam is given by u 5 A 1 COS zt 1 08 a 1 Cos 5579 4 20 Figure 4 6 Numerical setup Figure 4 7 presents the deformed beam at 3 different times during the simulation time steps 0 1000 and 2000 2000 1000 Figure 4 7 Deformed beam at 3 different times displacement are magnified by a factor 10 4 3 2 Explicit Time Integration The explicit dynamic time integration scheme is based on the Newmark f scheme with a 0 see equations 4 10 4
44. g s modulus o nu Real Poisson s ratio o Plane_Stress bool Plane stress simplification only 2D problems B 2 Linear elastic anisotropic Keyword elastic_anisotropic Description here 4 4 1 Parameters o rho Real Density n1 Vector lt Real gt Direction of the main material axis n2 Vector lt Real gt Direction of the second material axis if applicable n3 Vector lt Real gt Direction of the third material axis if applicable C Real Coefficient ij of the material tensor C all the 36 values in Voigt notation can be entered o o o o B 3 Linear elastic orthotropic Keyword elastic orthotropic Description here Parameters rho Real Density n1 Vector lt Real gt Direction of the main material axis n2 Vector lt Real gt Direction of the second material axis if applicable n3 Vector lt Real gt Direction of the third material axis if applicable E1 Real Young s modulus n1 2 o Real Young s modulus n2 Real Young s modulus n3 Real Poisson s ratio 12 nu13 Real Poisson s ratio 13 Real Poisson s ratio 23 G O OO 00 0 0 0 tri N 65 66 APPENDIX B MATERIAL PARAMETERS o G12 Real Shear modulus 12 o G13 Real Shear modulus 13 o G23 Real Shear modulus 23 B 4 Neohookean finite strains Keyword neohookean Description here 4 4 2 Parameters rho Real Density E Real Young s modulus nu Real Poisson s ratio
45. ge parameter D can be defined as D min 2 1 4 76 4 5 ADDING A NEW CONSTITUTIVE LAW 35 which varies from 0 undamaged condition and 1 fully damaged condition This variable can only increase because damage is an irreversible process A simple penalty contact model has been incorporated in the cohesive law so that normal tractions can be returned in case of compression Tn aAn if An 0 4 77 where a is a stiffness parameter that defaults to zero The relative contact energy is equivalent to reversible energy but in compression The material name of the linear decreasing cohesive law is material_cohesive_linear and its parameters with their respective default values are o Sigma_c 0 o beta 0 o G cI 0 o G_cII 0 o kappa 1 o penalty 0 A random number generator can be used to assign a random oc to each facet following a given distribution see Section 4 4 The bilinear constitutive law works exactly the same way as the linear one except for the additional parameter delta_0 that by default is zero Two examples for the extrinsic and in trinsic cohesive elements and also an example to assign different properties to intergranular and transgranular cohesive elements can be found in the folder examples cohesive_element Exponential Cohesive Law Ortiz and Pandolfi proposed this cohesive law in 1999 4 The traction opening equation for this law is as follows T eoe de 4 78 C This equation is plot
46. he dynamic insertion of extrinsic cohesive elements can be utilized in the following way SolidMechanicsModelCohesive model mesh model initFull SolidMechanicsModelCohesiveOptions explicit lumped mass true model limitInsertion x 1 1 20 CHAPTER 4 SOLID MECHANICS MODEL model updateAutomaticInsertion in which this time the method limitInsertion prevents the cohesive elements to be inserted out of the range 1 1 in the x direction In order to check stress and automatically insert elements it is necessary to call the function checkCohesiveStress in the main loop where solveStep is model checkCohesiveStress model solveStep At any time during the simulation it is possible to access the following energies with the relative function Real Ed model getEnergy dissipated Real Er model getEnergy reversible Real Ec model getEnergy contact Statics The only cohesive law that is applicable in this case is the exponential one see section 4 4 6 However unloading reloading cycles are not supported yet In this case cohesive elements have to be inserted before creating the SolidMechanicsModelCohesive model Mesh mesh spatial dimension mesh read implicit mesh msh CohesiveElementInserter inserter mesh inserter setLimit y 0 9 1 1 inserter insertIntrinsicElements SolidMechanicsModelCohesive model mesh model initFull SolidMechanicsModelCohesiveOptions static Also in
47. he boundary T of the problem domain Fig 1 1 t on t Vxel 4 4 4 1 MODEL SETUP 17 Different ways of imposing these boundary conditions exist A basic way is to loop over nodes or elements at the boundary and apply local values A more advanced method consists of using the notion of the boundary of the mesh In the following both ways are presented Starting with the basic approach as mentioned the Dirichlet boundary conditions can be applied by looping over the nodes and assigning the required values Figure 4 2 shows a beam with a fixed support on the left side On the right end of the beam a load is applied At the fixed support the displacement has a given value For this example the displacements in both the x and the y direction are set to zero Implementing this displacement boundary condition is similar to the implementation of initial displacement conditions described above However in order to impose a displacement boundary condition for all time steps the corresponding nodes need to be marked as boundary nodes as shown in the following code Array lt bool gt blocked model getBlockedDOFs const Array lt Real gt pos mesh getNodes UInt nb nodes mesh getNbNodes for UInt i 0 i lt nb nodes i if Math are float equal pos i _x 0 blocked i _x true block displacement in x direction blocked i _y true block displacement in y direction disp i _x 0 fixed displaceme
48. he next sections 15 16 CHAPTER 4 SOLID MECHANICS MODEL L 1 Figure 4 1 Problem domain Q with boundary in three dimensions The Dirchelet and the Neu mann regions of the boundary are denoted with T and T4 respecitvely 4 1 ModelSetup 41 1 Setting Initial Conditions Foraunique solution of the equations of motion initial displacements and velocities for all degrees of freedom must be specified u u 0 0 Q t uo t v 0 4 2 The solid mechanics model can be initialized as follows model initFull This function initializes the internal arrays and sets them to zero Initial displacements and velocities that are not equal to zero can be prescribed by running a loop over the total number of nodes Here the initial displacement in x direction and the initial velocity in y direction for all nodes is set to 0 1 and 1 respectively Array Real amp disp model getDisplacement Array Real amp velo model getVelocity for UInt i 0 i lt mesh getNbNodes i disp i 0 0 1 velo i 1 1 4 1 2 Setting Boundary Conditions This section explains how to impose Dirichlet or Neumann boundary conditions A Dirichlet boundary condition specifies the values that the displacement needs to take for every point x at the boundary T of the problem domain Fig 1 1 u Vxel 4 3 A Neumann boundary condition imposes the value of the gradient of the solution at t
49. heir derivatives on these respective nodes Also all the Gaussian quadrature points within each element are assigned together with the weight that is applied on these points The graphical representations of all the element types can be found in Chapter A 1 1D Shape Functions A 1 1 Segment 2 Element properties Node i Coord Shape function Nj Derivative ON 0 1 i 10 4 i 2 1 1 14 1 Gaussian quadrature points Coord 0 Weight 2 A 1 2 Segment 3 Element properties Node i Coord Shape function Nj Derivative 9N 0 1 1 ee pol 2 1 1E E 1 E 3 0 1 amp 2 Gaussian quadrature points Coord 1 vV3 193 Weight 1 1 59 60 APPENDIX A SHAPE FUNCTIONS A 2 2D Shape Functions A 2 1 Triangle 3 Element properties Node i Coord 17 Shape function N Derivative IN dE ON on 1 0 0 1 n 1 1 2 1 0 1 0 3 0 1 n 0 1 Gaussian quadrature points Coord m 4 1 Weight i A 2 2 Triangle 6 Element properties Node i Coord n Shape function N Derivative IN 0 ON on 1 0 0 ale 1 4 1 n 1 4 1 n 2 1 0 1 2 4 1 0 3 0 1 n 1 2n 0 4n 1 4 so 4E 1 1 4 1 2E m 4 5 1 1 4En 4n 4 6 o 3 An l amp n 74 40 29 Coord n 2 2 3 i i Weight i A 2 2D SHAPE FUNCTIONS 61 A 2 3 Quadrangle 4 Element properties Node i C
50. icity B 5 Visco elasticity is characterized by strain rate dependent behavior Moreover when such a material undergoes a deformation it dissipates energy This dissipation results in a hysteresis loop in the stress strain curve at every loading cycle see FigureH 12a In principle it can be applied to many materials since all materials exhibit a visco elastic behavior if subjected to particular conditions such as high temperatures The standard rheological linear solid model see Sections 10 2 and O a b Figure 4 12 a Characteristic stress strain behavior of a visco elastic material with hysteresis loop and b schematic representation of the standard rheological linear solid visco elastic model 10 3 of 8 has been implemented in Akantu This model results from the combination of a spring mounted in parallel with a spring and a dashpot connected in series as illustrated in Figure 4 12b The advantage of this model is that it allows to account for creep or stress relaxation The equation that relates the stress to the strain is in 1D de t E Ey do t E 7 dott n Le dt n n where n is the viscosity The equilibrium condition is unique and is attained in the limit as t oo At this stage the response is elastic and depends on the Young s modulus E The mandatory parameters for the material file are the following rho density E Young s modulus nu Poisson s ratio Plane_Stress if set to z
51. it inherits its properties from DefaultMaterialSelector This material selector first checks if an element is a cohesive element or a facet If there is then the first cohesive material in the material file is assigned to that element Otherwise this material is assigned to the element facet to which a cohesive element might be inserted later Hence this material selector works for the extrinsic cohesive elements which are dynamically inserted throughout the analysis If the element of interest is not a cohesive element or an element facet then DefaultMaterialSelector is used by default Apart from the Akantu s default material selectors users can always develop their own classes in the main code to tackle various multi material assignment situations 4 1 4 Insertion of Cohesive Elements Dynamics As far as dynamic simulations are concerned cohesive elements are currently compatible only with the explicit time integration scheme see section 4 3 2 They do not have to be inserted when the mesh is generated but during the simulation Intrinsic cohesive elements can be introduced at the beginning of the simulation as follows SolidMechanicsModelCohesive model mesh model initFull model limitInsertion _x 1 1 model insertIntrinsicElements where the insertion is limited to the facets whose barycenter s x coordinate is in the range 1 1 Additional restrictions with respect to y and z directions can be added as well Similarly t
52. ize steel dat argc argv create mesh Mesh mesh dim read mesh mesh read hertz 3D msh create model model type model mesh SolidMechanicsModelOptions opt static initialize material 56 CHAPTER 9 CONTACT model initFull opt Then we create the contact data which can be used to solve the contact problem Note that some of the parameters required by the contact object can be coded in the implementation file these can also be passed as arguments create data structure that holds contact data contact_type cd argc argv model optimal value of penalty multiplier cd Alpha 0 05 cd Uzawa tol 1 e 2 cd Newton tol 1 e 2 Next we find the area that corresponds to each slave node For this we use the fact that if we apply a unit distributed load over the contact surface the resulting force vector at each slave node has magnitude that is equal to the area get physical names from Gmsh file mesh createGroupsFromMeshData std string physical names get areas for the nodes of the circle this is done by applying a unit pressure to the contact surface elements model applyBC BC Neumann FromHigherDim Matrix lt Real gt eye 3 1 contact surface Array lt Real gt amp areas model getForce loop over contact surface nodes and store node areas ElementGroup amp eg mesh getElementGroup contact_surface cout lt lt Adding areas to slave nodes
53. lastic Neo Hookean constitutive law results from an extension of the linear elastic relationship Hooke s Law for large deformation Thus the model predicts nonlinear stress strain behavior for bodies undergoing large deformations O E Figure 4 11 Neo hookean Stress strain curve Asillustrated in Figure 4 11 the behavior is initially linear and the mechanical behavior is very close to the corresponding linear elastic material This constitutive relationship which accounts for compressibility is a modified version of the one proposed by Ronald Rivlin 7 The strain energy stored in the material is given by 1 1 W O 540 In uoIn Zu tr C 3 4 43 where Ao and uo are respectively Lam s first parameter and the shear modulus at the initial configuration J is the jacobian of the deformation gradient F We J det F Finally C is the right Cauchy Green deformation tensor 30 CHAPTER 4 SOLID MECHANICS MODEL Since this kind of material is used for large deformation problems a finite deformation frame work should be used Therefore the Cauchy stress 0 should be computed through the second Piola Kirchhoff stress tensor S o ESF 4 44 Finally the second Piola Kirchhoff stress tensor is given by OV 1 1 S 2236 AolnJC uo I C 4 45 The parameters to indicate in the material file are the same as those for the elastic case E Young s modulus nu Poisson s ratio 44 3 Visco Elast
54. lement type Order nodes quad points _tetrahedron_4 linear 4 1 _tetrahedron_10 quadratic 10 hexahedron 8 cubic 8 8 Table 3 3 Some basic properties of the three 3D isoparametric elements in Akantu 3 2 Cohesive Elements The cohesive elements that have been implemented in Akantu are based on the work of Ortiz and Pandolfi 4 Their main properties are reported in Table 3 4 Element type Facet type Order nodes quad points _cohesive_ld_2 point 1 linear 2 1 _cohesive_2d_4 _segment_2 linear 4 1 _cohesive_2d_6 segment 3 quadratic 6 2 cohesive 3d 6 triangle 3 linear 6 1 cohesive 3d 12 triangle 6 quadratic 12 3 Table 3 4 Some basic properties of the cohesive elements in Akantu Cohesive element insertion can be either realized at the beginning of the simulation or it can be carried out dynamically during the simulation The first approach is called intrinsic the second one extrinsic When an element is present from the beginning a bilinear or exponential cohesive law should be used instead of a linear one A bilinear law works exactly like a linear one except for an additional parameter 69 separating an initial linear elastic part from the linear irreversible one For additional details concerning cohesive laws see Section 4 4 6 Extrinsic cohesive elements are dynamically inserted between two standard elements when Oe gt 0 With oe 4102 3 1 in which oy is the tensile normal traction and
55. let exist FlagOnly FixedValue and IncrementValue The functor FlagOnly is used if a point is fixed in a given direction Therefore the input parameter to this functor is only the fixed direction The FixedValue functor is used when a displacement value is applied in a fixed direction The IncrementValue applies an increment to the displacement in a given direction The following 18 CHAPTER 4 SOLID MECHANICS MODEL code shows the utilization of three functors for the top bottom and side surface of the mesh which were already defined in the Gmsh file mesh createGroupsFromMeshData std string physical names model applyBC BC Dirichlet FixedValue 13 0 y Top model applyBC BC Dirichlet FlagOnly x Bottom model applyBC BC Dirichlet IncrementValue 13 0 x Side To apply a Neumann boundary condition the applied traction or stress should be specified before In case of specifying the traction on the surface the functor FromTraction of Neumann boundary conditions is called Otherwise the functor FromStress should be called which gets the stress tensor as an input parameter Array lt Real gt surface traction 3 surface traction 0 surface traction 1 surface traction 2 1 Array lt Real gt surface stress 3 3 0 surface stress 80 0 0 surface stress 1 1 0 surface stress 2 2 1 model applyBC BC Neumann FromTraction surface traction Bottom model applyBC BC
56. lier_dumper_ structures used to dump multipliers pressure_dumper_ structures used to dump pressure options_ flags_ uiter_ niter_ options map flags map Uzawa and Newton iteration counters The interface of the ContactData object contains three methods to solve for each contact step which is overloaded depending on the parameters passed Their signatures are as follows void solveContactStep void solveContactStep search type search template class PostAssemblyFunctor gt void solveContactStep search type search const PostAssemblyFunctor amp fn The second method allows the user to provide a pointer to an object that is used to search slave master pairs This can be done for example when due to the deformed configuration current slave master pairs are no longer valid The last method in the snippet above allows the user to provide a functor that is called after the assembly of the contact contributions to the stiffness matrix and the force vector The last takes place within the method computeTangentAndResidual Hertz Example Here we outline step by step the use of the implicit contact solver to obtain the solution of Hertzian contact The complete implementation can be found in examples contact hertz 3D cc The following class is used as the object that will perform the search for new contact elements when a slave node is found to lie outside its master element The class derives from a SearchBase class and impl
57. lt lt endl for auto nit eg node_begin nit eg node_end nit compute area contributing to the slave node Real a 0 for UInt i 0 i lt dim i a pow areas nit i 2 cd addArea nit sqrt a set force value to zero areas clear Note that we clear the force vector after the assignment of areas In the next step we prescribe the boundary conditions that do not change in time apply boundary conditions for the rigid plane model applyBC BC Dirichlet FixedValue 6 BC _x bottom body model applyBC BC Dirichlet FixedValue 0 BC _y bottom body model applyBC BC Dirichlet FixedValue 0 BC z bottom body block z disp in extreme points of top surface model getBoundary 1 2 true model getBoundary 2 2 true block x disp in extreme points of top surface model getBoundary Q 3 0 true model getBoundary 4 0 true 9 1 IMPLICIT CONTACT SOLVER 57 Next we add the slave master pairs for the analysis We use a bounding box to consider only a fraction of the slave nodes in the model Those slave nodes that are not within the bounding box are not considered in the analysis set up bounding box to include slave nodes that lie inside it Real 11 1 Real 12 0 2 Real 13 1 point_type c1 11 2 12 2 13 2 point type c2 11 2 12 2 13 2 bbox_type bb c1 c2 search policy for slave master pairs Assignator a model
58. mation to assign different materials With the proper physical name or tagname and index different materials can be assigned as demonstrated in the examples below MeshDataMaterialSelector lt std string gt mat selector mat selector new MeshDataMaterialSelector std string physical names model model setMaterialSelector mat selector In this example the physical names specified in a GMSH geometry file will by used to match the material names in the input file Another example would be MeshDataMaterialSelector UInt mat selector mat selector new MeshDataMaterialSelector UInt tag 1 model model setMaterialSelector mat selector where tag 1 of the mesh file is used as the classifier index and the elements that have index value equal to one will be assigned to the second material in the material file There are four different material selectors pre defined in Akantu MaterialSelector and DefaultMaterialSelector is used to assign a material to regular elements by default For the regular elements as in the example above MeshDataMaterialSelector can be used to assign different materials to different elements However for cohesive elements it is not possible to assign multiple cohesive materials automatically and user has to write a custom class Akantu has a pre defined material selector to assign the first cohesive material by default to the cohesive elements which is called DefaultMaterialCohesiveSelector and
59. n 2D the memory has Vector s of 2 components or the 2 order tensors with 2 x 2 components This memory can be dealt with addDumpFieldVector which always dumps Vectors with 3 components or addDumpFieldTensor which dumps 2 order tensors with 3 x 3 components respectively The routines addDumpFieldVector and addDumpFieldTensor were introduced because of Paraview which mostly manipulate 3D data Those fields which are stored by quadrature point are modified to be seen in the VTK file as elemental data To do this the default is to average the values of all the quadrature points The list of fields depends on the models for SolidMechanicsModel see table 7 1 The user can also register external fields which have the same mesh as the mesh from the model as support To do this an object oftype Field has to be created o For nodal fields Vector lt T gt vect nb_nodes nb_component Field field new DumperIOHelper NodalField lt T gt vect model addDumpFieldExternal my field field o For elemental fields ElementTypeMapArray lt T gt arr Field field new DumperIOHelper ElementalField lt T gt arr spatial displacement 49 50 CHAPTER 7 INPUT OUTPUT key type support displacement Vector lt Real gt nodes velocity Vector lt Real gt nodes acceleration Vector lt Real gt nodes force Vector lt Real gt nodes residual Vector lt Real gt nodes boundary Vector lt bool gt nodes mass Vector lt Real gt nodes
60. n can be prescribed by setting the blocked_dofs value of a degree of freedom to true A Neumann boundary condition can be applied by setting the blocked_dofs value of a degree of freedom to false The dis placement velocity and acceleration are computed for all degrees of freedom for which the blocked_dofs value is set to false For the remaining degrees of freedom the imposed values zero by default after initialization are kept displacement contains the displacements of all degrees of freedom It can be either a computed displacement for free degrees of freedom or an imposed displacement in case of blocked ones u in the following velocity contains the velocities of all degrees of freedom As displacement it contains computed or imposed velocities depending on the nature of the degrees of freedom in the following acceleration contains the accelerations of all degrees of freedom As displacement it contains computed or imposed accelerations depending on the nature of the degrees of freedom ii in the following force contains the external forces applied on the nodes f in the following ext residual contains the difference between external and internal forces On blocked degrees of freedom residual contains the support reactions rin the following Itshould be mentioned that at equilibrium residual should be zero on free degrees of freedom Some examples to help to understand how to use this model will be presented in t
61. n select nodes having a positive X coordinate with the following code 8 CHAPTER 2 GETTING STARTED Array lt Real gt amp nodes mesh getNodes NodeGroup amp group mesh createNodeGroup XpositiveNode Array lt Real gt const_vector_iterator it nodes begin spatial_dimension Array lt Real gt const_vector_iterator end nodes end spatial dimension UInt index 6 for it end Hit index const Vector lt Real gt amp position it if position 0 gt 0 group add index 2 6 2 The ElementGroup object A group of elements is stored in ElementGroup objects Since a group can contain elements of various types the ElementGroup object stores indexes in a ElementTypeMapArray lt UInt gt object Then elements can be added to the group by calling addElement For instance selecting the elements for which the barycenter of the nodes has a positive X coordinate can be made with ElementGroup amp group mesh createElementGroup XpositiveElement Mesh type_iterator it mesh firstType Mesh type_iterator end mesh lastTypeQ Vector lt Real gt barycenter spatial dimension for it end it UInt nb element mesh getNbElement it for UInt e 0 e lt nb element e ElementType type it mesh getBarycenter e type barycenter storage if barycenter 80 gt 0 group add type e Chapter 3 Elements The base for every Finite Elements computation is its mesh and
62. nd Lame coefficient shear modulus Real mu resistance to damage Real Yd damage threshold Real Sd Bulk modulus Real kpa damage internal variable InternalField lt Real gt damage m In order to enable to print the material parameters at any point in the user s example file using the standard output stream by typing for Ulnt m Q m lt model getNbMaterials m std cout model getMaterial m std endl the standard output stream operator has to be redefined This should be done at the end of the header file class LocalMaterialDamage public Material declare here the interace of your material FE VAS uA Rr MM Rp E SS V DER LS pn ACID EUR CR op e SRI EE VE ee AP een mu SO e MEE br Re REM t A inline functions DA PE E WEE UM A A A NA AI EE RIA AA ERAS AS A EA AS EA er 2 standard output stream operator inline std ostream amp operator lt lt std ostream amp stream const LocalMaterialDamage amp _this _this printself stream return stream However the user still needs to register the material parameters that should be printed out The registration is done during the call of the constructor Like all definitions the implementation of the constructor has to be written in the material_XXX cc file However the declaration has to be provided in the material_XXX hh file 38 CHAPTER 4 SOLID MECHANICS MODEL class LocalMaterialDamage public Materi
63. nd t_B are boolean defining if A and B should be transposed or not t_A t_B result false false C aAB false true C aAB true false C aA B true true C aA B o A eigs d V this method computes the eigenvalues and eigenvectors of A and store the results in d and V such that dG A and V i vj with Av Ajvj and A1 gt gt A gt gt AN 2 6 Manipulating group of nodes and or elements Akantu provides the possibility to manipulate subgroups of elements and nodes Any ElementGroup and or NodeGroup must be managed by a GroupManager Such a manager has the role to associate group objects to names This is a useful feature in particular for the application of the boundary conditions as will be demonstrated in section 4 1 2 To most general group manager is the Mesh class which inheritates from the GroupManager class For instance the following code shows how to request an element group to a mesh request creation of a group of nodes NodeGroup amp my node group mesh createNodeGroup my node group request creation of a group of elements ElementGroup amp my element group mesh createElementGroup my element group DE fill and use the groups 2 6 1 The NodeGroup object A group of nodes is stored in NodeGroup objects They are quite simple objects which store the indexes of the selected nodes in a Array lt UInt gt Nodes are selected by adding them when calling NodeGroup add For instance you ca
64. ns can be written as 4 3 DYNAMIC METHODS 23 1 Predictor 0 NT Wiad Un At D Mn 4 13 HD n Ati 4 14 sus 4 15 2 Solve cM dC eK w foren 7 fing 7 Chua 7 Mi fhri 4 16 3 Corrector il cw 4 17 il dw 4 18 Hl u ew 4 19 where i is the Newton Raphson iteration counter and c d and e are parameters depending on the method used to solve the equations w e d C in acceleration Oi afa BAt 1 in velocity bu At 1 aA in displacement u 1 At AP 4 3 1 Implicit Time Integration To solve a problem with an implicit time integration scheme firsta SolidMechanicsModel object has to be created and initialized Then the initial and boundary conditions have to be set Every thing is similar to the example in the static case Section 42 however in this case the implicit dynamic scheme is selected at the initialization of the model SolidMechanicsModel model mesh model initFull SolidMechanicsModelOptions implicit dynamic Boundary conditions see Section Because a dynamic simulation is conducted an integration time step At has to be specified In the case of implicit simulations Akantu implements a trapezoidal rule by default That is to say a 1 2 and f 1 2 which is unconditionally stable Therefore the value of the time step can be chosen arbitrarily within reason model setTimeStep time_step Since the system has to be solved for a given amount o
65. nt in x direction disp i _y 0 fixed displacement in y direction LLL Figure 4 2 Beam with fixed support For the more advanced approach one needs the notion of a boundary in the mesh Therefore the boundary should be created before boundary condition functors can be applied Generally the boundary can be specified from the mesh file or the geometry For the first case the function createGroupsFromMeshData is called This function can read any types of mesh data which are provided in the mesh file If the mesh file is created with Gmsh the function takes one input strings which is either tag_0 tag_1 or physical_names The first two tags are assigned by Gmsh to each element which shows the physical group that they belong to In Gmsh it is also possible to consider strings for different groups of elements These elements can be separated by giving a string physical_names to the function createGroupsFromMeshData Boundary conditions can also be created from the geometry by calling createBoundaryGroupFromGeometry This function gathers all the elements on the boundary of the geometry To apply the required boundary conditions the function applyBC needs to be called on a SolidMechanicsModel This function gets a Dirichlet or Neumann functor and a string which specifies the desired boundary on which the boundary conditions is to be applied The functors specify the type of conditions to apply Three built in functors for Dirich
66. nted within Akantu The explicit integration scheme is Forward Euler 5 o Predictor T 1 Ty AIT m t o Update residual R44 Qs x KT 1 o Corrector Thi C R44 The explicit integration scheme is conditionally stable The time step has to be smaller than the stable time step and it can be obtained in Akantu as follows time step model getStableTimeStep The stable time step is defined as peo Aterit 2A A 6 5 where Ax is the characteristic length e g the inradius in the case of linear triangle element p is the density x is the conductivity tensor and c is the specific heat capacity It is necessary to impose a time step which is smaller than the stable time step for instance by multiplying the stable time step by a safety factor smaller than one 6 2 USING THE HEAT TRANSFER MODEL 47 const Real safety_time_factor 0 1 Real applied_time_step time_step safety_time_factor model setTimeStep applied_time_step The following loop allows for each time step to update the temperature residual and temperature_rate fields following the previously described integration scheme for Ulnt s 1 s 1 applied_time_step lt total_time s model explicitPred model updateResidual model explicitCorr An example of explicit dynamic heat propagation is presented in examples heat_transfer explicit_heat_transfer cc This example consists of a square 2D plate of 1m ha
67. of them In what follows we show how it works on a SolidMechanics model 8 1 Initializing the Parallel Context The user must initialize Akantu by forwarding the arguments passed to the program by using the function initialize and close Akantu instances at the end of the program by calling the finalize function Note This step does not change from the sequential case as it was stated in Section 2 3 It only gives a additional motivation in the parallel MPI context The initialize function builds a StaticCommunicator object responsible for handling the interprocess communications later on The StaticCommunicator can for instance be used to ask the total number of declared processors available for computations as well as the process rank through the functions getNbProc and whoAmI respectively An example of the initializing sequence and basic usage of the StaticCommunicator is int main int argc char argv initialize material dat argc argv StaticCommunicator amp comm StaticCommunicator getStaticCommunicator Int psize comm getNbProc Int prank comm whoAmI finalize 51 52 CHAPTER 8 PARALLEL COMPUTATION 8 2 Partitioning the Mesh The mesh is partitioned after the correct initialization of the processes playing a role in the computation We assume that a Mesh object is constructed as presented in Section 2 4 Then a partition must be computed by using an appropriate mesh partitioner At present
68. ome of the most basic properties of these elements are listed in Table 2 It is important to note that the first element is linear the next two quadratic and the last one cubic Furthermore the last element type _quadrangle_8 is not a Lagrangian but a serendipity element n 7 3 3 6 5 amp 1 2 1 4 2 a triangle 3 b triangle 6 n 7 4 3 1 2 c _quadrangle_4 d quadrangle 8 Figure 3 2 Schematic overview of the four 2D element types in Akantu In each element the node numbering as used in Akantu is indicated and also the quadrature points are highlighted gray circles Element type Order nodes quad points _triangle_3 linear 3 1 _triangle_6 quadratic 6 3 _quadrangle_4 quadratic 4 4 _quadrangle_8 cubic 8 9 Table 3 2 Some basic properties of the four 2D isoparametric elements in Akantu 3D In Akantu there are three types of isoparametric elements defined in 3D These element types are called _tetrahedron_4 _tetrahedron_10 and _hexahedron_8 and all of them are depicted schematically in FigureB 3 As with the 1D and 2D elements some of the most basic properties of these elements are listed in TableB 3 3 2 COHESIVE ELEMENTS 11 a _tetrahedron_4 b _tetrahedron_10 c _hexahedron_8 Figure 3 3 Schematic overview of the three 3D element types in Akantu In each element the node numbering as used in Akantu is indicated and also the quadrature points are highlighted gray spheres E
69. oord y Shape function N Derivative IN od amp ON dn 1 1 1 i 1 amp 1 n 11 n 1a 9 2 1 1 1 amp 1 n 4 1 n 10 8 a ie 0 40 9 40 8 x VERI 1 8 17 ia 9 10 9 Gaussian quadrature points En ama Ge Ge ea aa Weight 1 1 1 1 A 2 4 Quadrangle 8 Element properties Node i Coord n Shape function N Derivative AN 0 IN 0n 1 1 1 10 80 9C1 i 9 4 mM E n 1a 5 29 2 1 1 ia0 20 9C1 5 i 0 90 5 i0 26 29 3 1 1 Id 9 A mMC1l E n i 9Q m 1a42 720 4 31 1 i 830 9Cc1 5 0 mQ 0 10 5 6 29 5 0 1 i 2 1 0 9 1a0 2 6 1 0 1 9 1 7 10 5 n 1 E 7 0 1 3 1 amp 1 n en 30 8 8 1 0 1 1 17 404 P 7 1 E Gaussian quadrature points Weight 64 81 25 81 25 81 25 81 25 81 coord 6 9 048 Jbo e 48 NE Weight 40 81 40 81 40 81 40 81 62 APPENDIX A SHAPE FUNCTIONS A 3 3D Shape Functions A 3 1 Tetrahedron 4 Element properties Node i Coord 17 C Shape function N Derivative ON dE ON 0n ON 0T 1 0 0 0 1 1 1 1 1 2 1 0 0 1 0 0 3 0 1 0 n 0 1 0 4 0 0 1 C 0 0 1 Gaussian quadrature points Coord n C 4 1 Weight A 3 2 Tetrahedron 10 Element properties Node i Coord n Shape function Nj Derivative N ONi On ON 00 4 amp 4n 46 8 1 0 0 0 1 9 C
70. ost for explicit and implicit simulations Real getStableTimeStep Real h const Element amp element In the following a detailed description of these functions is provided o initMaterial This method is called after the material file is fully read and the elements corresponding to each material are assigned Some of the frequently used constant parame ters are calculated in this method For example the Lam constants of elastic materials can be considered as such parameters o computeStress In this method the stresses are computed based on the constitutive law as a function of the strains of the quadrature points For example the stresses for the elastic material are calculated based on the following formula o Atr e I 2ue 4 81 Therefore this method contains a loop on all quadrature points assigned to the material using the two macros MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END MATERIAL STRESS QUADRATURE POINT LOOP BEGIN element type sigma f grad u MATERIAL STRESS QUADRATURE POINT LOOP END Note The strain vector in Akantu contains the values of Vu Le it is really the displacement gradient o computeTangentStiffness This method is called when the tangent to the stress strain curve is desired see Fig 4 16 For example it is called in the implicit solver when the stiffness matrix for the regular elements is assembled based on the following formula
71. partitions Real elements stress Matrix lt Real gt quadrature points strain Matrix lt Real gt quadrature points material internals variable quadrature points Table 7 1 List of dumpable fields for SolidMechanicsModel model addDumpFieldExternal my field field 7 2 Cohesive elements data Cohesive elements and their relative data can be easily dumped thanks to a specific dumper contained in SolidMechanicsModelCohesive In order to use it one has just to add the string cohesive elements when calling each method already illustrated Here is an example on how to dump displacement and damage model setBaseNameToDumper cohesive elements cohesive elements output model addDumpFieldVectorToDumper cohesive elements displacement model addDumpFieldToDumper cohesive elements damage model dump cohesive elements Chapter 8 Parallel Computation This section explains how to launch a parallel computation The strategy adopted by Akantu uses a mesh partitioning where elements are mapped to processors Mesh partitions are then distributed to available processors by adequate routines as will be described below The sequence of additional operations to be performed by the user are o Initializing the parallel context o Partitioning the mesh o Distributing mesh partitions After these steps the Model object proceeds with the interprocess communication automat ically without the user having to explicitly take care
72. pi dev Under Mac OS X the installation requires the following steps gt sudo port install mpich devel Dependencies AKANTU_SCOTCH AKANTU_NLOPT This package enable the use of the optimization alogorithm library Since there are no packaging for the common operating system by default Akantu compiles it as a third party project This behavior can be modified with the option AKANTU_USE_THIRD_PARTY_NLOPT If the automated download fails please download nlopt 2 4 2 tar gz and place it in akantu source third party download AKANTU MUMPS This package enables the MUMPS parallel direct solver for sparce matrices This is necessary to solve static or implicit problems Under Ubuntu 14 04 LTS the installation can be performed using the commands sudo apt get install libmumps seq dev for sequential sudo apt get install libmumps dev for parallel Under Mac OS X the installation requires the following steps gt sudo port install mumps If you activate the advanced option AKANTU_USE_THIRD_PARTY_MUMPS the make system of akantu can automatically compile MUMPS For this you will have to download MUMPS from http mumps enseeiht fr or http graal ens lyon fr MUMPS and place it in lt akantu source gt third party AKANTU_LAPACK This package provides access to a LAPACK implementation Under Ubuntu 14 04 LTS the installation can be performed using the following command gt sudo apt get install libatlas base dev AK
73. ry The strong form of the dynamic heat equation can be expressed as peot V K VT b 6 1 with T the scalar temperature field c the specific heat capacity p the mass density x the conduc tivity tensor and b the heat generation per unit of volume The discretized weak form with a finite number of elements is Vi CoN N dQ T KVN VN dO T N qa nd bN dO 6 2 PCoN j j j q Q j Q T Q j with i and j the node indices n the normal field to the surface T oO To simplify we can define the capacity and the conductivity matrices as c pony ndo and Kij vx vao 6 3 Q Q and the system to solve can be written T O _k T 6 4 with Q the consistent heat generated 6 2 Using the Heat Transfer Model Currently the HeatTransferModel object uses dynamic analysis with an explicit time integration scheme It can simply be created like this HeatTransferModel model mesh spatial dimension while an existing mesh has been used see 2 4 Then the model object can be initialized with model initFull where a material file name has to be provided This function will load the material properties and allocate initialize the nodes and element Arrays More precisely the heat transfer model contains 4 Arrays 45 46 CHAPTER 6 HEAT TRANSFER MODEL temperature contains the nodal temperature T zero by default after the initialization temperature_rate contains the variations of temperature T zero by default aft
74. s performed using the same Newton Raphson algorithm as described in Section 4 2 Note To date StructuralMechanicsModel handles only constitutively and geometrically linear problems the algorithm is therefore guaranteed to converge in two iterations model updateResidual model solve o model updateResidual assembles the internal forces and removes them from the external forces o model solve solves the Equation 5 1 The increment vector of the model will contain the new increment of displacements and the displacement_rotation vector is also updated to the new displacements At the end of the analysis the final solution is stored in the displacement_rotation vec tor A full example of how to solve a structural mechanics problem is presented in the code examples structural_mechanics bernoulli_beam_2_example cc This example is composed of a 2D beam clamped at the left end and supported by two rollers as shown in Figure 5 1 The problem is defined by the applied load q 6 KN m moment M 3 6 kN m moments of inertia I 250000 cm and I 128 000 cm and lengths L 10m and L 8 m The resulting rotations at node two and three are p2 0 001 167 and p3 0 000 771 44 CHAPTER 5 STRUCTURAL MECHANICS MODEL pk s A Figure 5 1 2D beam example Chapter 6 Heat Transfer Model The heat transfer model is a specific implementation of the Model interface dedicated to handle the dynamic heat equation 6 1 Theo
75. splacement is zero scc increment which also means that the equilibrium is found For a linear elastic problem the solution is obtained in one iteration and therefore the maximum number of iterations can be set to one But for a non linear case one needs to iterate as long as the norm of the residual exceeds the tolerance threshold and therefore the maximum number of iterations has to be higher e g 100 model solveStep scm newton raphson tangent modified scc residual 1e 4 100 At the end of the analysis the final solution is stored in the displacement vector A full example of how to solve a static problem is presented in the code examples static static cc This example is composed of a 2D plate of steel blocked with rollers on the left and bottom sides as shown in Figure The nodes from the right side of the sample are displaced by 0 01 of the length of the plate Figure 4 4 Numerical setup The results of this analysis is depicted in Figure 4 5 22 CHAPTER 4 SOLID MECHANICS MODEL Figure 4 5 Solution of the static analysis Left the initial condition right the solution deforma tion magnified 50 times 4 3 Dynamic Methods Different ways to solve the equations of motion are implemented in the solid mechanics model The complete equations that should be solved are M C Ku f 4 9 ext where M C and K are the mass damping and stiffness matrices respectively In the previous section it has already
76. ssigned to the problem one has to specify the material characteristics constitutive behavior and material properties in a text file e g material dat as follows material constitutive_law lt optional flavor gt name value rho value where constitutive_law is the adopted constitutive law followed by the material properties listed one by line in the bracket e g name and density rho Some constitutive laws can also have an optional flavor For example a non local constitutive law can be flavored by a weight function The file needs to be loaded in Akantu using the initialize method of Akantu as shown in Section initialize material dat argc argv In order to conveniently store values at each quadrature in a material point Akantu pro vides a special data structure the InternalField The internal fields are inheriting from the ElementTypeMapArray Furthermore itprovides several functions for initialization auto resizing and auto removal of quadrature points Sometimes it is also desired to generate random distributions of internal parameters An example might be the critical stress at which the material fails To generate such a field in the material input file a random quantity needs be added to the base value sigma_c base sigma_c base uniform min max sigma c base weibull A m All parameters are real numbers For the uniform distribution minimum and maximum values have to be specified Random param
77. stiffness matrix Material properties are defined using the StructuralMaterial structure described in Ta ble 5 1 Such a definition could for instance look like StructuralMaterial matl mat E 3e10 mat I 0 0025 mat A 0 01 Field Description E Young s modulus A Cross section area I Second cross sectional moment of inertia for 2D elements Iy I around beam y axis for 3D elements Iz I around beam z axis for 3D elements GJ Polar moment of inertia of beam cross section for 3D elements Table 5 1 Material properties for structural elements as defined by the structure StructuralMaterial Materials can be added to the model s element material vector using model addMaterial mat1 They are successively numbered and then assigned to specific elements for UInt i 0 i lt nb element mat 1 i model getElementMaterial _bernoulli_beam_2 1 0 1 5 1 2 Setting Boundary Conditions As explained before the Dirichlet boundary conditions are applied through the array blocked_dofs Two options exist to define Neumann conditions If anodal force is applied it has to be directly set in the array force_momentum For loads distributed along the beam length the method computeForcesFromFunction integrates them into nodal forces The method takes as in put a function describing the distribution of loads along the beam and a BoundaryFunctionType specifing if the function is expressed in the local coordinates _bft_trac
78. t they are initialized during the initFull The method initFull also initializes all appropriate vectors to zero Once the model is created and initialized the boundary conditions can be set as explained in Section 4 1 2 Boundary conditions will prescribe the external forces for some free degrees of freedom f and displacements for some others At this point of the analysis the function solveStep can be called model solveStep lt _scm_newton_raphson_tangent_modified _scc_residual gt le 4 1 This function is templated by the solving method and the convergence criterion and takes two arguments the tolerance and the maximum number of iterations which are 1 x 10 and 1 for this example The modified Newton Raphson method is chosen to solve the system In this method the equilibrium equation is modified in order to apply a Newton Raphson convergence algorithm Kitla tl r 4 6 foa z Jint ng feat Ku ur ul Suit where u is the increment of displacement to be added from one iteration to the other and i is the Newton Raphson iteration counter By invoking the solveStep method in the first step the global stiffness matrix K from Equation is automatically assembled A Newton Raphson iteration is subsequently started K is updated according to the displacement computed at the previous iteration and one loops until the forces are balanced scc residual ie r lt scc residual One can also iterate until the increment of di
79. tallation requires the following steps o Install Xcode o Install the command line tools o Install the MacPorts project which allows to automatically download and install opensource packages Then the following commands should be typed in a terminal gt sudo port install cmake gcc48 boost AKANTU_CORE_CXX11 This option activates some features of the C 11 standard This is usable with GCC gt 4 7 or Intel gt 13 AKANTU_DOCUMENTATION_DOXYGEN This generates the Doxygen documantation of the source code It depends on o Doxygen an automated source code documentations system o Optional Graphviz to generate the dependencies graph P p 8 P grap Under Ubuntu 14 04 LTS the installation of the dependencies can be performed using the following command gt sudo apt get install doxygen gt sudo apt get install graphviz 69 70 APPENDIX C PACKAGE DEPENDENCIES AKANTU_DOCUMENTATION_MANUAL This package alows to compile the user manual in the build folder build doc manual manual pdf Under Ubuntu 14 04 LTS the installation of the dependencies can be performed using the following command gt sudo apt get install install rubber texlive texlive science texlive latex extra AKANTU_HEAT_TRANSFER This package activates the heat transfer model within Akantu It has no additional dependencies AKANTU_STRUCTURAL_MECHANICS This package activates the compilation for the Struc tural Mechanics engine of Akantu Dependencies AKANT
80. ted in Figure 115 The term 0T 06 of equation after the necessary derivation can expressed as oT AT OT T 5125 s s t 2 m 4 79 where Ic p 6 5 1 o T 6 4 RS if gt max 4 80 if lt max n gt 0 45 Adding a New Constitutive Law There are several constitutive laws in Akantu as described in the previous Section 4 4 It is also possible to use a user defined material for the simulation These materials are referred to as local materials since they are local to the example of the user and not part of the Akantu library To define a new local material two files material_XXX hh and material_XXX cc have to be provided where XXX is the name of the new material The header file material_XXX hh defines 36 CHAPTER 4 SOLID MECHANICS MODEL Traction pas 0 fi L l l l 0 0 5 1 5 2 2 5 3 Opening m Figure 4 15 Exponential cohesive law the interface of your custom material Its implementation is provided in the material_XXX cc The new law must inherit from the Material class or any other existing material class It is therefore necessary to include the interface of the parent material in the header file of your local material and indicate the inheritance in the declaration of the class y VEREIN A IR A a Rer ME AI E IR O ES NOS NT RECO CEU AUN lA s include material hh Va DE NE A O E NE WO aE e E EN Te ma ARS all ila ER he t ome NO VE EN EU TS NES 24 ifndef __AKANTU_MAT
81. the velocity and acceleration fields at t 1 At Un 1 UP 41 zo 4 25 ii n O 4 26 The use of an explicit time integration scheme is illustrated by the example examples explicit explicit_dynamic cc This example models the propagation of a wave in a steel beam The beam and the applied displacement in the x direction are shown in Figure 4 8 The length and height of the beam are L 10m and h 1m respectively The material is linear elastic homogeneous and isotropic density 7800 kg m Young s modulus 210 GPa and Poisson s ratio 0 3 The imposed displacement follow a Gaussian function with a maximum amplitude of A 0 01m The potential kinetic and total energies are computed The safety factor is equal to 0 8 26 CHAPTER 4 SOLID MECHANICS MODEL u gt Y Figure 4 8 Numerical setup 44 Constitutive Laws In order to compute an element s response to deformation one needs to use an appropriate constitutive relationship The constitutive law is used to compute the element s stresses from the element s strains In the finite element discretization the constitutive formulation is applied to every quadrature point of each element When the implicit formulation is used the tangent matrix has to be computed The chosen materials for the simulation have to be specified in the mesh file or as an alter native they can be assigned using the element_material vector For every material a
82. this case the element insertion can be limited to a given range thanks to the method setLimit The first input parameter of this method indicates the direction while the other two indicate the extreme values of the range 0 9 1 1 In order to compute the energies the same functions illustrated for dynamics in the last section can be used 4 2 Static Analysis The SolidMechanicsModel class can handle different analysis methods the first one being pre sented is the static case In this case the equation to solve is Ku f 4 5 ext where K is the global stiffness matrix u the displacement vector and f the vector of external forces applied to the system ext To solve such a problem the static solver of the SolidMechanicsModel object is used First a model has to be created and initialized To create the model a mesh which can be read from a file is needed as explained in Section Once an instance of a SolidMechanicsModel is obtained the easiest way to initialize it is to use the initFull method by giving the SolidMechanicsModelOptions These options specify the type of analysis to be performed and whether the materials should be initialized or not SolidMechanicsModel model mesh model initFull SolidMechanicsModelOptions static 4 2 STATIC ANALYSIS 21 Here a static analysis is chosen by passing the argument _static to the method By default the Boolean for no initialization of the materials is set to false so tha
83. time the only partitioner available is MeshPartitionScotch which implements the function partitionate using the Scotch program This is achieved by the following code Mesh mesh spatial dimension MeshPartition partition NULL if prank 0 mesh read my mesh msh partition new MeshPartitionScotch mesh spatial dimension partition partitionate psize Note Only the processor of rank O should load the mesh file to partition it Nevertheless the Mesh object must by declared for all processors since the mesh distribution will store mesh pieces to that object 8 3 Distributing Mesh Partitions The distribution of the mesh is done automatically by the SolidMechanicsModel through the initParallel method Thus after creating a SolidMechanicsModel with our mesh as the initial parameter the initParallel method must be called receiving the partition as a parameter SolidMechanicsModel model mesh model initParallel partition After that point everything remains as in the sequential case from the user point of view This allows the user to care only about his simulation without concern for the parallelism An example of an explicit dynamic 2D bar in compression in a parallel context can be found in examples parallel_2d 84 Launching a Parallel Program Using MPla parallel run can be launched from a shell using the command mpirun np procs program name parameterl parameter2 Chapter 9 Contact 9 1 Impli
84. tion_local or in the global system of coordinates _bft_traction static void lin_load double position double load Real normal UInt surface_id memset load 0 sizeof Real 3 load 1 position 0 position 0 250 int main int argc char argv model computeForcesFromFunction lt _bernoulli_beam_2 gt lin_load 5 2 STATIC ANALYSIS 43 _bft_traction_local 5 2 Static Analysis The StructuralMechanicsModel class can perform static analyses of structures In this case the equation to solve is the same as for the SolidMechanicsModel used for static analyses Ku fext 5 1 where K is the global stiffness matrix u the generalized displacement vector and fe the vector of generalized external forces applied to the system To solve such a problem the static solver of the StructuralMechanicsModel object is used First a model has to be created and initialized StructuralMechanicsModel model mesh model initFullQ o model initFull initializes all internal vectors to zero Once the model is created and initialized the boundary conditions can be set as explained in Section Boundary conditions will prescribe the external forces or moments for the free degrees of freedom fext and displacements or rotations for the others To completely define the system represented by equation 5 1 the global stiffness matrix K must be assembled model assembleStiffnessMatrix The computation of the static equilibrium i
85. to make dissipated fracture energy path dependent The concept of free potential energy is dropped and a new independent parameter x 34 CHAPTER 4 SOLID MECHANICS MODEL T a Linear b Bilinear Figure 4 14 Irreversible cohesive laws for explicit simulations is introduced Gon K 4 69 Ga 4 69 where Geg and Gey are the necessary works of separation per unit area to open completely a cohesive zone under mode I and mode II respectively Their model yields to the following equation for cohesive tractions T in case of crack opening 2 Oc gt r gt 1 T 1 4 7 T Ene ann i I 3 4 70 where o is the material strength along the fracture the critical effective displacement after which cohesive tractions are zero complete decohesion A and A are the tangential and normal components of the opening displacement vector A respectively The parameter is a weight that indicates how big the tangential opening contribution is The effective opening displacement is 2 Vai AZ 4 71 In case of unloading or reloading lt max tractions are calculated as Oc T A 1 4 72 E s 5c 4 72 2 en 1 max 4 73 K max c so that they vary linearly between the origin and the maximum attained tractions As shown in Figure in this law the dissipated and reversible energies are 1 Ediss 59c max 4 74 1 E To 4 75 Moreover a dama
86. ver it is recommended to impose a time step that is smaller than the stable time step for instance by multiplying the stable time step by a safety factor smaller than one const Real safety time factor 0 8 Real applied time step time step safety time factor model setTimeStep applied time step The initial displacement and velocity fields are by default equal to zero if not given specifically by the user see 4 1 1 Like in implicit dynamics a time loop is used in which the displacement velocity and ac celeration fields are updated at each time step The values of these fields are obtained from the Newmark f equations with f 1 2 and a 0 In Akantu these computations at each time step are invoked by calling the function solveStep for UInt s 1 s 1 applied time step lt total time s model solveStep The method solveStep wraps the four following functions o model explicitPred allows to compute the displacement field at t 1 and a part of the velocity field at t 1 denoted by P 1 which will be used later in the method model explicitCorr The equations are sad Un 1 Un At n SY 4 22 Pn n Atin 4 23 o model updateResidual and model updateAcceleration compute the acceleration increment d m SAC West 4 24 intn 1 Note The internal force f 4 is computed from the displacement u 44 based on the constitutive law o model explicitCorr computes
87. ving an initial temperature of 100 K every where but a none centered hot point maintained at 300K Figure 6 1 presents the geometry of this case The material used is a linear fictitious elastic material with a density of 8940 kg m a conductivity of 401 Wm K and a specific heat capacity of 385 K kg The time step used is 0 12 s Figure 6 1 Initial temperature field left and after 15000 time steps 30 minutes right The lines represent iso surfaces Chapter 7 Input Output 7 1 Generic data In this chapter we address ways to get the internal data in human readable formats The models in Akantu handle data associated to the mesh but this data can be split into several Arrays For example the data stored per element type in a ElementTypeMapArray is composed of as many Arrays as types in the mesh In order to get this data in a visualization software the models contain a object to dump VTK files These files can be visualized in software such as ParaView 15 Visit or Mayavi 17 The internal dumper of the model can be configured to specify which data fields are to be written This is done with the addDumpField method By default all the files are generated in a folder called paraview model setBaseName output prefix for all generated files model addDumpField displacement model addDumpField stress model dump The fields are dumped with the number of components of the memory For example i

Akantu User's Guide - Computational Solid Mechanics Laboratory

Contents

Download Pdf Manuals

Related Search

Related Contents