diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 3240258..f5a25f0 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -67,6 +67,9 @@ export Problem, AbstractProblem, FieldProblem, BoundaryProblem, include("problems_elasticity.jl") export Elasticity +include("materials_plasticity.jl") +export plastic_von_mises + include("problems_dirichlet.jl") export Dirichlet diff --git a/src/elements.jl b/src/elements.jl index 59b9c0b..5925559 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -71,7 +71,7 @@ julia> el([0.0, 0.0], 0.0, 1) julia> el([0.0, 0.0], 0.0, 2) 2x8 Array{Float64,2}: - 0.25 0.0 0.25 0.0 0.25 0.0 0.25 0.0 + 0.25 0.0 0.25 0.0 0.25 0.0 0.25 0.0 0.0 0.25 0.0 0.25 0.0 0.25 0.0 0.25 """ diff --git a/src/materials_plasticity.jl b/src/materials_plasticity.jl new file mode 100644 index 0000000..d1c000a --- /dev/null +++ b/src/materials_plasticity.jl @@ -0,0 +1,119 @@ +using ForwardDiff +using NLsolve + +""" +Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values +""" +function find_root!(f, df, x; max_iter=50, norm_acc=1e-9) + converged = false + for i=1:max_iter + dx = -df(x) \ f(x) + x += dx + norm(dx) < norm_acc && (converged = true; break) + end + converged || error("No convergence in radial return!") + return x +end + +""" +Equivalent tensile stress. + +More info can be found from: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion + Section: Reduced von Mises equation for different stress conditions +""" +function equivalent_stress(stress, ::Type{Val{:type_3d}}) + stress_ten = [stress[1] stress[6] stress[5]; + stress[6] stress[2] stress[4]; + stress[5] stress[4] stress[3]] + stress_dev = stress_ten - 1/3 * trace(stress_ten) * eye(3) + s = vec(stress_dev) + return sqrt(3/2 * dot(s, s)) +end + +""" +http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm + +von mises: plane stress +https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf +""" +function equivalent_stress(stress, ::Type{Val{:type_2d}}) + s1, s2, t12 = stress + # Calculating principal stresses + # http://www.engineersedge.com/material_science/principal_vonmises_stress__13418.htm + se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2) + se2 = (s1 + s2)/2 - sqrt(((s1 - s2)/2)^2 + t12^2) + + return sqrt(se1^2 -se1*se2 + se2^2) +end + +""" +https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf +""" +function yield_function(stress, stress_y, ::Type{Val{:von_mises}}, type_) + equivalent_stress(stress, type_) - stress_y +end + +function radial_return(params, dstrain, D, stress_y, stress_base, yield_surface_, type_) + + # Creating wrapper for gradient + vm_wrap(stress_) = yield_function(stress_, stress_y, yield_surface_, type_) + dfds = x -> ForwardDiff.gradient(vm_wrap, x) + + # Stress rate and total strain + dstress = params[1:end-1] + stress_tot = stress_base + dstress + + # Calculating plastic strain rate + dstrain_p = params[end] * dfds(stress_tot) + + # Calculating equations + function_1 = dstress - D * (dstrain - dstrain_p) + function_2 = vm_wrap(stress_tot) + [vec(function_1); function_2] +end + +function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, time, dt, type_) + # Test stress + dstress = vec(D * dstrain_vec) + stress_trial = stress_last + dstress + stress_y = params["yield_stress"] + + yield_curr = x -> yield_function(x, stress_y, yield_surface_, type_) + + # Calculating and checking for yield + yield = yield_curr(stress_trial) + if isless(yield, 0.0) + stress_new[:] = stress_trial[:] + Dtan[:,:] = D[:,:] + else + # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ \ f and initial values + f = stress_ -> radial_return(stress_, dstrain_vec, D, stress_y, stress_last, yield_surface_, type_) + df = x -> ForwardDiff.jacobian(f, x) + + # Calculating root (two options) + vals = [vec(stress_trial - stress_last); 0.0] + + #results = nlsolve(not_in_place(f), vals).zero + results = find_root!(f, df, vals) + + # extracting results + dstress = results[1:end-1] + plastic_multiplier = results[end] + + # Updating stress + stress_new[:] = stress_last + dstress + + # Calculating plastic strain + dfds_ = x -> ForwardDiff.gradient(yield_curr, x) + dep = plastic_multiplier * dfds_(vec(stress_new)) + + # Equations for consistent tangent matrix can be found from: + # http://homes.civil.aau.dk/lda/continuum/plast.pdf + # equations: 152 & 153 + D2g = x -> ForwardDiff.hessian(yield_curr, x) + Dc = (D^-1 + plastic_multiplier * D2g(stress_new))^-1 + dfds = dfds_(stress_new) + Dtan[:,:] = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1] + + end +end diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index b7dc871..1a12c5a 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -71,6 +71,35 @@ typealias Elasticity2DVolumeElements Union{Tri3, Tri6, Quad4, Quad8, Quad9} typealias Elasticity3DSurfaceElements Union{Poi1, Tri3, Tri6, Quad4, Quad8, Quad9} typealias Elasticity3DVolumeElements Union{Tet4, Wedge6, Hex8, Tet10, Hex20, Hex27} +function initialize_internal_params!(params, ip, ::Type{Val{:type_2d}}) + param_keys = keys(params) + all_keys = ip.fields.keys + ip_fields = filter(x->isdefined(all_keys, x), collect(1:length(all_keys))) + + if !("params_initialized" in ip_fields) + for key in param_keys + update!(ip, key, 0.0 => params[key]) + end + update!(ip, "stress", 0.0 => [0.0,0.0,0.0]) + update!(ip, "strain", 0.0 => [0.0,0.0,0.0]) + update!(ip, "prev_time", 0.0 => 0.0) + update!(ip, "params_initialized", 0.0 => true) + end +end + +function get_keys(element) + all_keys = element.fields.keys + idx = filter(x->isdefined(all_keys, x), collect(1:length(all_keys))) + map(x -> all_keys[x], idx) +end + +function initialize_internal_params!(params, ip_id, ::Type{Val{:type_3d}}) + if !(ip_id in keys(params)) + params[ip_id] = Dict{Any, Any}() + params[ip_id]["last_stress"] = [0.0,0.0,0.0,0.0,0.0,0.0] + params[ip_id]["last_strain"] = [0.0,0.0,0.0,0.0,0.0,0.0] + end +end """ Elasticity equations for 2d cases. """ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}}) @@ -83,6 +112,7 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, Km = zeros(dim*nnodes, dim*nnodes) Kg = zeros(dim*nnodes, dim*nnodes) f = zeros(dim*nnodes) + Dtan = zeros(3,3) for ip in get_integration_points(element) @@ -90,7 +120,6 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, w = ip.weight*detJ N = element(ip, time) dN = element(ip, time, Val{:Grad}) - # kinematics gradu = element("displacement", ip, time, Val{:Grad}) @@ -129,15 +158,45 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, nu 1.0 0.0 0.0 0.0 (1.0-nu)/2.0] elseif props.formulation == :plane_strain - D = E/((1+nu)*(1-2*nu)) .* [ - 1-nu nu 0 - nu 1-nu 0 - 0 0 (1-2*nu)/2] + D = E/((1.0+nu)*(1.0-2.0*nu)) .* [ + 1.0-nu nu 0.0 + nu 1.0-nu 0.0 + 0.0 0.0 (1.0-2.0*nu)/2.0] else error("unknown plane formulation: $(props.formulation)") end + # calculate stress - stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec) + element_keys = get_keys(element) + + if "plasticity" in element_keys + plastic_def = element("plasticity")[ip.id] + + calculate_stress! = plastic_def["type"] + yield_surface_ = plastic_def["yield_surface"] + params = plastic_def["params"] + + initialize_internal_params!(params, ip, Val{:type_2d}) + + if time == 0.0 + error("Given step time = $(time). Please select time > 0.0") + end + + t_last = ip("prev_time", time) + update!(ip, "prev_time", time => t_last) + + dt = time - t_last + + stress_last = ip("stress", t_last) + strain_last = ip("strain", t_last) + + dstrain_vec = strain_vec - strain_last + stress_vec = [0.0, 0.0, 0.0] + calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, time, dt, Val{:type_2d}) + else + stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec) + Dtan[:,:] = D[:,:] + end :strain in props.store_fields && update!(ip, "strain", time => strain_vec) :stress in props.store_fields && update!(ip, "stress", time => stress_vec) @@ -145,7 +204,8 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, :stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2]) :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[3]) - Km += w*BL'*D*BL + Km += w*BL'*Dtan*BL + # stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]] # cauchy_stress = F'*stress*F/det(F) @@ -383,7 +443,6 @@ end """ Elasticity equations, 3d nonlinear. """ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}}) - props = problem.properties dim = get_unknown_field_dimension(problem) nnodes = length(element) @@ -450,7 +509,28 @@ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, 0.0 0.0 0.0 0.5-nu 0.0 0.0 0.0 0.0 0.0 0.0 0.5-nu 0.0 0.0 0.0 0.0 0.0 0.0 0.5-nu] - stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec) + + element_keys = get_keys(element) + + if "plasticity" in element_keys + plastic_def = element.dev["plasticity"] + calculate_stress! = plastic_def["stress"] + params = plastic_def["params"] + yield_surface_ = plastic_def["yield_surface"] + (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:type_3d}) + dstrain_vec = strain_vec - strain_last + stress_vec = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + Dtan = [0.0 0.0 0.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 0.0] + calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, Val{:type_3d}) + else + stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec) + Dtan = D + end :strain in props.store_fields && update!(ip, "strain", time => strain_vec) :stress in props.store_fields && update!(ip, "stress", time => stress_vec) @@ -461,8 +541,7 @@ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, :stress23 in props.store_fields && update!(ip, "stress23", time => stress_vec[5]) :stress13 in props.store_fields && update!(ip, "stress13", time => stress_vec[6]) - Km += w*BL'*D*BL - + Km += w*BL'*Dtan*BL # material stiffness end if props.geometric_stiffness diff --git a/src/vonmises.jl b/src/vonmises.jl deleted file mode 100644 index ecc477f..0000000 --- a/src/vonmises.jl +++ /dev/null @@ -1,304 +0,0 @@ -using ForwardDiff - -""" -Create a isotropic Hooke material matrix C - -More information: http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_isotropic.cfm - https://en.wikipedia.org/wiki/Hooke's_law - http://www.ce.berkeley.edu/~sanjay/ce231mse211/symidentity.pdf -Parameters ----------- - E: Float - Elastic modulus - ν: Float - Poisson constant - -Returns -------- - Array{Float64, (6,6)} -""" -function stiffnessTensor(E, ν) - a = 1 - ν - b = 1 - 2*ν - c = 1 + ν - multiplier = E / (b * c) - return Float64[a ν ν 0 0 0; - ν a ν 0 0 0; - ν ν a 0 0 0; - 0 0 0 b 0 0; - 0 0 0 0 b 0; - 0 0 0 0 0 b].*multiplier -end - -# Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values -function find_root!(f, df, x; max_iter=50, norm_acc=1e-10) - converged = false - for i=1:max_iter - dx = df(x) \ -f(x) - x += dx - norm(dx) < norm_acc && (converged = true; break) - end - converged || error("no convergence!") - x -end - -type State - C :: Array{Float64, 2} - stress_y :: Float64 - stress :: Array{Float64, 1} - strain :: Array{Float64, 1} -end - -""" -Equivalent tensile stress. - -More info can be found from: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion - Section: Reduced von Mises equation for different stress conditions - -Parameters ----------- - σ: Array{Float64, 6} - Stress in Voigt notation - -Returns -------- - Float -""" -function stress_eq(stress) - stress_ten = [stress[1] stress[6] stress[5]; - stress[6] stress[2] stress[4]; - stress[5] stress[4] stress[3]] - stress_dev = stress_ten - 1/3 * trace(stress_ten) * eye(3) - s = vec(stress_dev) - return sqrt(3/2 * dot(s, s)) -end - - -""" -Von Mises Yield criterion - -More info can be found from: http://csm.mech.utah.edu/content/wp-content/uploads/2011/10/9tutorialOnJ2Plasticity.pdf - -Parameters ----------- - σ: Array{Float64, 6} - Stress in Voigt notation - k: Float64 - Material constant, Yield limit - -Returns -------- - Float -""" -function vonMisesYield(stress, stress_y) - stress_eq(stress) - stress_y -end - -""" -Function for NLsolve. Inside this function are the equations which we want to find root. -Ψ is the yield function below. Functions defined here: - - dσ - C (dϵ - dλ*dΨ/dσ) = 0 - σₑ(σ) - k = 0 - -Parameters ----------- - params: Array{Float64, 7} - Array containing values from solver - dϵ: Array{Float64, 6} - Strain rate vector in Voigt notation - C: Array{Float64, (6, 6)} - Material tensor - k: Float - Material constant, yield limit - Δt: Float - time increment - σ_begin:Array{Float64, 6} - Stress vector in Voigt notation - -Returns -------- - Array{Float64, 7}, return values for solver -""" -function vonMisesRoot(params, dstrain, C, stress_y, stress_base) - - # Creating wrapper for gradient - vm_wrap(stress_) = vonMisesYield(stress_, stress_y) - dfds = ForwardDiff.gradient(vm_wrap) - - # Stress rate and total strain - dstress = params[1:6] - stress_tot = vec(stress_base) + params[1:6] - - # Calculating plastic strain rate - dstrain_p = params[end] * dfds(stress_tot) - - # Calculating equations - function_1 = dstress - C * (dstrain - dstrain_p) - function_2 = vm_wrap(stress_tot) - [vec(function_1); function_2] -end - - - -""" -Stress for ideal plastic von Mises material model - -Parameters ----------- - dϵ: Array{Float64, 6} - Strain rate vector in Voigt notation - Δt: Float - time increment - σ: Array{Float64, 6} - Last stress vector in Voigt notation - C: Array{Float64, (6, 6)} - Material tensor - k: Float - Material constant, yield limit - -Returns -------- - Tuple - Plastic strain rate dϵᵖ and new stress vector σ -""" -function calculate_stress!(dstrain, mat::State, ::Type{Val{:vonMises}}) - stress = mat.stress - C = mat.C - stress_y = mat.stress_y - # Test stress - stress_tria = stress + C * dstrain - - # Calculating and checking for yield - yield = vonMisesYield(stress_tria, stress_y) - if isless(yield, 0.0) - mat.stress = vec(stress_tria) - else - # Yielding happened - # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values - initial_guess = Float64[vec(stress_tria - stress); 0.1] - f(stress_) = vonMisesRoot(stress_, dstrain, C, stress_y, stress) - df = ForwardDiff.jacobian(f) - - # Calculating root - result = nlsolve(not_in_place(f, df), initial_guess).zero - mat.stress += result[1:6] - end -end - -function calculate_stress(dstrain, stress, C, stress_y, - ::Type{Val{:vonMises}}, - ::Type{Val{:ElasticPlasticProblem}}) - # Test stress - stress_tria = stress + C * dstrain - - # Calculating and checking for yield - yield = vonMisesYield(stress_tria, stress_y) - if isless(yield, 0.0) - # stress[i] = stress_tria[i] - return 0.0 - else - # Yielding happened - # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values - x = [vec(stress_tria - stress); 0.0] - f(stress_) = vonMisesRoot(stress_, dstrain, C, stress_y, stress) - df = ForwardDiff.jacobian(f) - - # Calculating root - # result = nlsolve(not_in_place(f, df), initial_guess).zero - max_iter = 10 - converged = false - for i=1:5 - dx = df(x) \ -f(x) - x += dx - # println(x) - norm(dx) < 1e-10 && (converged = true; break) - end - converged || error("no convergence!") - # stress[:] += x[1:6] - return x[end] - end -end - -################################################################################## -# ----- AFTER THIS POINT: VON MISES : PLANE STRESS IMPLEMENTATION ----- # -################################################################################## - -""" -http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm -""" -function stiffnessTensorPlaneStress(E, ν) - a = 1 - ν^2 - b = 1 - ν - multiplier = E / a - return Float64[1 ν 0; - ν 1 0; - 0 0 b].*multiplier -end - -# von mises: plane stress -# https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf -function stress_eq_plane_stress(stress) - s1, s2, t12 = stress - # Calculating principal stresses - # http://www.engineersedge.com/material_science/principal_vonmises_stress__13418.htm - se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2) - se2 = (s1 + s2)/2 - sqrt(((s1 - s2)/2)^2 + t12^2) - return sqrt(se1^2 -se1*se2 + se2^2) -end - -# https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf -function vonMisesYieldPlaneStress(stress, stress_y) - stress_eq_plane_stress(stress) - stress_y -end - -function vonMisesRootPlaneStress(params, dstrain, C, stress_y, stress_base) - - # Creating wrapper for gradient - vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y) - dfds = ForwardDiff.gradient(vm_wrap) - - # Stress rate and total strain - dstress = params[1:3] - stress_tot = vec(stress_base) + params[1:3] - - # Calculating plastic strain rate - dstrain_p = params[end] * dfds(stress_tot) - - # Calculating equations - function_1 = dstress - C * (dstrain - dstrain_p) - function_2 = vm_wrap(stress_tot) - [vec(function_1); function_2] -end - -function calculate_stress(dstrain, stress, C, stress_y, - ::Type{Val{:vonMises}}, - ::Type{Val{:PlaneStressElasticPlasticProblem}}) - # Test stress - dstress = C * dstrain - stress_tria = stress + dstress - - # Calculating and checking for yield - yield = vonMisesYieldPlaneStress(stress_tria, stress_y) - if isless(yield, 0.0) - return dstress, zeros(3) - else - info("yielded") - # Yielding happened - # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values - x = [vec(stress_tria - stress); 0.0] - f(stress_) = vonMisesRootPlaneStress(stress_, dstrain, C, stress_y, stress) - df = ForwardDiff.jacobian(f) - # Calculating root - results = find_root!(f, df, x) - dstress = results[1:3] - stress_tot = stress + dstress - plastic_multiplier = results[end] - vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y) - dfds = ForwardDiff.gradient(vm_wrap) - dep = plastic_multiplier * dfds(vec(stress_tot)) - info("II ", stress_tot) - info(vm_wrap(stress_tot)) - return dstress, dep - end -end diff --git a/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl b/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl index d4d3754..268fe87 100644 --- a/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl +++ b/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl @@ -70,4 +70,3 @@ using JuliaFEM.Testing end =# end - diff --git a/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl new file mode 100644 index 0000000..5eee961 --- /dev/null +++ b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl @@ -0,0 +1,59 @@ +# 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.Testing + +# @testset "2d nonlinear elasticity: test nonhomogeneous boundary conditions and stress calculation" begin + + # field problem + block = Problem(Elasticity, "BLOCK", 2) + block.properties.formulation = :plane_stress + block.properties.finite_strain = true + block.properties.geometric_stiffness = true + + nodes = Dict{Int, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [1.0, 0.0], + 3 => [1.0, 1.0], + 4 => [0.0, 1.0]) + + element = Element(Quad4, [1, 2, 3, 4]) + update!(element, "geometry", nodes) + update!(element, "youngs modulus", 288.0) + update!(element, "poissons ratio", 1/3) + + plastic_parameters = Dict{Any, Any}("type" => JuliaFEM.ideal_plasticity!, + "yield_surface" => Val{:von_mises}, + "params" => Dict("yield_stress" => 175.0)) + to_integ_points = Dict() + map(x-> to_integ_points[x] = plastic_parameters, get_connectivity(element)) + update!(element, "plasticity", to_integ_points) + push!(block, element) + + # boundary conditions + bc = Problem(Dirichlet, "bc", 2, "displacement") + bel1 = Element(Seg2, [1, 2]) + bel2 = Element(Seg2, [3, 4]) + bel3 = Element(Seg2, [4, 1]) + update!([bel1, bel2, bel3], "geometry", nodes) + update!(bel1, "displacement 2", 0.0) + update!(bel2, "displacement 2", 0.5) + update!(bel3, "displacement 1", 0.0) + push!(bc, bel1, bel2, bel3) + + solver = NonlinearSolver("solve block problem") + solver.time = 1.0 + push!(solver, block, bc) + solver() + + # from code aster + eps_expected = [-2.08333312468287E-01, 6.25000000000000E-01, 0.0] + sig_expected = [ 4.50685020821470E-06, 4.62857140373777E+02, 0.0] + u3_expected = [-2.36237356855269E-01, 5.00000000000000E-01] + + u3 = reshape(block.assembly.u, 2, 4)[:, 3] + info("u3 = $u3") + #@test isapprox(u3, u3_expected, atol=1.0e-5) +# end diff --git a/test/test_elasticplastic_3d_linear_with_surface_load.jl b/test/test_elasticplastic_3d_linear_with_surface_load.jl new file mode 100644 index 0000000..38f3102 --- /dev/null +++ b/test/test_elasticplastic_3d_linear_with_surface_load.jl @@ -0,0 +1,95 @@ +# 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.Testing + +#@testset "test continuum 3d linear elasticity with surface load" begin + nodes = Dict{Int64, Node}( + 1 => [0.0, 0.0, 0.0], + 2 => [1.0, 0.0, 0.0], + 3 => [1.0, 1.0, 0.0], + 4 => [0.0, 1.0, 0.0], + 5 => [0.0, 0.0, 1.0], + 6 => [1.0, 0.0, 1.0], + 7 => [1.0, 1.0, 1.0], + 8 => [0.0, 1.0, 1.0]) + + element1 = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8]) + update!([element1], "geometry", nodes) + update!([element1], "youngs modulus", 200e3) + update!([element1], "poissons ratio", 0.3) + + plastic_parameters = Dict{Any, Any}("type" => JuliaFEM.ideal_plasticity!, + "yield_surface" => Val{:von_mises}, + "params" => Dict("yield_stress" => 175.0)) + to_integ_points = Dict() + map(x-> to_integ_points[x] = plastic_parameters, get_connectivity(element)) + update!(element, "plasticity", to_integ_points) + + elasticity_problem = Problem(Elasticity, "solve continuum block", 3) + elasticity_problem.properties.finite_strain = false + elasticity_problem.properties.geometric_stiffness = false + push!(elasticity_problem, element1) + + bc = Element(Quad4, [1,4,8,5]) + update!([bc], "geometry", nodes) + bc["displacement 1"] = 0.0 + bc["displacement 2"] = 0.0 + bc["displacement 3"] = 0.0 + boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement") + push!(boundary_problem, bc) + + disp = Element(Quad4, [2,3,7,6]) + update!([disp], "geometry", nodes) + disp["displacement 1"] = 0.002 + boundary_motion = Problem(Dirichlet, "displacement bc", 3, "displacement") + push!(boundary_motion, disp) + + solver = NonlinearSolver("solve block problem") + push!(solver, elasticity_problem) + push!(solver, boundary_problem) + push!(solver, boundary_motion) + solver() + + disp = element1("displacement", [1.0, 1.0, 1.0], 0.0) + info("displacement at tip: $disp") + u_expected = 2.0 * [-1/3, -1/3, 1.0] + # @test isapprox(disp, u_expected) +#end + +# function solve_rod_model_elasticity(eltype) +# fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med" +# mesh = aster_read_mesh(fn, eltype) +# element_sets = join(keys(mesh.element_sets), ", ") +# info("element sets: $element_sets") +# p1 = Problem(Elasticity, "rod", 3) +# p2 = Problem(Elasticity, "trac", 3) +# p3 = Problem(Dirichlet, "fixed", 3, "displacement") +# p4 = Problem(Dirichlet, "fixed", 3, "displacement") +# p5 = Problem(Dirichlet, "fixed", 3, "displacement") +# p1.elements = create_elements(mesh, "ROD") +# p2.elements = create_elements(mesh, "FACE2") +# p3.elements = create_elements(mesh, "FACE1") +# p4.elements = create_elements(mesh, "FACE3") +# p5.elements = create_elements(mesh, "FACE5") +# update!(p1, "youngs modulus", 96.0) +# update!(p1, "poissons ratio", 1/3) +# update!(p2, "displacement traction force 1", 96.0) +# update!(p3, "displacement 1", 0.0) +# update!(p4, "displacement 2", 0.0) +# update!(p5, "displacement 3", 0.0) +# solver = LinearSolver(p1, p2, p3, p4, p5) +# solver() +# u_max = maximum(p1.assembly.u) +# info("$eltype, u_max = $u_max") +# return u_max +# end +# @testset "compare 3d rod to CA solution" begin +# @test isapprox(solve_rod_model_elasticity("Tet4"), 0.2) +# @test isapprox(solve_rod_model_elasticity("Tet10"), 0.2) +# @test isapprox(solve_rod_model_elasticity("Hex8"), 0.2) +# @test isapprox(solve_rod_model_elasticity("Hex20"), 0.2) +# @test isapprox(solve_rod_model_elasticity("Hex27"), 0.2) +# end diff --git a/test/test_von_mises_material.jl b/test/test_von_mises_material.jl index 27012b9..313fe79 100644 --- a/test/test_von_mises_material.jl +++ b/test/test_von_mises_material.jl @@ -1,5 +1,6 @@ -#using PyPlot +using PyPlot +using JuliaFEM using JuliaFEM.Testing #using JuliaFEM.MaterialModels: stiffnessTensor, calculate_stress, State #using JuliaFEM.MaterialModels: stiffnessTensorPlaneStress @@ -13,7 +14,14 @@ function test_von_mises_3D_basic() E = 200.0e3 nu = 0.3 ν = 0.3 - C = stiffnessTensor(E, ν) + nu = 0.3 + C = E/((1.0+nu)*(1.0-2.0*nu)) * [ + 1.0-nu nu nu 0.0 0.0 0.0 + nu 1.0-nu nu 0.0 0.0 0.0 + nu nu 1.0-nu 0.0 0.0 0.0 + 0.0 0.0 0.0 0.5-nu 0.0 0.0 + 0.0 0.0 0.0 0.0 0.5-nu 0.0 + 0.0 0.0 0.0 0.0 0.0 0.5-nu] strain_tot = zeros(Float64, (steps, 6)) strain_tot2 = zeros(Float64, (steps, 6)) @@ -49,34 +57,23 @@ function test_von_mises_3D_basic() a[3, 2] = b[4] end - mat = State(C, stress_y, zeros(Float64, 6), zeros(Float64, 6)) - info("Starting calculation") tic() - #= - for i=1:steps - strain_new = reshape(strain_tot[i, :, :], (6, 1)) - dstrain = strain_new - mat.strain - calculate_stress!(dstrain, mat, Val{:vonMises}) - mat.strain += vec(dstrain) - push!(ss, mat.stress[1]) - push!(ee, mat.strain[1]) - - fill_tensor(eig_stress, mat.stress) - eig_vals[i, :] = sort(eigvals(eig_stress)) - end - =# - stress = zeros(Float64, 6) + params = Dict("yield_stress" => stress_y) + stress_new = zeros(Float64, 6) + stress_last = zeros(Float64, 6) strain = zeros(Float64, 6) + Dtan = zeros(6,6) for i=1:steps strain_new = reshape(strain_tot[i, :, :], (6, 1)) dstrain = strain_new - strain - calculate_stress!(dstrain, stress, C, stress_y, Val{:vonMises}) - strain = vec(strain_new) + JuliaFEM.plastic_von_mises!(stress_new, stress_last, dstrain, C, params, Dtan, Val{:type_3d}) + strain[:] = vec(strain_new)[:] push!(ss, stress[1]) push!(ee, strain[1]) - fill_tensor(eig_stress, stress) + fill_tensor(eig_stress, stress_new) eig_vals[i, :] = sort(eigvals(eig_stress)) + stress_last[:] = stress_new[:] end toc() @@ -119,8 +116,8 @@ function test_von_mises_3D_basic() info("Calculation finished") - #PyPlot.plot(ee, ss) - #= + # plot3D(ee, ss) + plot3D(eig_vals[:, 1], eig_vals[:, 2], eig_vals[:, 3], color="red") PyPlot.title("Stress path and von Mises yield surface") PyPlot.xlabel("Eig Stress 1") @@ -128,18 +125,21 @@ function test_von_mises_3D_basic() PyPlot.zlabel("Eig Stress 3") PyPlot.grid() PyPlot.show() - =# + end function test_von_mises_planestress_basic() steps = 1000 - strain_max = 0.003 - num_cycles = 5 - E = 200.0e3 - nu = 0.3 + strain_max = 0.004 + num_cycles = 1. + E = 200000. + nu = 0.3 ν = 0.3 - C = stiffnessTensorPlaneStress(E, ν) + C = E/((1+nu)*(1-2*nu)) .* [ + 1-nu nu 0 + nu 1-nu 0 + 0 0 (1-2*nu)/2] strain_tot = zeros(Float64, (steps, 3)) @@ -151,7 +151,7 @@ function test_von_mises_planestress_basic() strain_last = zeros(Float64, (3)) strain_p = zeros(Float64, (3)) stress = zeros(Float64, (3, 1)) - stress_y = 200.0 + stress_y = 400 ss = Float64[] ee = Float64[] @@ -161,50 +161,39 @@ function test_von_mises_planestress_basic() eig_stress = zeros(Float64, (3, 3)) eig_vals = zeros(Float64, (steps, 3)) - #mat = State(C, stress_y, zeros(Float64, 6), zeros(Float64, 6)) info("Starting calculation") tic() - #= - for i=1:steps - strain_new = reshape(strain_tot[i, :, :], (6, 1)) - dstrain = strain_new - mat.strain - calculate_stress!(dstrain, mat, Val{:vonMises}) - mat.strain += vec(dstrain) - push!(ss, mat.stress[1]) - push!(ee, mat.strain[1]) - fill_tensor(eig_stress, mat.stress) - eig_vals[i, :] = sort(eigvals(eig_stress)) - end - =# - stress = zeros(Float64, 3) + stress_new = zeros(Float64, 3) + stress_last = zeros(Float64, 3) strain = zeros(Float64, 3) + strain_last = zeros(Float64, 3) + params = Dict("yield_stress" => stress_y) + #Dtan = C + Dtan = zeros(3,3) for i=1:steps - strain_new = reshape(strain_tot[i, :, :], (3, 1)) + println("last stress: ", round(stress_last, 2)) + strain_new = vec(strain_tot[i, :, :]) dstrain = strain_new - strain - stress_inc, lambda = calculate_stress(dstrain, - stress, - C, - stress_y, - Val{:vonMises}, - Val{:PlaneStressElasticPlasticProblem}) - stress += stress_inc - strain = vec(strain_new) - s1, s2, t12 = stress + println("analytical stress: ", round((C * strain_new)', 2)) + + JuliaFEM.plastic_von_mises!(stress_new, stress_last, dstrain, C, params, Dtan, Val{:type_2d}) + strain[:] = vec(strain_new)[:] + s1, s2, t12 = stress_new se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2) se2 = (s1 + s2)/2 - sqrt(((s1 - s2)/2)^2 + t12^2) push!(ss, se1) push!(ee, se2) + stress_last[:] = stress_new[:] end - toc() function vm_upper(a, c) vals = f(a[1], a[2], c) vm(vals[1], vals[2], 200) end - vm(a,b) = sqrt(a^2 - a*b + b^2) - 200 + vm(a,b) = sqrt(a^2 - a*b + b^2) - stress_y f(m,c) = [600*cos(c) 600*sin(c)].*m x_vals = [] max_iter = 100 @@ -229,15 +218,12 @@ function test_von_mises_planestress_basic() push!(x_vals, s11) push!(y_vals, s22) end - #= - PyPlot.plot(x_vals, y_vals) - PyPlot.plot(ee, ss) - PyPlot.grid() - PyPlot.show() - =# + + plot(x_vals, y_vals) + plot(ee, ss) + show() end -# test_von_mises_3D_basic() - -#test_von_mises_planestress_basic() +test_von_mises_3D_basic() +# test_von_mises_planestress_basic()