Two inclusions mechanical analysis

Here we present a study case of a thick plate with 2 inclusions. One softer and the other stiffer than the base material. We then compute the energy deformation on one inclusion.

A Plate with two inclusion.

Fig. 2 Example mechanical simulation and post processing

  1
  2from BasicTools.FE.UnstructuredFeaSym import UnstructuredFeaSym
  3
  4#Main class to
  5problem = UnstructuredFeaSym()
  6
  7# the mecanical problem
  8from BasicTools.FE.SymPhysics import MecaPhysics
  9mecaPhysics = MecaPhysics()
 10
 11# Definition of the degree of the spaces [1 or 2]
 12mecaPhysics.SetSpaceToLagrange(P=1)
 13
 14# add weak form terms to the tanget matrix
 15mecaPhysics.AddBFormulation( "Bulk",mecaPhysics.GetBulkFormulation(1.0,0.3)  )
 16mecaPhysics.AddBFormulation( "Inclusion1",mecaPhysics.GetBulkFormulation(5.0,0.3)  )
 17youngModulusInclusionII = 0.5
 18mecaPhysics.AddBFormulation( "Inclusion2",mecaPhysics.GetBulkFormulation(0.5,0.3)  )
 19
 20# add weak form term to the rhs
 21mecaPhysics.AddLFormulation( "Right", mecaPhysics.GetForceFormulation([1,0,0],1)  )
 22
 23problem.physics.append(mecaPhysics)
 24
 25# the boundary conditions
 26from BasicTools.FE.KR.KRBlock import KRBlockVector
 27dirichlet = KRBlockVector()
 28dirichlet.AddArg("u").On('Left').Fix0().Fix1().Fix2().To()
 29
 30problem.solver.constraints.AddConstraint(dirichlet)
 31
 32# Read The mesh
 33from BasicTools.IO.GmshReader import ReadGmsh
 34mesh = ReadGmsh("TwoInclusions.msh")
 35mesh.ConvertDataForNativeTreatment()
 36print(mesh)
 37
 38problem.SetMesh(mesh)
 39
 40# we compute the numbering (maping from the mesh to the linear system)
 41problem.ComputeDofNumbering()
 42
 43from BasicTools.Helpers.Timer import Timer
 44with Timer("Assembly "):
 45    k,f = problem.GetLinearProblem()
 46
 47#compute the constraints to add to the system
 48problem.ComputeConstraintsEquations()
 49
 50
 51with Timer("Solve"):
 52    problem.Solve(k,f)
 53
 54problem.PushSolutionVectorToUnkownFields()
 55
 56from BasicTools.FE.Fields.FieldTools import GetPointRepresentation
 57# Recover a point representation of the displacement
 58problem.mesh.nodeFields["sol"] = GetPointRepresentation(problem.unkownFields)
 59
 60from BasicTools.FE.Fields.FEField import FEField
 61#Creation of a fake fields to export the rhs member
 62rhsFields = [ FEField(mesh=mesh,space=None,numbering=problem.numberings[i]) for i in range(3) ]
 63from BasicTools.FE.Fields.FieldTools import VectorToFEFieldsData
 64VectorToFEFieldsData(f,rhsFields)
 65problem.mesh.nodeFields["RHS"] = GetPointRepresentation(rhsFields)
 66
 67print("Done solve")
 68print("Compute of the strain energy only on the second inclusion (integral in each element) ")
 69
 70import BasicTools.FE.SymWeakForm as SWF
 71from BasicTools.FE.MaterialHelp import HookeIso
 72from BasicTools.Containers.Filters import ElementFilter
 73
 74symdep = SWF.GetField("u",3)
 75K = HookeIso(youngModulusInclusionII,0.3,dim=3)
 76symCellData = SWF.GetField("cellData",1)
 77symCellDataT = SWF.GetTestField("cellData",1)
 78
 79print("Post process")
 80
 81EnerForm = SWF.ToVoigtEpsilon(SWF.Strain(symdep)).T*K*SWF.ToVoigtEpsilon(SWF.Strain(symdep))*symCellDataT
 82
 83print("Post process Eval")
 84ff = ElementFilter(mesh=problem.mesh, tag="Inclusion2")
 85from BasicTools.FE.DofNumbering import  ComputeDofNumbering
 86from BasicTools.FE.Spaces.FESpaces import LagrangeSpaceP0
 87from BasicTools.FE.Integration import IntegrateGeneral
 88p0Numbering = ComputeDofNumbering(mesh,LagrangeSpaceP0)
 89
 90energyDensityField = FEField(name="cellData",mesh=problem.mesh,numbering=p0Numbering,space=LagrangeSpaceP0)
 91
 92m,energyDensity = IntegrateGeneral(mesh=problem.mesh, wform=EnerForm,  constants={},
 93                        fields=problem.unkownFields, unkownFields = [ energyDensityField ], elementFilter=ff)
 94print("energyDensity",energyDensity)
 95energyDensityField.data = energyDensity
 96
 97problem.mesh.elemFields["Energy"] = energyDensity
 98
 99import numpy as np
100print("Strain energy on the second inclusion:", np.sum(energyDensity) )
101from BasicTools.IO import XdmfWriter as XW
102writer = XW.XdmfWriter('TwoInclusions_Output.xdmf')
103writer.SetHdf5(False)
104writer.Open()
105writer.Write(mesh,PointFields=list(mesh.nodeFields.values()), PointFieldsNames=list(mesh.nodeFields.keys()),
106                CellFields=list(mesh.elemFields.values()), CellFieldsNames=list(mesh.elemFields.keys()))
107writer.Close()
108
109