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. Then we compute the energy deformation on one inclusion.

A Plate with two inclusion.

Fig. 1 Example of mechanical simulation and post processing

  1
  2from Muscat.FE.UnstructuredFeaSym import UnstructuredFeaSym
  3
  4#Main class to
  5problem = UnstructuredFeaSym()
  6
  7# the mecanical problem
  8from Muscat.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 Muscat.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 Muscat.IO.GmshReader import ReadGmsh
 34mesh = ReadGmsh("TwoInclusions.msh")
 35print(mesh)
 36
 37problem.SetMesh(mesh)
 38
 39# we compute the numbering (maping from the mesh to the linear system)
 40problem.ComputeDofNumbering()
 41
 42from Muscat.Helpers.Timer import Timer
 43with Timer("Assembly "):
 44    k,f = problem.GetLinearProblem()
 45
 46#compute the constraints to add to the system
 47problem.ComputeConstraintsEquations()
 48
 49
 50with Timer("Solve"):
 51    problem.Solve(k,f)
 52
 53problem.PushSolutionVectorToUnknownFields()
 54
 55from Muscat.FE.Fields.FieldTools import GetPointRepresentation
 56# Recover a point representation of the displacement
 57problem.mesh.nodeFields["sol"] = GetPointRepresentation(problem.unknownFields)
 58
 59from Muscat.FE.Fields.FEField import FEField
 60#Creation of a fake fields to export the rhs member
 61rhsFields = [ FEField(mesh=mesh,space=None,numbering=problem.numberings[i]) for i in range(3) ]
 62from Muscat.FE.Fields.FieldTools import VectorToFEFieldsData
 63VectorToFEFieldsData(f,rhsFields)
 64problem.mesh.nodeFields["RHS"] = GetPointRepresentation(rhsFields)
 65
 66print("Done solve")
 67print("Compute of the strain energy only on the second inclusion (integral in each element) ")
 68
 69import Muscat.FE.SymWeakForm as SWF
 70from Muscat.FE.MaterialHelp import HookeIso
 71from Muscat.Containers.Filters.FilterObjects import ElementFilter
 72
 73symdep = SWF.GetField("u",3)
 74K = HookeIso(youngModulusInclusionII,0.3,dim=3)
 75symCellData = SWF.GetField("cellData",1)
 76symCellDataT = SWF.GetTestField("cellData",1)
 77
 78print("Post process")
 79
 80EnerForm = SWF.ToVoigtEpsilon(SWF.Strain(symdep)).T*K*SWF.ToVoigtEpsilon(SWF.Strain(symdep))*symCellDataT
 81
 82print("Post process Eval")
 83ff = ElementFilter(eTag="Inclusion2")
 84from Muscat.FE.DofNumbering import  ComputeDofNumbering
 85from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP0
 86from Muscat.FE.Integration import IntegrateGeneral
 87p0Numbering = ComputeDofNumbering(mesh,LagrangeSpaceP0)
 88
 89energyDensityField = FEField(name="cellData",mesh=problem.mesh,numbering=p0Numbering,space=LagrangeSpaceP0)
 90
 91m,energyDensity = IntegrateGeneral(mesh=problem.mesh, wform=EnerForm,  constants={},
 92                        fields=problem.unknownFields, unknownFields = [ energyDensityField ], elementFilter=ff)
 93print("energyDensity",energyDensity)
 94energyDensityField.data = energyDensity
 95
 96problem.mesh.elemFields["Energy"] = energyDensity
 97
 98import numpy as np
 99print("Strain energy on the second inclusion:", np.sum(energyDensity) )
100# To visualize this xdmf file you can use ParaView (downloadable from https://www.paraview.org/)
101from Muscat.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
109print("Done")
110