{ "cells": [ { "cell_type": "markdown", "source": [ "# Computational homogenization\n", "\n", "![](rve_homogenization.png)\n", "\n", "*Figure 1*: von Mises stress in an RVE with 5 stiff inclusions embedded in a softer matrix\n", "material that is loaded in shear. The problem is solved by using homogeneous Dirichlet\n", "boundary conditions (left) and (strong) periodic boundary conditions (right)." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "## Introduction\n", "\n", "In this example we will solve the Representative Volume Element (RVE) problem for\n", "computational homogenization of linear elasticity and compute the effective/homogenized\n", "stiffness of an RVE with 5 stiff circular inclusions embedded in a softer matrix material\n", "(see Figure 1).\n", "\n", "It is possible to obtain upper and lower bounds on the stiffness analytically, see for\n", "example [Rule of mixtures](https://en.wikipedia.org/wiki/Rule_of_mixtures). An upper\n", "bound is obtained from the Voigt model, where the *strain* is assumed to be the same in\n", "the two constituents,\n", "\n", "$$\n", "\\mathsf{E}_\\mathrm{Voigt} = v_\\mathrm{m} \\mathsf{E}_\\mathrm{m} +\n", "(1 - v_\\mathrm{m}) \\mathsf{E}_\\mathrm{i}\n", "$$\n", "\n", "where $v_\\mathrm{m}$ is the volume fraction of the matrix material, and where\n", "$\\mathsf{E}_\\mathrm{m}$ and $\\mathsf{E}_\\mathrm{i}$ are the individual stiffness for\n", "the matrix material and the inclusions, respectively. The lower bound is obtained from\n", "the Reuss model, where the *stress* is assumed to be the same in the two constituents,\n", "\n", "$$\n", "\\mathsf{E}_\\mathrm{Reuss} = \\left(v_\\mathrm{m} \\mathsf{E}_\\mathrm{m}^{-1} +\n", "(1 - v_\\mathrm{m}) \\mathsf{E}_\\mathrm{i}^{-1} \\right)^{-1}.\n", "$$\n", "\n", "However, neither of these assumptions are, in general, very close to the \"truth\" which is\n", "why it is of interest to computationally find the homogenized properties for a given RVE.\n", "\n", "The canonical version of the RVE problem can be formulated as follows:\n", "For given homogenized field $\\bar{\\boldsymbol{u}}$, $\\bar{\\boldsymbol{\\varepsilon}} =\n", "\\boldsymbol{\\varepsilon}[\\bar{\\boldsymbol{u}}]$, find $\\boldsymbol{u} \\in\n", "\\mathbb{U}_\\Box$, $\\boldsymbol{t} \\in \\mathbb{T}_\\Box$ such that\n", "\n", "$$\n", "\\frac{1}{|\\Omega_\\Box|} \\int_{\\Omega_\\Box}\\boldsymbol{\\varepsilon}[\\delta\\boldsymbol{u}]\n", ": \\mathsf{E} : \\boldsymbol{\\varepsilon}[\\boldsymbol{u}]\\ \\mathrm{d}\\Omega\n", "- \\frac{1}{|\\Omega_\\Box|} \\int_{\\Gamma_\\Box}\\delta \\boldsymbol{u} \\cdot\n", "\\boldsymbol{t}\\ \\mathrm{d}\\Gamma = 0 \\quad\n", "\\forall \\delta \\boldsymbol{u} \\in \\mathbb{U}_\\Box,\\quad (1\\mathrm{a})\\\\\n", "- \\frac{1}{|\\Omega_\\Box|} \\int_{\\Gamma_\\Box}\\delta \\boldsymbol{t} \\cdot\n", "\\boldsymbol{u}\\ \\mathrm{d}\\Gamma = - \\bar{\\boldsymbol{\\varepsilon}} :\n", "\\left[ \\frac{1}{|\\Omega_\\Box|} \\int_{\\Gamma_\\Box}\\delta \\boldsymbol{t} \\otimes\n", "[\\boldsymbol{x} - \\bar{\\boldsymbol{x}}]\\ \\mathrm{d}\\Gamma \\right]\n", "\\quad \\forall \\delta \\boldsymbol{t} \\in \\mathbb{T}_\\Box, \\quad (1\\mathrm{b})\n", "$$\n", "\n", "where $\\boldsymbol{u} = \\bar{\\boldsymbol{\\varepsilon}} \\cdot [\\boldsymbol{x} -\n", "\\bar{\\boldsymbol{x}}] + \\boldsymbol{u}^\\mu$, where $\\Omega_\\Box$ and $|\\Omega_\\Box|$\n", "are the domain and volume of the RVE, where $\\Gamma_\\Box$ is the boundary, and where\n", "$\\mathbb{U}_\\Box$, $\\mathbb{T}_\\Box$ are set of \"sufficiently regular\" functions\n", "defined on the RVE.\n", "\n", "This system is not solvable without introducing extra restrictions on $\\mathbb{U}_\\Box$,\n", "$\\mathbb{T}_\\Box$. In this example we will consider the common cases of Dirichlet\n", "boundary conditions and (strong) periodic boundary conditions.\n", "\n", "**Dirichlet boundary conditions**\n", "\n", "We can introduce the more restrictive sets of $\\mathbb{U}_\\Box$:\n", "\n", "$$\n", "\\begin{align*}\n", "\\mathbb{U}_\\Box^\\mathrm{D} &:= \\left\\{\\boldsymbol{u} \\in \\mathbb{U}_\\Box|\\ \\boldsymbol{u}\n", "= \\bar{\\boldsymbol{\\varepsilon}} \\cdot [\\boldsymbol{x} - \\bar{\\boldsymbol{x}}]\n", "\\ \\mathrm{on}\\ \\Gamma_\\Box\\right\\},\\\\\n", "\\mathbb{U}_\\Box^{\\mathrm{D},0} &:= \\left\\{\\boldsymbol{u} \\in \\mathbb{U}_\\Box|\\ \\boldsymbol{u}\n", "= \\boldsymbol{0}\\ \\mathrm{on}\\ \\Gamma_\\Box\\right\\},\n", "\\end{align*}\n", "$$\n", "\n", "and use these as trial and test sets to obtain a solvable RVE problem pertaining to\n", "Dirichlet boundary conditions. Eq. $(1\\mathrm{b})$ is trivially fulfilled, the second\n", "term of Eq. $(1\\mathrm{a})$ vanishes, and we are left with the following problem:\n", "Find $\\boldsymbol{u} \\in \\mathbb{U}_\\Box^\\mathrm{D}$ that solve\n", "\n", "$$\n", "\\frac{1}{|\\Omega_\\Box|} \\int_{\\Omega_\\Box}\\boldsymbol{\\varepsilon}[\\delta\\boldsymbol{u}]\n", ": \\mathsf{E} : \\boldsymbol{\\varepsilon}[\\boldsymbol{u}]\\ \\mathrm{d}\\Omega = 0\n", "\\quad \\forall \\delta \\boldsymbol{u} \\in \\mathbb{U}_\\Box^{\\mathrm{D},0}.\n", "$$\n", "\n", "Note that, since $\\boldsymbol{u} = \\bar{\\boldsymbol{\\varepsilon}} \\cdot [\\boldsymbol{x} -\n", "\\bar{\\boldsymbol{x}}] + \\boldsymbol{u}^\\mu$, this problem is equivalent to solving for\n", "$\\boldsymbol{u}^\\mu \\in \\mathbb{U}_\\Box^{\\mathrm{D},0}$, which is what we will do in\n", "the implementation.\n", "\n", "**Periodic boundary conditions**\n", "\n", "The RVE problem pertaining to periodic boundary conditions is obtained by restricting\n", "$\\boldsymbol{u}^\\mu$ to be periodic, and $\\boldsymbol{t}$ anti-periodic across the\n", "RVE. Similarly as for Dirichlet boundary conditions, Eq. $(1\\mathrm{b})$ is directly\n", "fulfilled, and the second term in Eq. $(1\\mathrm{a})$ vanishes, with these restrictions,\n", "and we are left with the following problem:\n", "Find $\\boldsymbol{u}^\\mu \\in \\mathbb{U}_\\Box^{\\mathrm{P},0}$ such that\n", "\n", "$$\n", "\\frac{1}{|\\Omega_\\Box|} \\int_{\\Omega_\\Box}\\boldsymbol{\\varepsilon}[\\delta\\boldsymbol{u}]\n", ": \\mathsf{E} : (\\bar{\\boldsymbol{\\varepsilon}} + \\boldsymbol{\\varepsilon}\n", "[\\boldsymbol{u}^\\mu])\\ \\mathrm{d}\\Omega = 0\n", "\\quad \\forall \\delta \\boldsymbol{u} \\in \\mathbb{U}_\\Box^{\\mathrm{P},0},\n", "$$\n", "\n", "where\n", "\n", "$$\n", "\\mathbb{U}_\\Box^{\\mathrm{P},0} := \\left\\{\\boldsymbol{u} \\in \\mathbb{U}_\\Box|\n", "\\ [\\![ \\boldsymbol{u} ]\\!]_\\Box = \\boldsymbol{0}\n", "\\ \\mathrm{on}\\ \\Gamma_\\Box^+\\right\\}\n", "$$\n", "\n", "where $[\\![ \\bullet ]\\!]_\\Box = \\bullet(\\boldsymbol{x}^+) -\n", "\\bullet(\\boldsymbol{x}^-)$ defines the \"jump\" over the RVE, i.e. the difference between\n", "the value on the image part $\\Gamma_\\Box^+$ (coordinate $\\boldsymbol{x}^+$) and the\n", "mirror part $\\Gamma_\\Box^-$ (coordinate $\\boldsymbol{x}^-$) of the boundary.\n", "To make sure this restriction holds in a strong sense we need a periodic mesh.\n", "\n", "Note that it would be possible to solve for the total $\\boldsymbol{u}$ directly by\n", "instead enforcing the jump to be equal to the jump in the macroscopic part,\n", "$\\boldsymbol{u}^\\mathrm{M}$, i.e.\n", "\n", "$$\n", "[\\![ \\boldsymbol{u} ]\\!]_\\Box =\n", "[\\![ \\boldsymbol{u}^\\mathrm{M} ]\\!]_\\Box =\n", "[\\![ \\bar{\\boldsymbol{\\varepsilon}} \\cdot [\\boldsymbol{x} - \\bar{\\boldsymbol{x}}]\n", "]\\!]_\\Box =\n", "\\bar{\\boldsymbol{\\varepsilon}} \\cdot [\\boldsymbol{x}^+ - \\boldsymbol{x}^-].\n", "$$\n", "\n", "**Homogenization of effective properties**\n", "\n", "In general it is necessary to compute the homogenized stress and the stiffness on the fly,\n", "but since we in this example consider linear elasticity it is possible to compute the\n", "effective properties once and for all for a given RVE configuration. We do this by\n", "computing sensitivity fields for every independent strain component (6 in 3D, 3 in 2D).\n", "Thus, for a 2D problem, as in the implementation below, we compute sensitivities\n", "$\\hat{\\boldsymbol{u}}_{11}$, $\\hat{\\boldsymbol{u}}_{22}$, and\n", "$\\hat{\\boldsymbol{u}}_{12} = \\hat{\\boldsymbol{u}}_{21}$ by using\n", "\n", "$$\n", "\\bar{\\boldsymbol{\\varepsilon}} = \\begin{pmatrix}1 & 0\\\\ 0 & 0\\end{pmatrix}, \\quad\n", "\\bar{\\boldsymbol{\\varepsilon}} = \\begin{pmatrix}0 & 0\\\\ 0 & 1\\end{pmatrix}, \\quad\n", "\\bar{\\boldsymbol{\\varepsilon}} = \\begin{pmatrix}0 & 0.5\\\\ 0.5 & 0\\end{pmatrix}\n", "$$\n", "\n", "as the input to the RVE problem. When the sensitivies are solved we can compute the\n", "entries of the homogenized stiffness as follows\n", "\n", "$$\n", "\\mathsf{E}_{ijkl} = \\frac{\\partial\\ \\bar{\\sigma}_{ij}}{\\partial\\ \\bar{\\varepsilon}_{kl}}\n", "= \\bar{\\sigma}_{ij}(\\hat{\\boldsymbol{u}}_{kl}),\n", "$$\n", "\n", "where the homogenized stress, $\\bar{\\boldsymbol{\\sigma}}(\\boldsymbol{u})$, is computed\n", "as the volume average of the stress in the RVE, i.e.\n", "\n", "$$\n", "\\bar{\\boldsymbol{\\sigma}}(\\boldsymbol{u}) :=\n", "\\frac{1}{|\\Omega_\\Box|} \\int_{\\Omega_\\Box} \\boldsymbol{\\sigma}\\ \\mathrm{d}\\Omega =\n", "\\frac{1}{|\\Omega_\\Box|} \\int_{\\Omega_\\Box}\n", "\\mathsf{E} : \\boldsymbol{\\varepsilon}[\\boldsymbol{u}]\\ \\mathrm{d}\\Omega.\n", "$$" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "## Commented program\n", "\n", "Now we will see how this can be implemented in `Ferrite`. What follows is a program\n", "with comments in between which describe the different steps." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "using Ferrite, SparseArrays, LinearAlgebra" ], "metadata": {}, "execution_count": 1 }, { "cell_type": "markdown", "source": [ "We first load the mesh file [`periodic-rve.msh`](periodic-rve.msh)\n", "([`periodic-rve-coarse.msh`](periodic-rve-coarse.msh) for a coarser mesh). The mesh is\n", "generated with [`gmsh`](https://gmsh.info/), and we read it in as a `Ferrite` grid using\n", "the [`FerriteGmsh`](https://github.com/Ferrite-FEM/FerriteGmsh.jl) package:" ], "metadata": {} }, { "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Info : Reading 'periodic-rve-coarse.msh'...\n", "Info : 38 entities\n", "Info : 112 nodes\n", "Info : 222 elements\n", "Info : Done reading 'periodic-rve-coarse.msh'\n" ] }, { "output_type": "execute_result", "data": { "text/plain": "Grid{2, Triangle, Float64} with 186 Triangle cells and 112 nodes" }, "metadata": {}, "execution_count": 2 } ], "cell_type": "code", "source": [ "using FerriteGmsh\n", "# grid = saved_file_to_grid(\"periodic-rve.msh\")\n", "grid = saved_file_to_grid(\"periodic-rve-coarse.msh\")" ], "metadata": {}, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "Next we construct the interpolation and quadrature rule, and combining them into\n", "cellvalues as usual:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "dim = 2\n", "ip = Lagrange{dim, RefTetrahedron, 1}()\n", "qr = QuadratureRule{dim, RefTetrahedron}(2)\n", "cellvalues = CellVectorValues(qr, ip);" ], "metadata": {}, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "We define a dof handler with a displacement field `:u`:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "dh = DofHandler(grid)\n", "add!(dh, :u, 2)\n", "close!(dh);" ], "metadata": {}, "execution_count": 4 }, { "cell_type": "markdown", "source": [ "Now we need to define boundary conditions. As discussed earlier we will solve the problem\n", "using (i) homogeneous Dirichlet boundary conditions, and (ii) periodic Dirichlet boundary\n", "conditions. We construct two different constraint handlers, one for each case. The\n", "`Dirichlet` boundary condition we have seen in many other examples. Here we simply\n", "define the condition that the field, `:u`, should have both components prescribed to `0`\n", "on the full boundary:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "ch_dirichlet = ConstraintHandler(dh)\n", "dirichlet = Dirichlet(\n", " :u,\n", " union(getfaceset.(Ref(grid), [\"left\", \"right\", \"top\", \"bottom\"])...),\n", " (x, t) -> [0, 0],\n", " [1, 2]\n", ")\n", "add!(ch_dirichlet, dirichlet)\n", "close!(ch_dirichlet)\n", "update!(ch_dirichlet, 0.0)" ], "metadata": {}, "execution_count": 5 }, { "cell_type": "markdown", "source": [ "For periodic boundary conditions we use the `PeriodicDirichlet` constraint type,\n", "which is very similar to the `Dirichlet` type, but instead of a passing a faceset we pass\n", "a vector with \"face pairs\", i.e. the mapping between mirror and image parts of the\n", "boundary. In this example the `\"left\"` and `\"bottom\"` boundaries are mirrors, and the\n", "`\"right\"` and `\"top\"` boundaries are the mirrors." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "ch_periodic = ConstraintHandler(dh);\n", "periodic = PeriodicDirichlet(\n", " :u,\n", " [\"left\" => \"right\", \"bottom\" => \"top\"],\n", " [1, 2]\n", ")\n", "add!(ch_periodic, periodic)\n", "close!(ch_periodic)\n", "update!(ch_periodic, 0.0)" ], "metadata": {}, "execution_count": 6 }, { "cell_type": "markdown", "source": [ "This will now constrain any degrees of freedom located on the mirror boundaries to\n", "the matching degree of freedom on the image boundaries. Internally this will create\n", "a number of `AffineConstraint`s of the form `u_i = 1 * u_j + 0`:\n", "```julia\n", "a = AffineConstraint(u_m, [u_i => 1], 0)\n", "```\n", "where `u_m` is the degree of freedom on the mirror and `u_i` the matching one on the\n", "image part. `PeriodicDirichlet` is thus simply just a more convenient way of\n", "constructing such affine constraints since it computes the degree of freedom mapping\n", "automatically.\n", "\n", "To simplify things we group the constraint handlers into a named tuple" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "ch = (dirichlet = ch_dirichlet, periodic = ch_periodic);" ], "metadata": {}, "execution_count": 7 }, { "cell_type": "markdown", "source": [ "We can now construct the sparse matrix. Note that, since we are using affine constraints,\n", "which need to modify the matrix sparsity pattern in order to account for the constraint\n", "equations, we construct the matrix for the periodic case by passing both the dof handler\n", "and the constraint handler." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "K = (\n", " dirichlet = create_sparsity_pattern(dh),\n", " periodic = create_sparsity_pattern(dh, ch.periodic),\n", ");" ], "metadata": {}, "execution_count": 8 }, { "cell_type": "markdown", "source": [ "We define the fourth order elasticity tensor for the matrix material, and define the\n", "inclusions to have 10 times higher stiffness" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "λ, μ = 1e10, 7e9 # Lamé parameters\n", "δ(i,j) = i == j ? 1.0 : 0.0\n", "Em = SymmetricTensor{4, 2}(\n", " (i,j,k,l) -> λ * δ(i,j) * δ(k,l) + μ * (δ(i,k) * δ(j,l) + δ(i,l) * δ(j,k))\n", ")\n", "Ei = 10 * Em;" ], "metadata": {}, "execution_count": 9 }, { "cell_type": "markdown", "source": [ "As mentioned above, in order to compute the apparent/homogenized stiffness we will solve\n", "the problem repeatedly with different macroscale strain tensors to compute the sensitvity\n", "of the homogenized stress, $\\bar{\\boldsymbol{\\sigma}}$, w.r.t. the macroscopic strain,\n", "$\\bar{\\boldsymbol{\\varepsilon}}$. The corresponding unit strains are defined below,\n", "and will result in three different right-hand-sides:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "εᴹ = [\n", " SymmetricTensor{2,2}([1.0 0.0; 0.0 0.0]), # ε_11 loading\n", " SymmetricTensor{2,2}([0.0 0.0; 0.0 1.0]), # ε_22 loading\n", " SymmetricTensor{2,2}([0.0 0.5; 0.5 0.0]), # ε_12/ε_21 loading\n", "];" ], "metadata": {}, "execution_count": 10 }, { "cell_type": "markdown", "source": [ "The assembly function is nothing strange, and in particular there is no impact from the\n", "choice of boundary conditions, so the same function can be used for both cases. Since\n", "we want to solve the system 3 times, once for each macroscopic strain component, we\n", "assemble 3 right-hand-sides." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "function doassemble!(cellvalues::CellVectorValues, K::SparseMatrixCSC, dh::DofHandler, εᴹ)\n", "\n", " n_basefuncs = getnbasefunctions(cellvalues)\n", " ndpc = ndofs_per_cell(dh)\n", " Ke = zeros(ndpc, ndpc)\n", " fe = zeros(ndpc, length(εᴹ))\n", " f = zeros(ndofs(dh), length(εᴹ))\n", " assembler = start_assemble(K)\n", "\n", " for cell in CellIterator(dh)\n", "\n", " E = cellid(cell) in getcellset(dh.grid, \"inclusions\") ? Ei : Em\n", " reinit!(cellvalues, cell)\n", " fill!(Ke, 0)\n", " fill!(fe, 0)\n", "\n", " for q_point in 1:getnquadpoints(cellvalues)\n", " dΩ = getdetJdV(cellvalues, q_point)\n", " for i in 1:n_basefuncs\n", " δεi = shape_symmetric_gradient(cellvalues, q_point, i)\n", " for j in 1:n_basefuncs\n", " δεj = shape_symmetric_gradient(cellvalues, q_point, j)\n", " Ke[i, j] += (δεi ⊡ E ⊡ δεj) * dΩ\n", " end\n", " for (rhs, ε) in enumerate(εᴹ)\n", " σᴹ = E ⊡ ε\n", " fe[i, rhs] += ( - δεi ⊡ σᴹ) * dΩ\n", " end\n", " end\n", " end\n", "\n", " cdofs = celldofs(cell)\n", " assemble!(assembler, cdofs, Ke)\n", " f[cdofs, :] .+= fe\n", " end\n", " return f\n", "end;" ], "metadata": {}, "execution_count": 11 }, { "cell_type": "markdown", "source": [ "We can now assemble the system. The assembly function modifies the matrix in-place, but\n", "return the right hand side(s) which we collect in another named tuple." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "rhs = (\n", " dirichlet = doassemble!(cellvalues, K.dirichlet, dh, εᴹ),\n", " periodic = doassemble!(cellvalues, K.periodic, dh, εᴹ),\n", ");" ], "metadata": {}, "execution_count": 12 }, { "cell_type": "markdown", "source": [ "The next step is to solve the systems. Since application of boundary conditions, using\n", "the `apply!` function, modifies both the matrix and the right hand sides we can\n", "not use it directly in this case since we want to reuse the matrix again for the next\n", "right hand sides. We could of course re-assemble the matrix for every right hand side,\n", "but that would not be very efficient. Instead we will use the `get_rhs_data`\n", "function, together with `apply_rhs!` in a later step. This will extract the\n", "necessary data from the matrix such that we can apply it for all the different right\n", "hand sides. Note that we call `apply!` with just the matrix and no right hand side." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "rhsdata = (\n", " dirichlet = get_rhs_data(ch.dirichlet, K.dirichlet),\n", " periodic = get_rhs_data(ch.periodic, K.periodic),\n", ")\n", "\n", "apply!(K.dirichlet, ch.dirichlet)\n", "apply!(K.periodic, ch.periodic)" ], "metadata": {}, "execution_count": 13 }, { "cell_type": "markdown", "source": [ "We can now solve the problem(s). Note that we only use `apply_rhs!` in the loops below.\n", "The boundary conditions are already applied to the matrix above, so we only need to\n", "modify the right hand side." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "u = (\n", " dirichlet = Vector{Float64}[],\n", " periodic = Vector{Float64}[],\n", ")\n", "\n", "for i in 1:size(rhs.dirichlet, 2)\n", " rhs_i = @view rhs.dirichlet[:, i] # Extract this RHS\n", " apply_rhs!(rhsdata.dirichlet, rhs_i, ch.dirichlet) # Apply BC\n", " u_i = cholesky(Symmetric(K.dirichlet)) \\ rhs_i # Solve\n", " apply!(u_i, ch.dirichlet) # Apply BC on the solution\n", " push!(u.dirichlet, u_i) # Save the solution vector\n", "end\n", "\n", "for i in 1:size(rhs.periodic, 2)\n", " rhs_i = @view rhs.periodic[:, i] # Extract this RHS\n", " apply_rhs!(rhsdata.periodic, rhs_i, ch.periodic) # Apply BC\n", " u_i = cholesky(Symmetric(K.periodic)) \\ rhs_i # Solve\n", " apply!(u_i, ch.periodic) # Apply BC on the solution\n", " push!(u.periodic, u_i) # Save the solution vector\n", "end" ], "metadata": {}, "execution_count": 14 }, { "cell_type": "markdown", "source": [ "When the solution(s) are known we can compute the averaged stress,\n", "$\\bar{\\boldsymbol{\\sigma}}$ in the RVE. We define a function that does this, and also\n", "returns the von Mise stress in every quadrature point for visualization." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "function compute_stress(cellvalues::CellVectorValues, dh::DofHandler, u, εᴹ)\n", " σvM_qpdata = zeros(getnquadpoints(cellvalues), getncells(dh.grid))\n", " σ̄Ω = zero(SymmetricTensor{2,2})\n", " Ω = 0.0 # Total volume\n", " for cell in CellIterator(dh)\n", " E = cellid(cell) in getcellset(dh.grid, \"inclusions\") ? Ei : Em\n", " reinit!(cellvalues, cell)\n", " for q_point in 1:getnquadpoints(cellvalues)\n", " dΩ = getdetJdV(cellvalues, q_point)\n", " εμ = function_symmetric_gradient(cellvalues, q_point, u[celldofs(cell)])\n", " σ = E ⊡ (εᴹ + εμ)\n", " σvM_qpdata[q_point, cellid(cell)] = sqrt(3/2 * dev(σ) ⊡ dev(σ))\n", " Ω += dΩ # Update total volume\n", " σ̄Ω += σ * dΩ # Update integrated stress\n", " end\n", " end\n", " σ̄ = σ̄Ω / Ω\n", " return σvM_qpdata, σ̄\n", "end;" ], "metadata": {}, "execution_count": 15 }, { "cell_type": "markdown", "source": [ "We now compute the homogenized stress and von Mise stress for all cases" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "σ̄ = (\n", " dirichlet = SymmetricTensor{2,2}[],\n", " periodic = SymmetricTensor{2,2}[],\n", ")\n", "σ = (\n", " dirichlet = Vector{Float64}[],\n", " periodic = Vector{Float64}[],\n", ")\n", "\n", "projector = L2Projector(ip, grid)\n", "\n", "for i in 1:3\n", " σ_qp, σ̄_i = compute_stress(cellvalues, dh, u.dirichlet[i], εᴹ[i])\n", " proj = project(projector, σ_qp, qr; project_to_nodes=false)\n", " push!(σ.dirichlet, proj)\n", " push!(σ̄.dirichlet, σ̄_i)\n", "end\n", "\n", "for i in 1:3\n", " σ_qp, σ̄_i = compute_stress(cellvalues, dh, u.periodic[i], εᴹ[i])\n", " proj = project(projector, σ_qp, qr; project_to_nodes=false)\n", " push!(σ.periodic, proj)\n", " push!(σ̄.periodic, σ̄_i)\n", "end" ], "metadata": {}, "execution_count": 16 }, { "cell_type": "markdown", "source": [ "The remaining thing is to compute the homogenized stiffness. As mentioned in the\n", "introduction we can find all the components from the average stress of the sensitivity\n", "fields that we have solved for\n", "\n", "$$\n", "\\mathsf{E}_{ijkl} = \\bar{\\sigma}_{ij}(\\hat{\\boldsymbol{u}}_{kl}).\n", "$$\n", "\n", "So we have now already computed all the components, and just need to gather the data in\n", "a fourth order tensor:" ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "E_dirichlet = SymmetricTensor{4,2}((i, j, k, l) -> begin\n", " if k == l == 1\n", " σ̄.dirichlet[1][i, j] # ∂σ∂ε_**11\n", " elseif k == l == 2\n", " σ̄.dirichlet[2][i, j] # ∂σ∂ε_**22\n", " else\n", " σ̄.dirichlet[3][i, j] # ∂σ∂ε_**12 and ∂σ∂ε_**21\n", " end\n", "end)\n", "\n", "E_periodic = SymmetricTensor{4,2}((i, j, k, l) -> begin\n", " if k == l == 1\n", " σ̄.periodic[1][i, j]\n", " elseif k == l == 2\n", " σ̄.periodic[2][i, j]\n", " else\n", " σ̄.periodic[3][i, j]\n", " end\n", "end);" ], "metadata": {}, "execution_count": 17 }, { "cell_type": "markdown", "source": [ "We can check that the result are what we expect, namely that the stiffness with Dirichlet\n", "boundary conditions is higher than when using periodic boundary conditions, and that\n", "the Reuss assumption is an lower bound, and the Voigt assumption a upper bound. We first\n", "compute the volume fraction of the matrix, and then the Voigt and Reuss bounds:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "0.64796265456868" }, "metadata": {}, "execution_count": 18 } ], "cell_type": "code", "source": [ "function matrix_volume_fraction(grid, cellvalues)\n", " V = 0.0 # Total volume\n", " Vm = 0.0 # Volume of the matrix\n", " for c in CellIterator(grid)\n", " reinit!(cellvalues, c)\n", " is_matrix = !(cellid(c) in getcellset(grid, \"inclusions\"))\n", " for qp in 1:getnquadpoints(cellvalues)\n", " dΩ = getdetJdV(cellvalues, qp)\n", " V += dΩ\n", " if is_matrix\n", " Vm += dΩ\n", " end\n", " end\n", " end\n", " return Vm / V\n", "end\n", "\n", "vm = matrix_volume_fraction(grid, cellvalues)" ], "metadata": {}, "execution_count": 18 }, { "outputs": [], "cell_type": "code", "source": [ "E_voigt = vm * Em + (1-vm) * Ei\n", "E_reuss = inv(vm * inv(Em) + (1-vm) * inv(Ei));" ], "metadata": {}, "execution_count": 19 }, { "cell_type": "markdown", "source": [ "We can now compare the different computed stiffness tensors. We expect\n", "$E_\\mathrm{Reuss} \\leq E_\\mathrm{PeriodicBC} \\leq E_\\mathrm{DirichletBC} \\leq\n", "E_\\mathrm{Voigt}$. A simple thing to compare are the eigenvalues of the tensors. Here\n", "we look at the first eigenvalue:" ], "metadata": {} }, { "outputs": [ { "output_type": "execute_result", "data": { "text/plain": "(2.05e10, 2.34e10, 2.82e10, 5.84e10)" }, "metadata": {}, "execution_count": 20 } ], "cell_type": "code", "source": [ "ev = (first ∘ eigvals).((E_reuss, E_periodic, E_dirichlet, E_voigt))\n", "round.(ev; digits=-8)" ], "metadata": {}, "execution_count": 20 }, { "cell_type": "markdown", "source": [ "Finally, we export the solution and the stress field to a VTK file. For the export we\n", "also compute the macroscopic part of the displacement." ], "metadata": {} }, { "outputs": [], "cell_type": "code", "source": [ "uM = zeros(ndofs(dh))\n", "\n", "vtk_grid(\"homogenization\", dh) do vtk\n", " for i in 1:3\n", " # Compute macroscopic solution\n", " apply_analytical!(uM, dh, :u, x -> εᴹ[i] ⋅ x)\n", " # Dirichlet\n", " vtk_point_data(vtk, dh, uM + u.dirichlet[i], \"_dirichlet_$i\")\n", " vtk_point_data(vtk, projector, σ.dirichlet[i], \"σvM_dirichlet_$i\")\n", " # Periodic\n", " vtk_point_data(vtk, dh, uM + u.periodic[i], \"_periodic_$i\")\n", " vtk_point_data(vtk, projector, σ.periodic[i], \"σvM_periodic_$i\")\n", " end\n", "end;" ], "metadata": {}, "execution_count": 21 }, { "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.8.5" }, "kernelspec": { "name": "julia-1.8", "display_name": "Julia 1.8.5", "language": "julia" } }, "nbformat": 4 }