Giter Club home page Giter Club logo

femmaterials.jl's Introduction

FEMMaterials.jl

Computational material model bindings to JuliaFEM. The JuliaFEM compability of the materials is to be tested in this package. Here is one usage example. Note this README file is generated edit readme_header.txt instead

3D Beam with ideal plastic material model from Materials.jl

using JuliaFEM, FEMMaterials, Materials, FEMBase, LinearAlgebra, Plots
import FEMMaterials: Continuum3D, MecaMatSo
pkg_dir = dirname(dirname(pathof(FEMMaterials)))

Let's read the discretized geometry and create boundary conditions

File plactic_beam.inp is created with 3rd party meshing tool. File contains surface sets BC1, BC2 and PRESSURE as well element set Body1

mesh = abaqus_read_mesh(joinpath(pkg_dir,"examples","data_3dbeam","plastic_beam.inp"))
beam_elements = create_elements(mesh, "Body1")
bc_elements_1 = create_nodal_elements(mesh, "BC1")
bc_elements_2 = create_nodal_elements(mesh, "BC2")
trac_elements = create_surface_elements(mesh, "PRESSURE")

for j in 1:3
    update!(bc_elements_1, "displacement $j", 0.0)
end
update!(bc_elements_2, "displacement 1", 0.0)
update!(bc_elements_2, "displacement 2", 0.0)
update!(trac_elements, "surface pressure", 0.0 => 0.00)
update!(trac_elements, "surface pressure", 1.0 => 2.70)
trac = Problem(Continuum3D, "traction", 3)
bc = Problem(Dirichlet, "fix displacement", 3, "displacement")
add_elements!(trac, trac_elements)
add_elements!(bc, bc_elements_1)
add_elements!(bc, bc_elements_2)

Next, we set the material properties for each element set

In this example we only have the one element set Body1 in variable beam_elements

update!(beam_elements, "youngs_modulus", 200.0e3)
update!(beam_elements, "poissons_ratio", 0.3)
update!(beam_elements, "yield_stress", 100.0)

beam = Problem(Continuum3D, "plastic beam", 3)
beam.properties.material_model = :IdealPlastic
add_elements!(beam, beam_elements)

And next, we setup the analysis

t0, t1, and dt are the start time, end time, and time step respectively.

analysis = Analysis(MecaMatSo, "solve problem")
analysis.properties.max_iterations = 50
analysis.properties.t0 = 0.0
analysis.properties.t1 = 1.0
analysis.properties.dt = 0.05

Writing the results needs to be setup as well

xdmf = Xdmf("3dbeam_results_output"; overwrite=true)
add_results_writer!(analysis, xdmf)

Finally, adding the problems together and running the analysis

All earlier defined problems are added together. Also result file need to be closed to flush everything from the writing buffer to the file.

add_problems!(analysis, beam, trac, bc)
run!(analysis)
close(xdmf)

The first post-processing step is to calculate maximum von Mises stresses

vmis contains all integration points stresses and vmis_ just the maximum. Let's plot the maximum von Mises stress as a function of time

tim = 0.0:0.05:1.0
vmis_ = []
for t in tim
    vmis = []
    for element in beam_elements
        for ip in get_integration_points(element)
            s11, s22, s33, s12, s23, s31 = ip("stress", t)
            push!(vmis, sqrt(1/2*((s11-s22)^2 + (s22-s33)^2 + (s33-s11)^2 + 6*(s12^2+s23^2+s31^2))))
        end
    end
    push!(vmis_, maximum(vmis))
end
plot(tim,vmis_)
png("max_vonmises_stress_as_a_function_of_time")

![max_vonmises_stress_as_a_function_of_time][vonmises] [vonmises]: https://raw.githubusercontent.com/JuliaFEM/FEMMaterials.jl/master/notebooks/max_vonmises_stress_as_a_function_of_time.png

The second post-processing step is to collect displacements

Here as an example node number 96 second degree of freedom displacement is extracted.

u2_96 = []
for t in tim
    push!(u2_96, beam("displacement", t)[96][2])
end
plot(tim,u2_96)
png("node_96_displacement_as_a_function_of_time")

![node_96_displacement_as_a_function_of_time][displacement] [displacement]: https://raw.githubusercontent.com/JuliaFEM/FEMMaterials.jl/master/notebooks/node_96_displacement_as_a_function_of_time.png

Page generated at 2019-10-07T16:54:37.559.

This page was generated using Literate.jl.

femmaterials.jl's People

Contributors

ahojukka5 avatar ivanyashchuk avatar jvaara avatar kristofferc avatar terofrondelius avatar

Watchers

avatar avatar avatar avatar

Forkers

ivanyashchuk juliatagbot

femmaterials.jl's Issues

MecaMatSo vector notations for symmetric tensors

MecaMatSo's assembler uses ordering 11 22 33 12 23 13
https://github.com/JuliaFEM/FEMMaterials.jl/blob/master/src/mecamatso.jl#L101

while Voigt notation is 11 22 33 23 13 12

MFront library uses Mandel notation

In Voigt ordering the strain-displacement matrix is then formed as

fill!(BL, 0.0)
for i=1:nnodes
    BL[1, 3*(i-1)+1] = dN[1,i]
    BL[2, 3*(i-1)+2] = dN[2,i]
    BL[3, 3*(i-1)+3] = dN[3,i]
    BL[6, 3*(i-1)+1] = dN[2,i]
    BL[6, 3*(i-1)+2] = dN[1,i]
    BL[4, 3*(i-1)+2] = dN[3,i]
    BL[4, 3*(i-1)+3] = dN[2,i]
    BL[5, 3*(i-1)+1] = dN[3,i]
    BL[5, 3*(i-1)+3] = dN[1,i]
end

MecaMatSo result fluctuation

The test_mecamatso.jl test has some fluctuation of the displacement field during plastic flow causing also the stress state to be not exactly uniaxial. The relative tolerance of stresses was increased to 1e-1 to make the test pass. This is possibly due to

  1. Lack of good initial guess for Newton's iterations
  2. The penalty method used to enforce the forced displacements

Below are some prints of the stress development during the first yield. The elastic border is reached at time 0.5 and the stress state should remain uniaxial (100 MPa in 3,3 - component and zero everywhere else)

julia> ip("stress", 0.5)
3×3 Array{Float64,2}:
  1.70221e-14  -8.47033e-16   -3.25782e-16
 -8.47033e-16   1.52458e-14   -7.81877e-16
 -3.25782e-16  -7.81877e-16  100.0

julia> ip("stress", 0.6)
3×3 Array{Float64,2}:
  0.157087    0.0533097   -0.058128
  0.0533097   0.0216236   -0.0310419
 -0.058128   -0.0310419  100.089

julia> ip("stress", 0.7)
3×3 Array{Float64,2}:
 -0.150201   -0.153023    0.0919305
 -0.153023   -0.0509997   0.106634
  0.0919305   0.106634   99.8987

julia> ip("stress", 0.8)
3×3 Array{Float64,2}:
  0.0745154   0.119496   -0.0302556
  0.119496    0.234576   -0.140247
 -0.0302556  -0.140247  100.154

julia> ip("stress", 0.9)
3×3 Array{Float64,2}:
 -0.363192    0.208532   -0.0921637
  0.208532   -0.78912     0.0801902
 -0.0921637   0.0801902  99.4223

julia> ip("stress", 1.0)
3×3 Array{Float64,2}:
  0.478094   -0.344551    0.0928852
 -0.344551    0.560487    0.04186
  0.0928852   0.04186   100.517

Perfomance tests

As @ahojukka5 pointed out in #7 performance is a critical thing in FEM material models. Thus, we need performance tests. I think that #6 could be the first test to add.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.