{ "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", "@law σ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": [ "@law function update(ε_in,r_in,d_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": [ "@law function σ(ε_in,r_in,d_in)\n", "\n", " _, _, d_out = update(ε_in,r_in,d_in)\n", "\n", " (1-d_out)*σe(ε_in)\n", "\n", "end\n", "\n", "@law function dσ(dε_in,ε_in,state)\n", "\n", " damaged, r_out, d_out = state\n", "\n", " if ! damaged\n", " return (1-d_out)*σe(dε_in)\n", "\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", "\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,trian,quad,order)\n", "\n", " a(u,v) = u*v\n", " l(v) = v*q\n", " t_Ω = AffineFETerm(a,l,trian,quad)\n", "\n", " V = TestFESpace(\n", " reffe=:Lagrangian, valuetype=Float64, order=order,\n", " triangulation=trian, conformity=:L2)\n", "\n", " U = TrialFESpace(V)\n", " op = AffineFEOperator(U,V,t_Ω)\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", " V = TestFESpace(\n", " reffe=:Lagrangian, valuetype=VectorValue{3,Float64}, order=order,\n", " model=model,\n", " labels = labeling,\n", " dirichlet_tags=[\"supports\"])\n", "\n", " U = TrialFESpace(V)\n", "\n", " trian = Triangulation(model)\n", " degree = 2*order\n", " quad = CellQuadrature(trian,degree)\n", "\n", " r = CellField(r_0,trian,quad)\n", " d = CellField(0.0,trian,quad)\n", "\n", " nls = NLSolver(show_trace=true, method=:newton)\n", " solver = FESolver(nls)\n", "\n", " function step(uh_in,factor)\n", "\n", " b = factor*b_max\n", " res(u,v) = ( ε(v) ⊙ σ(ε(u),r,d) ) - v⋅b\n", " function jac(u,du,v)\n", " state = update(ε(u),r,d)\n", " ε(v) ⊙ dσ(ε(du),ε(u),state)\n", " end\n", " t_Ω = FETerm(res,jac,trian,quad)\n", " op = FEOperator(U,V,t_Ω)\n", "\n", " uh_out, = solve!(uh_in,solver,op)\n", "\n", " update_state_variables!(update,quad,ε(uh_out),r,d)\n", "\n", " uh_out\n", " end\n", "\n", " factors = collect(1:nsteps)*(1/nsteps)\n", " uh = zero(V)\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 = step(uh,factor)\n", " dh = project(d,trian,quad,order)\n", " rh = project(r,trian,quad,order)\n", "\n", " writevtk(\n", " trian,\"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.5.3" }, "kernelspec": { "name": "julia-1.5", "display_name": "Julia 1.5.3", "language": "julia" } }, "nbformat": 4 }