mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-21 18:33:36 +00:00
minor changes to 2d contact formulation
This commit is contained in:
Jukka Aho
+10
-10
@@ -158,19 +158,19 @@ As a result element now have time invariant (variable) vector field "geometry" w
|
||||
|
||||
"""
|
||||
function update!(element::Element, field_name, data::Dict)
|
||||
element[field_name] = Field(data)
|
||||
#element[field_name] = [data[i] for i in get_connectivity(element)]
|
||||
#element[field_name] = Field(data)
|
||||
element[field_name] = [data[i] for i in get_connectivity(element)]
|
||||
end
|
||||
|
||||
function update!{K,V}(element::Element, field_name, data::Pair{Float64, Dict{K, V}})
|
||||
#time, field_data = data
|
||||
#element_data = V[field_data[i] for i in get_connectivity(element)]
|
||||
#update!(element, field_name, time => element_data)
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
element[field_name] = Field(data)
|
||||
end
|
||||
time, field_data = data
|
||||
element_data = V[field_data[i] for i in get_connectivity(element)]
|
||||
update!(element, field_name, time => element_data)
|
||||
#if haskey(element, field_name)
|
||||
# update!(element[field_name], data)
|
||||
#else
|
||||
# element[field_name] = Field(data)
|
||||
#end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name::AbstractString, datas::Union{Real, Vector, Pair{Float64, Union{Float64, Real, Vector{Any}}}}...)
|
||||
|
||||
@@ -19,6 +19,7 @@ type Contact <: BoundaryProblem
|
||||
use_forwarddiff :: Bool
|
||||
minimum_active_set_size :: Int
|
||||
distval :: Float64
|
||||
remove_from_set :: Bool # allow removal of non-potential contact pairs
|
||||
store_fields :: Vector{AbstractString}
|
||||
end
|
||||
|
||||
@@ -26,7 +27,7 @@ function Contact()
|
||||
default_fields = ["element area", "contact area", "weighted gap",
|
||||
"contact pressure", "active nodes", "inactive nodes", "stick nodes",
|
||||
"slip nodes", "complementarity condition", "contact error"]
|
||||
return Contact(-1, false, false, false, true, false, 0, 5.0, default_fields)
|
||||
return Contact(-1, false, false, false, true, false, 0, 5.0, false, default_fields)
|
||||
end
|
||||
|
||||
function get_unknown_field_name(problem::Problem{Contact})
|
||||
|
||||
+64
-23
@@ -21,6 +21,8 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
update!(slave_elements, "normal", time => normals)
|
||||
update!(slave_elements, "tangent", time => tangents)
|
||||
|
||||
Rn = 0.0
|
||||
|
||||
# 2. loop all slave elements
|
||||
for slave_element in slave_elements
|
||||
|
||||
@@ -56,8 +58,13 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
u2 = master_element("displacement", time)
|
||||
x2 = X2 + u2
|
||||
|
||||
norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) < props.distval || continue
|
||||
norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) < props.distval || continue
|
||||
if norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) > props.distval
|
||||
continue
|
||||
end
|
||||
|
||||
if norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) > props.distval
|
||||
continue
|
||||
end
|
||||
|
||||
# 3.1 calculate segmentation
|
||||
xi1a = project_from_master_to_slave(slave_element, X2[1], time)
|
||||
@@ -89,7 +96,11 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
# local mortar matrices
|
||||
fill!(De, 0.0)
|
||||
fill!(Me, 0.0)
|
||||
ge = zeros(field_dim*nsl)
|
||||
Ne = zeros(nsl, 2*nsl)
|
||||
Te = zeros(nsl, 2*nsl)
|
||||
He = zeros(nsl, 2*nsl)
|
||||
ce = zeros(nsl)
|
||||
ge = zeros(nsl)
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ*l
|
||||
@@ -97,40 +108,56 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
|
||||
N1 = vec(get_basis(slave_element, xi_s, time))
|
||||
Phi = Ae*N1
|
||||
|
||||
# project gauss point from slave element to master element in direction n_s
|
||||
X_s = N1*X1 # coordinate in gauss point
|
||||
n_s = N1*n1 # normal direction in gauss point
|
||||
t_s = N1*t1 # tangent condition in gauss point
|
||||
n_s /= norm(n_s)
|
||||
t_s /= norm(t_s)
|
||||
xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
|
||||
N2 = vec(get_basis(master_element, xi_m, time))
|
||||
X_m = N2*X2
|
||||
De += w*Phi*N1'
|
||||
Me += w*Phi*N2'
|
||||
|
||||
u_s = N1*u1
|
||||
u_m = N2*u2
|
||||
x_s = X_s + u_s
|
||||
x_m = X_m + u_m
|
||||
la_s = Phi*la1
|
||||
ge += w*vec((x_m-x_s)*Phi')
|
||||
|
||||
# virtual work
|
||||
De += w*Phi*N1'
|
||||
Me += w*Phi*N2'
|
||||
|
||||
# contact constraints
|
||||
Ne += w*reshape(kron(N1, n_s, Phi), 2, 4)
|
||||
Te += w*reshape(kron(N2, n_s, Phi), 2, 4)
|
||||
He += w*reshape(kron(N1, t_s, Phi), 2, 4)
|
||||
ge += w*Phi*dot(n_s, x_m-x_s)
|
||||
ce += w*N1*dot(n_s, -la_s)
|
||||
Rn += w*dot(n_s, -la_s)
|
||||
|
||||
contact_area += w
|
||||
contact_error += 1/2*w*dot(n_s, x_s-x_m)^2
|
||||
end
|
||||
|
||||
# add contribution to contact virtual work
|
||||
sdofs = get_gdofs(problem, slave_element)
|
||||
mdofs = get_gdofs(problem, master_element)
|
||||
nsldofs = length(sdofs)
|
||||
nmdofs = length(mdofs)
|
||||
D2 = zeros(nsldofs, nsldofs)
|
||||
M2 = zeros(nmdofs, nmdofs)
|
||||
|
||||
# add contribution to contact virtual work
|
||||
for i=1:field_dim
|
||||
D2[i:field_dim:end, i:field_dim:end] += De
|
||||
M2[i:field_dim:end, i:field_dim:end] += Me
|
||||
lsdofs = sdofs[i:field_dim:end]
|
||||
lmdofs = mdofs[i:field_dim:end]
|
||||
add!(problem.assembly.C1, lsdofs, lsdofs, De)
|
||||
add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
|
||||
end
|
||||
add!(problem.assembly.C1, sdofs, sdofs, D2)
|
||||
add!(problem.assembly.C1, sdofs, mdofs, -M2)
|
||||
add!(problem.assembly.C2, sdofs, sdofs, Q2'*D2)
|
||||
add!(problem.assembly.C2, sdofs, mdofs, -Q2'*M2)
|
||||
add!(problem.assembly.g, sdofs, Q2'*ge)
|
||||
|
||||
# add contribution to contact constraints
|
||||
add!(problem.assembly.C2, sdofs[1:field_dim:end], sdofs, Ne)
|
||||
add!(problem.assembly.C2, sdofs[1:field_dim:end], mdofs, -Te)
|
||||
add!(problem.assembly.D, sdofs[2:field_dim:end], sdofs, He)
|
||||
add!(problem.assembly.g, sdofs[1:field_dim:end], ge)
|
||||
add!(problem.assembly.c, sdofs[1:field_dim:end], ce)
|
||||
|
||||
end # master elements done
|
||||
|
||||
@@ -153,13 +180,22 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
is_slip = Dict{Int64, Int}()
|
||||
is_stick = Dict{Int64, Int}()
|
||||
|
||||
g = full(problem.assembly.g)
|
||||
la = problem.assembly.la
|
||||
ndofs = length(la)
|
||||
info("contact ndofs: $ndofs")
|
||||
info("Rn = $Rn")
|
||||
|
||||
C1 = sparse(problem.assembly.C1)
|
||||
C2 = sparse(problem.assembly.C2, ndofs, ndofs)
|
||||
D = sparse(problem.assembly.D, ndofs, ndofs)
|
||||
g = full(problem.assembly.g, ndofs, 1)
|
||||
c = full(problem.assembly.c, ndofs, 1)
|
||||
|
||||
# active / inactive node detection
|
||||
for j in S
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
weighted_gap[j] = g[dofs]
|
||||
|
||||
if length(la) != 0
|
||||
p = dot(normals[j], la[dofs])
|
||||
t = dot(tangents[j], la[dofs])
|
||||
@@ -167,6 +203,8 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
else
|
||||
contact_pressure[j] = [0.0, 0.0]
|
||||
end
|
||||
|
||||
# contact_pressure[j] = c[dofs]
|
||||
complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
|
||||
if complementarity_condition[j][1] < 0
|
||||
is_inactive[j] = 1
|
||||
@@ -178,6 +216,11 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
is_active[j] = 1
|
||||
is_slip[j] = 1
|
||||
is_stick[j] = 0
|
||||
_c1 = complementarity_condition[j][1]
|
||||
_c2 = c[dofs]
|
||||
_c3 = contact_pressure[j][1]
|
||||
_c4 = g[dofs]
|
||||
info("active $j: c1 = $_c1, c2 = $_c2, c3 = $_c3, c4 = $_c4")
|
||||
end
|
||||
end
|
||||
|
||||
@@ -212,12 +255,9 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
|
||||
# solve variational inequality
|
||||
|
||||
C1 = sparse(problem.assembly.C1)
|
||||
ndofs = size(C1, 1)
|
||||
C2 = sparse(problem.assembly.C2)
|
||||
D = spzeros(ndofs, ndofs)
|
||||
|
||||
# constitutive modelling in tangent direction, frictionless contact
|
||||
#=
|
||||
for j in S
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
if (is_active[j] == 1) && (is_slip[j] == 1)
|
||||
@@ -227,6 +267,7 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
D[dofs[2], dofs] = tangents[j]
|
||||
end
|
||||
end
|
||||
=#
|
||||
|
||||
# remove inactive nodes from assembly
|
||||
for j in S
|
||||
|
||||
@@ -68,9 +68,10 @@ function calculate_normals(elements, time, ::Type{Val{1}}; rotate_normals=false)
|
||||
tangents = Dict{Int64, Vector{Float64}}()
|
||||
for element in elements
|
||||
conn = get_connectivity(element)
|
||||
X1 = element("geometry", time)
|
||||
dN = get_dbasis(element, [0.0], time)
|
||||
tangent = vec(sum([kron(dN[:,i], X1[i]') for i=1:length(X1)]))
|
||||
#X1 = element("geometry", time)
|
||||
#dN = get_dbasis(element, [0.0], time)
|
||||
#tangent = vec(sum([kron(dN[:,i], X1[i]') for i=1:length(X1)]))
|
||||
tangent = vec(element([0.0], time, Val{:Jacobian}))
|
||||
for nid in conn
|
||||
if haskey(tangents, nid)
|
||||
tangents[nid] += tangent
|
||||
@@ -116,8 +117,8 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
|
||||
# 1. calculate nodal normals and tangents for slave element nodes j ∈ S
|
||||
normals, tangents = calculate_normals(slave_elements, time, Val{1};
|
||||
rotate_normals=props.rotate_normals)
|
||||
update!(slave_elements, "normal", normals)
|
||||
update!(slave_elements, "tangent", tangents)
|
||||
update!(slave_elements, "normal", time => normals)
|
||||
update!(slave_elements, "tangent", time => tangents)
|
||||
|
||||
# 2. loop all slave elements
|
||||
for slave_element in slave_elements
|
||||
|
||||
@@ -0,0 +1,45 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Preprocess
|
||||
using JuliaFEM.Postprocess
|
||||
using JuliaFEM.Testing
|
||||
|
||||
@testset "test 2d linear elasticity with surface + volume load" begin
|
||||
meshfile = "/geometry/2d_block/BLOCK_1elem.med"
|
||||
mesh = aster_read_mesh(Pkg.dir("JuliaFEM")*meshfile)
|
||||
|
||||
# field problem
|
||||
block = Problem(Elasticity, "BLOCK", 2)
|
||||
block.properties.formulation = :plane_strain
|
||||
block.properties.finite_strain = false
|
||||
block.properties.geometric_stiffness = false
|
||||
block.elements = create_elements(mesh, "BLOCK")
|
||||
update!(block.elements, "youngs modulus", 288.0)
|
||||
update!(block.elements, "poissons ratio", 1/3)
|
||||
|
||||
# traction
|
||||
traction = Problem(Elasticity, "TRACTION", 2)
|
||||
traction.properties.formulation = :plane_strain
|
||||
traction.properties.finite_strain = false
|
||||
traction.properties.geometric_stiffness = false
|
||||
traction.elements = create_elements(mesh, "TOP")
|
||||
update!(traction, "displacement traction force 2", 288.0*9/8)
|
||||
|
||||
# boundary conditions
|
||||
bc_sym_23 = Problem(Dirichlet, "symmetry bc 23", 2, "displacement")
|
||||
bc_sym_23.elements = create_elements(mesh, "LEFT")
|
||||
update!(bc_sym_23, "displacement 1", 0.0)
|
||||
bc_sym_13 = Problem(Dirichlet, "symmetry bc 13", 2, "displacement")
|
||||
bc_sym_13.elements = create_elements(mesh, "BOTTOM")
|
||||
update!(bc_sym_13, "displacement 2", 0.0)
|
||||
|
||||
solver = LinearSolver(block, traction, bc_sym_23, bc_sym_13)
|
||||
solver()
|
||||
|
||||
info("u = ", block.assembly.u)
|
||||
info("λ = ", block.assembly.la)
|
||||
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user