mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
Jukka Aho
215 lines
8.6 KiB
Julia
215 lines
8.6 KiB
Julia
# This file is a part of JuliaFEM.
|
|
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
|
|
|
# # 2D Hertz contact problem
|
|
|
|
# 
|
|
|
|
# In the example, a cylinder is pressed agains block with a force of 35 kN.
|
|
# A similar example can be found from NAFEMS report FENET D3613 (advanced
|
|
# finite element contact benchmarks).
|
|
#
|
|
# Solution for maximum pressure ``p_0`` and contact radius ``a`` is
|
|
# ```math
|
|
# p_{0} = \sqrt{\frac{FE}{2\pi R}}, \\
|
|
# a = \sqrt{\frac{8FR}{\pi E}},
|
|
# ```
|
|
# where
|
|
#
|
|
# ```math
|
|
# E = \frac{2E_{1}E_{2}}{E_{2}\left(1-\nu_{1}^{2}\right)+E_{1}\left(1-\nu_{2}^{2}\right)}.
|
|
# ```
|
|
#
|
|
# Substituting values, one gets accurate solution to be ``p_0 = 3585 \;\mathrm{MPa}`` and
|
|
# ``a = 6.21 \;\mathrm{mm}``.
|
|
|
|
using JuliaFEM, LinearAlgebra
|
|
|
|
# Simulation starts by reading the mesh. Model is constructed and meshed using
|
|
# SALOME, thus mesh format is .med. Mesh type is quite simple structure,
|
|
# containing things like `mesh.nodes`, `mesh.elements` and so on. Keep on mind,
|
|
# that Mesh contains only standard Julia types and we think it as a structure
|
|
# helping us to construct elements needed in simulation. In principle, we don't
|
|
# need to use `Mesh` in simulation anyway if we figure some other way to define
|
|
# the geometry for elements.
|
|
|
|
datadir = abspath(joinpath(pathof(JuliaFEM), "..", "..", "examples", "2d_hertz_contact"))
|
|
meshfile = joinpath(datadir, "hertz_2d_full.med")
|
|
mesh = aster_read_mesh(meshfile)
|
|
for (elset_name, element_ids) in mesh.element_sets
|
|
nel = length(element_ids)
|
|
println("Element set $elset_name contains $nel elements.")
|
|
end
|
|
for (nset_name, node_ids) in mesh.node_sets
|
|
nno = length(node_ids)
|
|
println("Node set $nset_name contains $nno nodes.")
|
|
end
|
|
nnodes = length(mesh.nodes)
|
|
println("Total number of nodes in mesh: $nnodes")
|
|
nelements = length(mesh.elements)
|
|
println("Total number of elements in mesh: $nelements")
|
|
|
|
# Next, define two bodies. Technically, we could have only one problem and add
|
|
# elements from both bodies to the same problem, but defining two different
|
|
# problems is recommended for clarity. Plain strain assumption is used.
|
|
# To make clear what is happening here: we first create a set of elements
|
|
# (elements are in vector called `upper_elements`), then we define new
|
|
# problem which type is `Elasticity`, give it some meaningful name (this time
|
|
# `cylinder`), and last value 2 means that problems does have two degrees of
|
|
# freedom per node.
|
|
|
|
upper_elements = create_elements(mesh, "CYLINDER")
|
|
update!(upper_elements, "youngs modulus", 70.0e3)
|
|
update!(upper_elements, "poissons ratio", 0.3)
|
|
upper = Problem(Elasticity, "cylinder", 2)
|
|
upper.properties.formulation = :plane_strain
|
|
add_elements!(upper, upper_elements)
|
|
|
|
lower_elements = create_elements(mesh, "BLOCK")
|
|
update!(lower_elements, "youngs modulus", 210.0e3)
|
|
update!(lower_elements, "poissons ratio", 0.3)
|
|
lower = Problem(Elasticity, "block", 2)
|
|
lower.properties.formulation = :plane_strain
|
|
add_elements!(lower, lower_elements)
|
|
|
|
# Next we define some boundary conditions: creating "boundary" problems goes
|
|
# in the same way than defining "field" problems, the only difference is that
|
|
# we add extra argument giving what field are we tring to fix. This time,
|
|
# we have 2 dofs / node and we fix displacement in direction 2.
|
|
|
|
bc_fixed_elements = create_elements(mesh, "FIXED")
|
|
update!(bc_fixed_elements, "displacement 2", 0.0)
|
|
bc_fixed = Problem(Dirichlet, "fixed", 2, "displacement")
|
|
add_elements!(bc_fixed, bc_fixed_elements)
|
|
|
|
# Defining symmetry boundary condition goes with the same idea
|
|
|
|
bc_sym_23_elements = create_elements(mesh, "SYM23")
|
|
update!(bc_sym_23_elements, "displacement 1", 0.0)
|
|
bc_sym_23 = Problem(Dirichlet, "symmetry line 23", 2, "displacement")
|
|
add_elements!(bc_sym_23, bc_sym_23_elements)
|
|
|
|
# Next we define point load. To define that, we first need to find some node
|
|
# near the top of cylinder, using function `find_nearest_node`. Then we create
|
|
# a new problem, again of type Elasticity. Like told already, we don't need to
|
|
# use `Mesh` if we have some other procedure to define the geometry of the
|
|
# element (and it's connectivity, of course). So we can directly create an
|
|
# element of type `Poi1`, meaning 1-node point element, update it's geometry
|
|
# and apply 35.0e3 kN load in negative y-direction:
|
|
|
|
nid = find_nearest_node(mesh, [0.0, 100.0])
|
|
load = Problem(Elasticity, "point load", 2)
|
|
load.properties.formulation = :plane_strain
|
|
load.elements = [Element(Poi1, [nid])]
|
|
update!(load.elements, "geometry", mesh.nodes)
|
|
update!(load.elements, "displacement traction force 2", -35.0e3)
|
|
|
|
# Next, we define another boudary problem, this time the type of problem is
|
|
# Contact2D, which is a mortar contact formulation for two dimensions.
|
|
# Elements are added using `add_slave_elements!` and `add_master_elements!`.
|
|
# Problems, in general, can have some properties defined, like the formulation
|
|
# in `Elasticity` (we also have `:plane_stress`). For contact, we need to swap
|
|
# normal direction for meshes created by SALOME because in Code Aster, element
|
|
# orientation is defined opposite to what is used in ABAQUS, and in JuliaFEM in
|
|
# general we follow the same conventions what are used in ABAQUS.
|
|
|
|
contact = Problem(Contact2D, "contact", 2, "displacement")
|
|
contact.properties.rotate_normals = true
|
|
contact_slave_elements = create_elements(mesh, "BLOCK_TO_CYLINDER")
|
|
contact_master_elements = create_elements(mesh, "CYLINDER_TO_BLOCK")
|
|
add_master_elements!(contact, contact_master_elements)
|
|
add_slave_elements!(contact, contact_slave_elements)
|
|
|
|
# After all problems are defined, we define some `Analysis`, which can be e.g.
|
|
# static analysis, dynamic analysis, modal analysis, linear perturbation
|
|
# analysis and so on. Here, the analysis type is `Nonlinear`, which is nonlinear
|
|
# quasistatic analysis. In the same manner as we do `add_elements!` to add
|
|
# elements to `Problem`, we use `add_problems!` to add problems to analysis.
|
|
# Because we are not restricted to some particular input and output formats,
|
|
# we "connect" a `ResultsWriter` to our analysis, this time we want to visualize
|
|
# results using ParaView, thus we write our results to Xdmf format, which uses
|
|
# well defined standards XML and HDF to store model data.
|
|
|
|
analysis = Analysis(Nonlinear)
|
|
add_problems!(analysis, upper, lower, bc_fixed, bc_sym_23, load, contact)
|
|
xdmf = Xdmf("2d_hertz_results"; overwrite=true)
|
|
add_results_writer!(analysis, xdmf)
|
|
|
|
# In last part, we run the analysis.
|
|
|
|
run!(analysis)
|
|
close(xdmf)
|
|
|
|
# # Results
|
|
|
|
# Results are stored in `2d_hertz_results.xmf` and `2d_hertz_results.h5` for
|
|
# visual inspection. We can also postprocess results programmatically because
|
|
# we are inside a real scripting / programming environment all the time. For
|
|
# example, we can integrate the resultant force in normal and tangential direction
|
|
# in contact surface to validate our result.
|
|
|
|
Rn = 0.0
|
|
Rt = 0.0
|
|
time = 0.0
|
|
for sel in contact_slave_elements
|
|
for ip in get_integration_points(sel)
|
|
global Rn, Rt
|
|
w = ip.weight*sel(ip, time, Val{:detJ})
|
|
n = sel("normal", ip, time)
|
|
t = sel("tangent", ip, time)
|
|
la = sel("lambda", ip, time)
|
|
Rn += w*dot(n, la)
|
|
Rt += w*dot(t, la)
|
|
end
|
|
end
|
|
|
|
println("2d hertz contact resultant forces: Rn = $Rn, Rt = $Rt")
|
|
|
|
using Test
|
|
@test isapprox(Rn, 35.0e3)
|
|
@test isapprox(Rt, 0.0)
|
|
|
|
# Visualization of the results can be done using ParaView:
|
|
# 
|
|
|
|
# For optimization loops, we want to programmatically find, for example, maximum
|
|
# contact pressure. We can, for example, get all the values in nodes:
|
|
|
|
lambda = contact("lambda", time)
|
|
normal = contact("normal", time)
|
|
p0 = 0.0
|
|
p0_acc = 3585.0
|
|
for (nid, n) in normal
|
|
lan = dot(n, lambda[nid])
|
|
println("$nid => $lan")
|
|
global p0
|
|
p0 = max(p0, lan)
|
|
end
|
|
p0 = round(p0, digits=2)
|
|
rtol = round(norm(p0-p0_acc)/max(p0,p0_acc)*100, digits=2)
|
|
println("Maximum contact pressure p0 = $p0, p0_acc = $p0_acc, rtol = $rtol %")
|
|
|
|
# To get rough approximation where does the contact open, we can find the element
|
|
# from slave contact surface, where contact pressure is zero in the other node
|
|
# and something nonzero in the other node.
|
|
|
|
a_rad = 0.0
|
|
for element in contact_slave_elements
|
|
la1, la2 = element("lambda", time)
|
|
p1, p2 = norm(la1), norm(la2)
|
|
a, b = isapprox(p1, 0.0), isapprox(p2, 0.0)
|
|
if (a && !b) || (b && !a)
|
|
X1, X2 = element("geometry", time)
|
|
println("Contact opening element geometry: X1 = $X1, X2 = $X2")
|
|
println("Contact opening element lambda: la1 = $la1, la2 = $la2")
|
|
x11, y11 = X1
|
|
x12, y12 = X2
|
|
global a_rad
|
|
a_rad = 1/2*abs(x11+x12)
|
|
break
|
|
end
|
|
end
|
|
println("Contact radius: $a_rad")
|
|
|
|
# This example briefly described some of the core features of JuliaFEM.
|