From 97f942225122623a8544355af7358c0a02e12e5a Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Fri, 17 Jun 2016 16:14:58 +0300 Subject: [PATCH] modal analysis with geometric stiffness --- src/elasticity.jl | 77 ++++++++++++++++++------------------- test/test_modal_analysis.jl | 6 ++- 2 files changed, 42 insertions(+), 41 deletions(-) diff --git a/src/elasticity.jl b/src/elasticity.jl index e373d2d..20a203c 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -9,10 +9,11 @@ type Elasticity <: FieldProblem # these are found from problem.properties for type Problem{Elasticity} formulation :: Symbol finite_strain :: Bool + geometric_stiffness :: Bool end function Elasticity() # formulations: plane_stress, plane_strain, continuum - return Elasticity(:continuum, false) + return Elasticity(:continuum, false, false) end function get_unknown_field_name(problem::Problem{Elasticity}) @@ -356,14 +357,16 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el if haskey(element, "displacement") gradu += element("displacement", ip, time, Val{:Grad}) end - strain = zeros(dim , dim) - strain += 1/2*(gradu' + gradu) + strain = 1/2*(gradu' + gradu) + F = eye(dim) if props.finite_strain F += gradu strain += 1/2*gradu'*gradu end + # material stiffness start + fill!(BL, 0.0) for i=1:nnodes BL[1, 3*(i-1)+1] = F[1,1]*dN[1,i] @@ -386,11 +389,12 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el BL[6, 3*(i-1)+3] = F[3,3]*dN[1,i] + F[3,1]*dN[3,i] end - # material stiffness start + strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; strain[1,2]; strain[2,3]; strain[1,3]] + update!(ip, "strain", time => strain_vec) + # calculate stress E = element("youngs modulus", ip, time) nu = element("poissons ratio", ip, time) - D = 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 @@ -398,51 +402,41 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el 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] - - # calculate stress - strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; strain[1,2]; strain[2,3]; strain[1,3]] stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec) - - # update values to integration point - update!(ip, "strain", time => strain_vec) update!(ip, "stress", time => stress_vec) Km += w*BL'*D*BL - if get_formulation_type(problem) == :incremental - f -= w*BL'*stress_vec - end - # material stiffness end - # geometric stiffness start + if props.geometric_stiffness # take geometric stiffness into account - fill!(BNL, 0.0) - for i=1:size(dN, 2) - BNL[1, 3*(i-1)+1] = dN[1,i] - BNL[2, 3*(i-1)+1] = dN[2,i] - BNL[3, 3*(i-1)+1] = dN[3,i] - BNL[4, 3*(i-1)+2] = dN[1,i] - BNL[5, 3*(i-1)+2] = dN[2,i] - BNL[6, 3*(i-1)+2] = dN[3,i] - BNL[7, 3*(i-1)+3] = dN[1,i] - BNL[8, 3*(i-1)+3] = dN[2,i] - BNL[9, 3*(i-1)+3] = dN[3,i] - end - S3 = zeros(3*dim, 3*dim) - S3[1,1] = stress_vec[1] - S3[2,2] = stress_vec[2] - S3[3,3] = stress_vec[3] - S3[1,2] = S3[2,1] = stress_vec[4] - S3[2,3] = S3[3,2] = stress_vec[5] - S3[1,3] = S3[3,1] = stress_vec[6] - S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3] + fill!(BNL, 0.0) + + for i=1:size(dN, 2) + BNL[1, 3*(i-1)+1] = dN[1,i] + BNL[2, 3*(i-1)+1] = dN[2,i] + BNL[3, 3*(i-1)+1] = dN[3,i] + BNL[4, 3*(i-1)+2] = dN[1,i] + BNL[5, 3*(i-1)+2] = dN[2,i] + BNL[6, 3*(i-1)+2] = dN[3,i] + BNL[7, 3*(i-1)+3] = dN[1,i] + BNL[8, 3*(i-1)+3] = dN[2,i] + BNL[9, 3*(i-1)+3] = dN[3,i] + end + + S3 = zeros(3*dim, 3*dim) + S3[1,1] = stress_vec[1] + S3[2,2] = stress_vec[2] + S3[3,3] = stress_vec[3] + S3[1,2] = S3[2,1] = stress_vec[4] + S3[2,3] = S3[3,2] = stress_vec[5] + S3[1,3] = S3[3,1] = stress_vec[6] + S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3] - if props.finite_strain Kg += w*BNL'*S3*BNL - end - # geometric stiffness end + end # external load start @@ -450,6 +444,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el T = element("displacement load", ip, time) f += w*vec(T*N) end + for i=1:dim if haskey(element, "displacement load $i") b = element("displacement load $i", ip, time) @@ -459,6 +454,10 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el # external load end + if get_formulation_type(problem) == :incremental + f -= w*BL'*stress_vec + end + end return Km, Kg, f diff --git a/test/test_modal_analysis.jl b/test/test_modal_analysis.jl index 61b9176..2d7f4f9 100644 --- a/test/test_modal_analysis.jl +++ b/test/test_modal_analysis.jl @@ -26,7 +26,7 @@ using JuliaFEM.Test "displacement 2" => 0.0, "displacement 3" => 0.0) p1 = Problem(Elasticity, 3) -# p1.properties.finite_strain = true + p1.properties.finite_strain = false p2 = Problem(Dirichlet, p1) push!(p1, e1) push!(p2, e2) @@ -36,7 +36,9 @@ using JuliaFEM.Test call(s1; debug=true) @test isapprox(s1.properties.eigvals, [4/3, 1/3]) + + p1.properties.geometric_stiffness = true s1.properties.geometric_stiffness = true - call(s1) + call(s1; debug=true) @test isapprox(s1.properties.eigvals, [5/3, 2/3]) end