From c3bc708f01ed2bdafbb6c92b2b175e498ac73925 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olli=20V=C3=A4in=C3=B6l=C3=A4?= Date: Fri, 28 Oct 2016 15:46:04 +0300 Subject: [PATCH] Ideal plasticity converged, both 2D and 3D --- src/materials_plasticity.jl | 6 +- src/problems_elasticity.jl | 83 +++++++++++++------ ...sticplastic_3d_linear_with_surface_load.jl | 18 ++-- test/test_von_mises_material.jl | 82 +++++++++++++++--- 4 files changed, 144 insertions(+), 45 deletions(-) diff --git a/src/materials_plasticity.jl b/src/materials_plasticity.jl index d1c000a..42e8c38 100644 --- a/src/materials_plasticity.jl +++ b/src/materials_plasticity.jl @@ -72,7 +72,7 @@ function radial_return(params, dstrain, D, stress_y, stress_base, yield_surface_ [vec(function_1); function_2] end -function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, time, dt, type_) +function ideal_plasticity!(stress_new, stress_last, dstrain_vec, pstrain, D, params, Dtan, yield_surface_, time, dt, type_) # Test stress dstress = vec(D * dstrain_vec) stress_trial = stress_last + dstress @@ -83,6 +83,7 @@ function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan # Calculating and checking for yield yield = yield_curr(stress_trial) if isless(yield, 0.0) + stress_new[:] = stress_trial[:] Dtan[:,:] = D[:,:] else @@ -103,6 +104,7 @@ function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan # Updating stress stress_new[:] = stress_last + dstress + # Calculating plastic strain dfds_ = x -> ForwardDiff.gradient(yield_curr, x) dep = plastic_multiplier * dfds_(vec(stress_new)) @@ -114,6 +116,6 @@ function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan Dc = (D^-1 + plastic_multiplier * D2g(stress_new))^-1 dfds = dfds_(stress_new) Dtan[:,:] = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1] - + pstrain[:] = plastic_multiplier * dfds end end diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 1a12c5a..83ae94b 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -71,7 +71,7 @@ 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}}) +function initialize_internal_params!(params, ip, type_) #::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))) @@ -80,8 +80,15 @@ function initialize_internal_params!(params, ip, ::Type{Val{:type_2d}}) 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]) + if type_ == Val{:type_2d} + update!(ip, "stress", 0.0 => [0.0,0.0,0.0]) + update!(ip, "strain", 0.0 => [0.0,0.0,0.0]) + elseif type_ == Val{:type_3d} + update!(ip, "stress", 0.0 => [0.0,0.0,0.0,0.0,0.0,0.0]) + update!(ip, "strain", 0.0 => [0.0,0.0,0.0,0.0,0.0,0.0]) + else + error("daa") + end update!(ip, "prev_time", 0.0 => 0.0) update!(ip, "params_initialized", 0.0 => true) end @@ -93,13 +100,13 @@ function get_keys(element) 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 +#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}}) @@ -513,21 +520,48 @@ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, 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}) + 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_3d}) + + 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, 0.0, 0.0, 0.0] + plastic_strain = [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, plastic_strain, D, params, Dtan, yield_surface_, time, dt, Val{:type_3d}) + + + # 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 + + # 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 @@ -540,6 +574,7 @@ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[4]) :stress23 in props.store_fields && update!(ip, "stress23", time => stress_vec[5]) :stress13 in props.store_fields && update!(ip, "stress13", time => stress_vec[6]) + :plastic_strain in props.store_fields && update!(ip, "plastic_strain", time => plastic_strain) Km += w*BL'*Dtan*BL # material stiffness end diff --git a/test/test_elasticplastic_3d_linear_with_surface_load.jl b/test/test_elasticplastic_3d_linear_with_surface_load.jl index 38f3102..66f77ea 100644 --- a/test/test_elasticplastic_3d_linear_with_surface_load.jl +++ b/test/test_elasticplastic_3d_linear_with_surface_load.jl @@ -16,22 +16,23 @@ using JuliaFEM.Testing 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) + element = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8]) + update!([element], "geometry", nodes) + update!([element], "youngs modulus", 200e3) + 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)) + "params" => Dict("yield_stress" => 400.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) + push!(elasticity_problem.properties.store_fields, :plastic_strain) + push!(elasticity_problem, element) bc = Element(Quad4, [1,4,8,5]) update!([bc], "geometry", nodes) @@ -48,12 +49,13 @@ using JuliaFEM.Testing push!(boundary_motion, disp) solver = NonlinearSolver("solve block problem") + solver.time = 1.0 push!(solver, elasticity_problem) push!(solver, boundary_problem) push!(solver, boundary_motion) solver() - disp = element1("displacement", [1.0, 1.0, 1.0], 0.0) + disp = element("displacement", [1.0, 1.0, 1.0], 1.0) info("displacement at tip: $disp") u_expected = 2.0 * [-1/3, -1/3, 1.0] # @test isapprox(disp, u_expected) diff --git a/test/test_von_mises_material.jl b/test/test_von_mises_material.jl index 313fe79..f397076 100644 --- a/test/test_von_mises_material.jl +++ b/test/test_von_mises_material.jl @@ -64,17 +64,17 @@ function test_von_mises_3D_basic() 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 - 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_new) - eig_vals[i, :] = sort(eigvals(eig_stress)) - stress_last[:] = stress_new[:] - end + #for i=1:steps + # strain_new = reshape(strain_tot[i, :, :], (6, 1)) + # dstrain = strain_new - strain + # 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_new) + # eig_vals[i, :] = sort(eigvals(eig_stress)) + # stress_last[:] = stress_new[:] + #end toc() # ================ Plotting =================== # @@ -118,6 +118,66 @@ function test_von_mises_3D_basic() info("Calculation finished") # plot3D(ee, ss) + + # plot the surface + xx = zeros(10, 10) + yy = zeros(10, 10) + + for i=1:10 + for j=1:10 + xx[i, j] = (i - 5) * 100 + yy[i, j] = (j - 5) * 100 + end + end + + # calculate corresponding z + z = zeros(10, 10) + for i=1:10 + for j=1:10 + z[i, j] = 1 + end + end + + # ================================================================== + # plot the surface + plot_surface(xx, yy, z, color="blue") + + stress_y = 200.0 + + 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) - stress_y + f(m,c) = [600*cos(c) 600*sin(c)].*m + x_vals = [] + max_iter = 100 + y_vals = [] + for i=0:0.1:(2*pi+0.3) + wf(x) = f(x, i) + t = 0.01 + step = 2 + merkki = -1 + s11, s22 = wf(t) + ii = 0 + while (abs(vm(s11, s22)) > 1e-7) && ii < max_iter + val = vm(s11, s22) + if sign(val) != merkki + merkki *= -1 + step *= -0.5 + end + t += step + s11, s22 = wf(t) + ii += 1 + end + push!(x_vals, s11) + push!(y_vals, s22) + end + plot(x_vals, y_vals, zeros(length(y_vals)), color="yellow") + axis("equal") + # ================================================================== + + 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")