{ "cells": [ { "cell_type": "markdown", "source": [ "# Tutorial 10: Isotropic damage model" ], "metadata": {} }, { "cell_type": "markdown", "source": [ " This tutorial is under construction, but the code below is already functional." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "using Gridap\n", "using LinearAlgebra" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "## Model definition\n", "\n", "Elastic branch" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "const E = 3.0e10 # Pa\n", "const ν = 0.3 # dim-less\n", "const λ = (E*ν)/((1+ν)*(1-2*ν))\n", "const μ = E/(2*(1+ν))\n", "σe(ε) = λ*tr(ε)*one(ε) + 2*μ*ε # Pa\n", "τ(ε) = sqrt(ε ⊙ σe(ε)) # Pa^(1/2)" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "Damage" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "const σ_u = 4.0e5 # Pa\n", "const r_0 = σ_u / sqrt(E) # Pa^(1/2)\n", "const H = 0.5 # dim-less\n", "\n", "function d(r)\n", " 1 - q(r)/r\n", "end\n", "\n", "function q(r)\n", " r_0 + H*(r-r_0)\n", "end" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "Update of the state variables" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "function new_state(r_in,d_in,ε_in)\n", " τ_in = τ(ε_in)\n", " if τ_in <= r_in\n", " r_out = r_in\n", " d_out = d_in\n", " damaged = false\n", " else\n", " r_out = τ_in\n", " d_out = d(r_out)\n", " damaged = true\n", " end\n", " damaged, r_out, d_out\n", "end" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "Constitutive law and its linearization" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "function σ(ε_in,r_in,d_in)\n", " _, _, d_out = new_state(r_in,d_in,ε_in)\n", " (1-d_out)*σe(ε_in)\n", "end\n", "\n", "function dσ(dε_in,ε_in,state)\n", " damaged, r_out, d_out = state\n", " if ! damaged\n", " return (1-d_out)*σe(dε_in)\n", " else\n", " c_inc = ((q(r_out) - H*r_out)*(σe(ε_in) ⊙ dε_in))/(r_out^3)\n", " return (1-d_out)*σe(dε_in) - c_inc*σe(ε_in)\n", " end\n", "end" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "max dead load" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "const b_max = VectorValue(0.0,0.0,-(9.81*2.5e3))" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "## L2 projection\n", "form Gauss points to a Lagrangian piece-wise discontinuous space" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "function project(q,model,dΩ,order)\n", " reffe = ReferenceFE(lagrangian,Float64,order)\n", " V = FESpace(model,reffe,conformity=:L2)\n", " a(u,v) = ∫( u*v )*dΩ\n", " l(v) = ∫( v*q )*dΩ\n", " op = AffineFEOperator(a,l,V,V)\n", " qh = solve(op)\n", " qh\n", "end" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "## Main function" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "function main(;n,nsteps)\n", "\n", " r = 12\n", " domain = (0,r,0,1,0,1)\n", " partition = (r*n,n,n)\n", " model = CartesianDiscreteModel(domain,partition)\n", "\n", " labeling = get_face_labeling(model)\n", " add_tag_from_tags!(labeling,\"supportA\",[1,3,5,7,13,15,17,19,25])\n", " add_tag_from_tags!(labeling,\"supportB\",[2,4,6,8,14,16,18,20,26])\n", " add_tag_from_tags!(labeling,\"supports\",[\"supportA\",\"supportB\"])\n", "\n", " order = 1\n", "\n", " reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},order)\n", " V = TestFESpace(model,reffe,labels=labeling,dirichlet_tags=[\"supports\"])\n", " U = TrialFESpace(V)\n", "\n", " degree = 2*order\n", " Ω = Triangulation(model)\n", " dΩ = Measure(Ω,degree)\n", "\n", " r = CellState(r_0,dΩ)\n", " d = CellState(0.0,dΩ)\n", "\n", " nls = NLSolver(show_trace=true, method=:newton)\n", " solver = FESolver(nls)\n", "\n", " function step(uh_in,factor,cache)\n", " b = factor*b_max\n", " res(u,v) = ∫( ε(v) ⊙ (σ∘(ε(u),r,d)) - v⋅b )*dΩ\n", " jac(u,du,v) = ∫( ε(v) ⊙ (dσ∘(ε(du),ε(u),new_state∘(r,d,ε(u)))) )*dΩ\n", " op = FEOperator(res,jac,U,V)\n", " uh_out, cache = solve!(uh_in,solver,op,cache)\n", " update_state!(new_state,r,d,ε(uh_out))\n", " uh_out, cache\n", " end\n", "\n", " factors = collect(1:nsteps)*(1/nsteps)\n", " uh = zero(V)\n", " cache = nothing\n", "\n", " for (istep,factor) in enumerate(factors)\n", "\n", " println(\"\\n+++ Solving for load factor $factor in step $istep of $nsteps +++\\n\")\n", "\n", " uh,cache = step(uh,factor,cache)\n", " dh = project(d,model,dΩ,order)\n", " rh = project(r,model,dΩ,order)\n", "\n", " writevtk(\n", " Ω,\"results_$(lpad(istep,3,'0'))\",\n", " cellfields=[\"uh\"=>uh,\"epsi\"=>ε(uh),\"damage\"=>dh,\n", " \"threshold\"=>rh,\"sigma_elast\"=>σe∘ε(uh)])\n", "\n", " end\n", "\n", "end" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "Run!" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "main(n=6,nsteps=20)" ], "metadata": {}, "execution_count": null }, { "cell_type": "markdown", "source": [ "## Results" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "Animation of the load history using for `main(n=8,nsteps=30)`\n", "![](../assets/isotropic_damage/damage.gif)" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "---\n", "\n", "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" ], "metadata": {} } ], "nbformat_minor": 3, "metadata": { "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.2" }, "kernelspec": { "name": "julia-1.6", "display_name": "Julia 1.6.2", "language": "julia" } }, "nbformat": 4 }