1

I'm following this tutorial in order to try and do an FEA of a model.msh that I have to see how it would deform given different external forces in different places.

There they define the weak form as

a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )*dΩ
l(v) = 0

and they state "The linear form is simply l(v) = 0 since there are not external forces in this example."

As mentioned, I would like to analyse the different deformation that different external forces would cause on my model, but I can't seem to find anywhere an example of this. Could someone help me on defining this linear form for external forces different than 0?

Thanks.

guin0x
  • 337
  • 2
  • 10

1 Answers1

1

Maybe this helps you out. It was written in a hurry, so please do not mind if you encounter spelling mistakes or other beauty issues :)

# define where the output shall go
output_Path ="Output/3_Gridap/1_Lin_FEA/FE_8"
mkpath(output_Path)
output_Name ="pde_6"

using Gridap

# please load the model that is shown in: https://gridap.github.io/Tutorials/dev/pages/t001_poisson/
model = DiscreteModelFromFile("Code/Meshes/Data/possion.json")

# just in case you want to see the model using paraview
writevtk(model,"$(output_Path)/md")

order = 1
reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},order)

V0 = TestFESpace(model,reffe;
  conformity=:H1,

  # to see which elements belongs to "bottom" open the model which is saved through "writevtk(model,"$(output_Path)/md")"
  dirichlet_tags=["bottom"],
  
  # activate/deactivate the boundary conditions
  dirichlet_masks=[
        (true, true, true),       # clamp the bottom
        ])

# define displacement
clamping(x) = VectorValue(0.0,0.0,0.0)

U = TrialFESpace(V0,[clamping])
const E = 7e+7
const ν = 0.33
const λ = (E*ν)/((1+ν)*(1-2*ν))
const μ = E/(2*(1+ν))
σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε

degree = 2*order
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)

# Neumann boundary conditions
# here we define the surface on which we want an external force
Γ_Tr = BoundaryTriangulation(model,tags=["triangle"])
dΓ_Tr = Measure(Γ_Tr,degree)

# a force shall be applied on the y-direction
f_Tr(x) = VectorValue(0.0, 1e+6, 0.0)

# mass forces due to gravity, the value is set quite high, such that an impact can be seen
mass_Forces(x) = VectorValue(0.0, -1e+7, 0.0)

# Weak form
a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )*dΩ

l(v) = ∫( v ⋅ mass_Forces )* dΩ + ∫( v ⋅ f_Tr )* dΓ_Tr 

op = AffineFEOperator(a,l,U,V0)
uh = solve(op)

writevtk(Ω,"$(output_Path)/$(output_Name)",
      
cellfields=[
            "uh" => uh,
            "epsi" => ε(uh),
            "sigma" => σ∘ε(uh)])
Javed Was
  • 56
  • 1
  • 5
  • hey, thanks for that... what is the purpose of the mass_Forces here? – guin0x Jan 28 '23 at 14:50
  • mass forces are body forces, i.e., volume forces like the force due to the density and gravitation. These are the Dirichlet boundary conditions. The surface forces are the Neumann boundary conditions. – Javed Was Jan 29 '23 at 14:17