Author(s): Jukka Aho
Abstract: Solving elasticity equations using JuliaFEM.
Given function spaces \begin{align} \boldsymbol{\mathcal{U}} & =\left\{ \boldsymbol{u}\in H^{1}\left(\Omega\right)|\boldsymbol{u}\left(\boldsymbol{X},t\right)=\hat{\boldsymbol{u}}\left(\boldsymbol{X},t\right)\text{ on }\Gamma_{\mathrm{u}}\right\} ,\\ \boldsymbol{\mathcal{V}} & =\left\{ \delta\boldsymbol{u}\in H^{1}\left(\Omega\right)|\delta\boldsymbol{u}\left(\boldsymbol{X}\right)=0\text{ on }\Gamma_{\mathrm{u}}\right\} , \end{align} find $\boldsymbol{u}\in\boldsymbol{\mathcal{U}}$ such that \begin{equation} \delta\mathcal{W}:=\int_{\Omega_{0}}\rho_{0}\ddot{\boldsymbol{u}}\cdot\delta\boldsymbol{u}\,\mathrm{d}V_{0}+\int_{\Omega_{0}}\boldsymbol{S}:\delta\boldsymbol{E}\,\mathrm{d}V_{0}-\int_{\Omega_{0}}\hat{\boldsymbol{b}}_{0}\cdot\delta\boldsymbol{u}\,\mathrm{d}V_{0}-\int_{\Gamma_{\sigma}}\hat{\boldsymbol{t}}_{0}\cdot\delta\boldsymbol{u}\,\mathrm{d}A_{0} =0 \qquad\forall\delta\boldsymbol{u}\in\boldsymbol{\mathcal{V}} \end{equation}
https://en.wikipedia.org/wiki/Strain_energy_density_function
Saint-Venant-Kirchhoff model https://en.wikipedia.org/wiki/Hyperelastic_material \begin{equation} \psi\left(\mathbf{E}\right)=\frac{\lambda}{2}\left[\mbox{tr}\left(\mathbf{E}\right)\right]^{2}+\mu\mbox{tr}\left(\mathbf{E}^2\right) \end{equation}
neo-Hookean material https://en.wikipedia.org/wiki/Neo-Hookean_solid \begin{equation} \psi=\frac{\mu}{2}\left(I_{c}-3\right)-\mu\ln\left(J\right)+\frac{\lambda}{2}\ln\left(J\right)^{2} \end{equation}
In [1]:
using JuliaFEM.API: Model, Material, add_material!, LoadCase, add_element_set!, ForceBC, DisplacementBC
using JuliaFEM.API: add_boundary_condition!, add_load_case!, add_solver!
using JuliaFEM.Interfaces: solve!
using JuliaFEM.Preprocess: parse_abaqus
In [2]:
@time begin
# Linear models
#model = open(JuliaFEM.Core.parse_abaqus, "../geometry/piston/piston_8789_P1.inp")
#model = open(JuliaFEM.Core.parse_abaqus, "../geometry/piston/piston_16436_P1.inp")
#model = open(JuliaFEM.Core.parse_abaqus, "../geometry/piston/piston_27343_P1.inp")
#model = open(JuliaFEM.Core.parse_abaqus, "../geometry/piston/piston_45510_P1.inp")
#model = open(JuliaFEM.Core.parse_abaqus, "../geometry/piston/piston_75470_P1.inp")
#model = open(JuliaFEM.Core.parse_abaqus, "../geometry/wrench/wrench_128903_P1.inp")
# Quadratic models
#model = open(JuliaFEM.Core.parse_abaqus, "../geometry/piston/piston_19611_P2.inp")
#model = open(JuliaFEM.Core.parse_abaqus, "../geometry/piston/piston_55950_P2.inp")
#model = open(JuliaFEM.Core.parse_abaqus, "/Temp/piston_107168_P2.inp")
#model = open(JuliaFEM.Core.parse_abaqus, "../geometry/piston/piston_345757_P2.inp")
mesh = open(parse_abaqus, "/Temp/piston_107168_P2.inp")
model = Model("piston model", mesh)
info("model loaded.")
end
In [3]:
mat = Material("steel")
mat["youngs modulus"] = 210.0e3
mat["poissons ratio"] = 0.3
add_material!(model, "PISTON", mat)
Out[3]:
In [4]:
simulation = LoadCase(:ElasticityProblem)
# simulation = Simulation(:ElasticityProblem, "elasticity equations")
# add_element_set!(simulation, "PISTON")
# add_boundary_condition!(simulation, bc1, bc2)
# add_simulation!(model, simulation)
# add_solver!(simulation, :DirectSolver)
# solve!(model, "elasticity equations", 0.0)
# solve!(model, "heat equations", 0.0)
# solve!(model, "solve equations 3", 0.0)
add_element_set!(simulation, "PISTON")
Out[4]:
In [5]:
traction = Vector{Float64}[[0.0, 0.0, -10.0] for i in 1:6]
bc1 = ForceBC("BC1", "displacement traction force" => traction)
bc2 = DisplacementBC("BC2", "displacement" => 0.0)
# JuliaFEM.API.add_boundary_condition!(prob, bc1, bc2)
add_boundary_condition!(simulation, bc1)
add_boundary_condition!(simulation, bc2)
Out[5]:
In [6]:
add_solver!(simulation, :DirectSolver)
add_load_case!(model, "elasticity equations", simulation)
Out[6]:
In [7]:
solver = JuliaFEM.Interfaces.get_solver(model, "elasticity equations", 0.0);
In [11]:
celements = solver.field_problems[1].elements;
celements[1]
Out[11]:
In [12]:
surface_elements = filter((e) -> isa(e, JuliaFEM.Core.Element{JuliaFEM.Core.Tri6}), celements)
info("$(length(surface_elements)) surface elements")
In [13]:
length(celements)
Out[13]:
In [7]:
solve!(model, "elasticity equations", 0.0)
Out[7]:
In [20]:
xdoc, xmodel = JuliaFEM.Postprocess.xdmf_new_model()
temporal_collection = JuliaFEM.Postprocess.xdmf_new_temporal_collection(xmodel)
grid = JuliaFEM.Postprocess.xdmf_new_grid(temporal_collection; time=0)
Out[20]:
Save geometry to xdmf file
In [21]:
nnodes = length(model["nodes"])
info("number of nodes in model: $nnodes")
X = zeros(3, nnodes)
for nid in keys(model["nodes"])
X[:, perm[nid]] = model["nodes"][nid]
end
nelements = length(model["elsets"]["PISTON"])
info("Number of elements in model: $nelements")
elmap = zeros(Int64, 5, nelements)
#elmap[1,:] = 0x0026
elmap[1,:] = 0x6
for (i, elid) in enumerate(model["elsets"]["PISTON"])
elmap[2:end,i] = Int64[perm[nid] for nid in model["elements"][elid]]
end
In [22]:
JuliaFEM.Postprocess.xdmf_new_mesh(grid, X, elmap)
Out[22]:
Save nodal data to model
In [23]:
using JuliaFEM.Core: get_connectivity
u = zeros(3, nnodes)
for element in values(elements)
isa(element, Element{Tet4}) || continue
connectivity = get_connectivity(element)
field = element["displacement"](0.0)
for (i, nid) in enumerate(connectivity)
u[:, nid] = field[i]
end
end
size(u)
Out[23]:
In [24]:
JuliaFEM.Postprocess.xdmf_new_field(grid, "Displacement", "nodes", u)
Out[24]:
In [25]:
JuliaFEM.Postprocess.xdmf_save_model(xdoc, "/tmp/piston_8789_P1.xmf")
Out[25]:
In [26]:
u[:, 1:10]
Out[26]:
In [27]:
minimum(u, 2)'
Out[27]:
In [28]:
maximum(u, 2)'
Out[28]:
In [30]:
piston_8789_P1_solution_norm = 38.52232721814439
Out[30]:
In [31]:
@assert isapprox(norm(vec(u)), piston_8789_P1_solution_norm)
In [ ]: