We impose the following boundary conditions. All components of the displacement vector are constrained to zero on the surface , which is marked in green in the figure. On the other hand, the first component of the displacement vector is prescribed to the value mm on the surface , which is marked in blue. No body or surface forces are included in this example. Formally, the PDE to solve is
The variable stands for the unknown displacement vector, the vector is the unit outward normal to the Neumann boundary and is the stress tensor defined as
where is the 2nd order identity tensor, and and are the Lamé parameters of the material. The operator is the symmetric gradient operator (i.e., the strain tensor). Here, we consider material parameters corresponding to aluminum with Young's modulus Pa and Poisson's ratio . From these values, the Lamé parameters are obtained as and .
We use a conventional Galerkin FE method with conforming Lagrangian FE spaces. For this formulation, the weak form is: find such that for all , where is the subset of functions in that fulfill the Dirichlet boundary conditions of the problem, whereas are functions in fulfilling on and on . The bilinear form of the problem is
On top of it's mesh-generation functionalities, Gridap provides a convenient ways to read and write discrete models from and to files. We import the model for this problem from GMSH
as follows:
using Gridap, GridapGmsh
using Gridap.Io
using DrWatson
msh_file_gmsh = projectdir("meshes/elasticity.msh")
model = GmshDiscreteModel(msh_file_gmsh)
This model contains the mesh and the physical tags of the model, which have been created directly through GMSH
. Another option would be to create the model using Gridap, then exporting to .json
format in the following way:
msh_file_json = projectdir("meshes/elasticity.json")
to_json_file(model,msh_file_json)
We can then load the saved model
model = DiscreteModelFromFile(msh_file_json)
You can easily inspect the generated discrete model in ParaView by writing it in vtk
format.
writevtk(model,datadir("elasticity_model"))
The previous line generates four different files model_0.vtu
, model_1.vtu
, model_2.vtu
, and model_3.vtu
containing the vertices, edges, faces, and cells present in the discrete model. Moreover, you can easily inspect which boundaries are defined within the model.
Open the resulting files with ParaView. Visualize the faces of the model and color them by each of the available fields. Identify the field names representing the boundaries and .
We will construct the vector-valued test FE space as follows:
order = 1
reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},order)
V0 = TestFESpace(model,reffe;
conformity=:H1,
dirichlet_tags=dirichlet_tags,
dirichlet_masks=dirichlet_masks)
The vector-valued interpolation is selected via the option valuetype=VectorValue{3,Float64}
, where we use the type VectorValue{3,Float64}
, which is the way Gridap represents vectors of three Float64
components.
In the next two exercises, we will fill in
the dirichlet_tags
optional argument to identify the Dirichlet regions, and
the dirichlet_masks
optional argument to specify which components of the displacement are constrained.
Let's go step-by-step.
Fill in the dirichlet_tags
using the tag names identified in Exercise 2 for and .
Hint: The general input format of dirichlet_tags
is a one-dimensional array of tag name strings ["tag_name_1",...,"tag_name_n"]
.
# Solution for exercise 2
Fill in the dirichlet_masks
to select the displacement components to constrain.
Hint: The general input format of dirichlet_masks
is a one-dimensional array with the same length and order as dirichlet_tags
. Each entry of dirichlet_masks
is a boolean tuple of length the problem dimension D
, i.e. (bool_x1,...,bool_xD)
. If bool_xi == true
, then the xi
component of the displacement is constrained, otherwise it is free. Beware of passing the dirichlet_masks
in the same order as the dirichlet_tags
.
Recall that we constrain only the first component on the boundary , whereas we constrain all components on .
# Solution for exercise 3
We can now instantiate the vector-valued test FE space.
order = 1
reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},order)
V0 = TestFESpace(model,reffe;
conformity=:H1,
dirichlet_tags=dirichlet_tags,
dirichlet_masks=dirichlet_masks)
Next, we construct the trial FE space. The Dirichlet boundary conditions must be described with two different functions, one for boundary and another one for .
Define the Dirichlet functions according to the problem statement.
Hint: The functions must be vector-valued with the format VectorValue(val_x1,...,val_xD)
.
# Solution for exercise 4
Define the trial FE space U
.
Hint: Pass the Dirichlet functions in a one-dimensional array, in the same order as dirichlet_tags
and dirichlet_masks
.
# Solution for exercise 5
We will now visually check the Dirichlet values are being correctly assigned on the Dirichlet boundaries.
Create a BoundaryTriangulation
of the Dirichlet boundaries.
# Solution for exercise 6
Next, we create a FE Function of U
with zero-valued free values.
vh = zero(U)
Plot vh
on the Dirichlet boundaries using writevtk
# Solution for exercise 7
The plot of the x-component of vh
should look as follows.
Why do the results at the contour of look weird?
The reason is because we are enforcing the Dirichlet values only in the interior of the Dirichlet regions. So, in the contour of the x-component of the DoF values must be zero. We can use the tag names of the model ending with _c
to enforce the Dirichlet conditions on the boundary closure too.
Read the details about how we define weak form, solve the problem and visualise the results.
The construction of the weak form needs to account for the constitutive law that relates strain and stress. The symmetric gradient operator is represented by the function ε
provided by Gridap (also available as symmetric_gradient
). However, function σ
representing the stress tensor is not predefined in the library and it has to be defined ad-hoc by the user, namely
const E = 70.0e9
const ν = 0.33
const λ = (E*ν)/((1+ν)*(1-2*ν))
const μ = E/(2*(1+ν))
σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε
Function σ
takes a strain tensor ε
(one can interpret this strain as the strain at an arbitrary integration point) and computes the associated stress tensor using the Lamé operator.
We build now the integration mesh and the corresponding measure
degree = 2*order
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)
From these objects and the constitutive law previously defined, we can write the weak form as follows
a(u,v) = ∫( (σ∘ε(u)) ⊙ ε(v) )dΩ
l(v) = 0
Note that we have composed function σ
with the strain field ε(u)
in order to compute the stress field associated with the trial function u
. The linear form is simply l(v) = 0
since there are not external forces in this example.
The remaining steps for solving the FE problem are rather standard.
op = AffineFEOperator(a,l,U,V0)
uh = solve(op)
Finally, we write the results to a file. Note that we also include the strain and stress tensors into the results file.
out_file = datadir("elasticity_sol")
writevtk(Ω,out_file,cellfields=["uh"=>uh,"epsi"=>ε(uh),"sigma"=>σ∘ε(uh)])
If you plot the deformed displacement in ParaView, it can be clearly observed (see next figure) that the surface is pulled in -direction and that the solid deforms accordingly.
Modify this model to solve a multi-material linear elasticity problem as done here.
Exercise done!
This exercise is adapted from the Gridap Tutorial 3: Linear elasticity.