{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Tutorial 11: 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",
""
],
"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.10.8"
},
"kernelspec": {
"name": "julia-1.10",
"display_name": "Julia 1.10.8",
"language": "julia"
}
},
"nbformat": 4
}