{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Continuum mechanics -- crack tip stress\n", "\n", "In this lesson we will look at the stress distribution around a crack tip. Here we switch from atomistic simulations to continuum mechanics with the finite element method (FEM) using the [FEniCS code](https://fenicsproject.org) under the hood. \n", "\n", "We won't spend very much time on the mathematics for FEM, but will examine some of the common input parameters and their effect.\n", "\n", "For physics, we will examine three fracture modes for a simple triangular crack in a cubic sample. This is shown schematically below courtesy of wikipedia for a rectangular prism:\n", "\n", "![wikipedia fracture mechanics](https://upload.wikimedia.org/wikipedia/commons/thumb/e/e7/Fracture_modes_v2.svg/2560px-Fracture_modes_v2.svg.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Setup" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:54:51.111449Z", "start_time": "2020-12-26T13:54:46.639369Z" } }, "outputs": [], "source": [ "%matplotlib notebook\n", "from pyiron_base import Project\n", "import pyiron_continuum\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:54:52.166156Z", "start_time": "2020-12-26T13:54:51.115069Z" } }, "outputs": [], "source": [ "pr = Project('fenics_linear_elasticity')\n", "pr.remove_jobs_silently(recursive=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# FEM basics\n", "\n", "FEM is an approach for solving partial differential equations (PDEs) numerically by discretizing them with a mesh and using the calculus of variations. \n", "\n", "The key physical concepts which specify the problem are the equation being solved (obviously), the spatial domain on which this equation is being solved, and the boundary conditions on that domain. The last of these is typically broken down into two varieties: \"Dirichlet\" boundary conditions, which specify the value of the field being solved for, and \"Neumann\" boundary conditions, which specify the gradient of that field dotted with the domain normal vector.\n", "\n", "Some of the important numerical concepts are the mesh generation (we will look in a bit of detail at the effect of mesh density) and the element order, which controls how the fields are interpolated within the individual elements (which we will touch on only briefly).\n", "\n", "There is much, much more depth on both the mathematical and numeric sides of FEM, but it is not my core area of research and from here on we will restrict ourselves to a concrete example." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear Elasticity\n", "\n", "The static solution to linear elasticity can be summarized with three equations:\n", "\n", "1. $\\nabla \\cdot \\sigma = -f$\n", "2. $\\sigma = \\lambda~\\mathrm{tr}(\\epsilon) \\mathrm{I} + 2 \\mu \\epsilon$\n", "3. $\\epsilon = \\frac{1}{2}\\left(\\nabla u + (\\nabla u)^\\mathrm{T} \\right)$\n", "\n", "Where the derivatives of the stress tensor, $\\sigma$ balance out body forces on the sample, $f$. The stress tensor is shown here using Lame's constants $\\lambda$ and $\\mu$ -- which we can also describe in terms of very familair bulk ($K$) and shear ($G$) moduli as $\\lambda = K - (2 G / 3)$ and $\\mu = G$ -- and the symmeterized strain rate tensor, $\\epsilon$, which in turn is constructed from gradients of the displacement field $u$.\n", "\n", "With FEM we convert this the so-called weak form which contains integrals and a test function `v`. In this form, we solve\n", "\n", "$\\int_\\Omega \\sigma(u):\\epsilon(v) dx = \\int_\\Omega f \\cdot v dx + \\int_{\\partial \\Omega_T} T \\cdot v ds$\n", "\n", "Where $\\Omega$ is the domain of our sample and $\\partial \\Omega_T$ is its boundaries subject to the traction $T$ (remaining boundaries are subject to Dirichlet conditions where the solution is given directly), and $:$ is a tensor contraction. \n", "\n", "For the problems we look at here both $f$ and $T$ will be strictly zero, i.e. we have a linear elastic material subject to boundaries which are either displacement controlled or free to relax with no external forces on them. So the physcis of our solution (displacement and strain) is controlled by our material properties ($K$ and $G$ moduli), the geometry of our sample, and how we control its deformation with Dirichlet boundary conditions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Geometry\n", "\n", "We'll create our sample by making a box, then subtracting out a triangular crack from this domain. The code to do this is nicely wrapped up in the function below, `set_domain_to_cracked_box`. The important thing for us is that we can easily play around with the crack length, width, and depth (in the x-, y-, z-directions of our box, respectively)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def set_domain_to_cracked_box(job, crack_length=1.0, crack_width=0.1, crack_depth=0.5):\n", " \"\"\"\n", " Sets the job's domain to a unit box with a triangular crack tip in it. The crack starts at the\n", " z=0 plane running along the x-axis.\n", "\n", " The domain is initially a box, then we subtract out a triangular wedge made from three tets.\n", " The mesh generation is not perfect and there can be small defects right at the crack tip.\n", "\n", " Args:\n", " job (Fenics): The job whose domain to set.\n", " crack_length (float): How long the crack is. (Default is 1, run the entire length of the box.)\n", " crack_width (float): How wide the mouth of the crack is. (Default is 0.1.)\n", " crack_depth (float): How deep the crack is. (Default is 0.5, reach to the center of the box.\n", " \"\"\"\n", " bulk = job.create.domain.box()\n", " p1 = (0.5 * (1 - crack_length), 0.5 * (1 - crack_width), 0.)\n", " p2 = (0.5 * (1 - crack_length), 0.5 * (1 + crack_width), 0.)\n", " p3 = (0.5 * (1 + crack_length), 0.5 * (1 + crack_width), 0.)\n", " p4 = (0.5 * (1 + crack_length), 0.5 * (1 - crack_width), 0.)\n", " p5 = (0.5 * (1 - crack_length), 0.5, crack_depth)\n", " p6 = (0.5 * (1 + crack_length), 0.5, crack_depth)\n", " crack1 = job.create.domain.tetrahedron(p1, p2, p3, p5)\n", " # crack2 = job.create.domain.tetrahedron(p1, p2, p3, p6)\n", " # crack3 = job.create.domain.tetrahedron(p1, p2, p4, p5)\n", " # crack4 = job.create.domain.tetrahedron(p1, p2, p4, p6)\n", " # crack5 = job.create.domain.tetrahedron(p1, p2, p5, p6)\n", " # crack6 = job.create.domain.tetrahedron(p1, p3, p4, p5)\n", " crack7 = job.create.domain.tetrahedron(p1, p3, p4, p6)\n", " crack8 = job.create.domain.tetrahedron(p1, p3, p5, p6)\n", " # crack9 = job.create.domain.tetrahedron(p2, p3, p4, p5)\n", " # crack10 = job.create.domain.tetrahedron(p2, p3, p4, p6)\n", " # crack11 = job.create.domain.tetrahedron(p3, p4, p5, p6)\n", " job.domain = bulk - crack1 - crack7 - crack8" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's create a new FEniCS FEM job and examine the behaviour of this function. In addition to the geometry, we will also increase the mesh density from its very low default value of `2`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:54:55.488036Z", "start_time": "2020-12-26T13:54:52.168872Z" } }, "outputs": [], "source": [ "job = pr.create.job.FenicsLinearElastic('fem')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:54:55.497768Z", "start_time": "2020-12-26T13:54:55.492227Z" } }, "outputs": [], "source": [ "set_domain_to_cracked_box(job)\n", "job.input.mesh_resolution = 10" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's visualize the mesh" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T18:56:08.720763Z", "start_time": "2020-12-26T18:56:08.705230Z" } }, "outputs": [], "source": [ "job.mesh" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 1\n", "\n", "Play around with the `crack_length`, `crack_width`, and `crack_depth` arguments until you comfortably understand the effect they have on the sample geometry. Then increase the `input.mesh_resolution` and note its impact on both the overall meshing, and especially on the numeric artefacts right at the tip.\n", "\n", "Note: After calling `set_domain_to_cracked_box` again and/or changing the `job.input.mesh_resolution`, you will need to call `job.generate_mesh()` to update the mesh to use the new parameters." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Solving for displacement and stress\n", "\n", "Let's solve the displacement and von Mises stress for mode 1: openening. This calculation will serve as an example of how the syntax works for the other calculations we'll perform." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running the calculation\n", "\n", "First, we'll instantiate the job and set the material properties, for which we'll use experimental values for Al reported [on wikipedia](https://en.wikipedia.org/wiki/Aluminium) in GPa." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:54:58.515123Z", "start_time": "2020-12-26T13:54:58.496606Z" } }, "outputs": [], "source": [ "mode1 = pr.create.job.FenicsLinearElastic('mode1', delete_existing_job=True)\n", "K_Al, G_Al = 76, 26\n", "mode1.input.bulk_modulus = K_Al\n", "mode1.input.shear_modulus = G_Al" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we'll apply a technical setting: switch the solver from its default value of direct solution to an iterative solver. This allows us to increase the mesh density to larger values without running out of memory. If you want to know more details about FEM solvers, you can read about them [here](https://www.simscale.com/blog/2016/08/how-to-choose-solvers-for-fem/). For our purposes it's sufficient to know that we're switching over to an iterative solver." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:54:58.567904Z", "start_time": "2020-12-26T13:54:58.520105Z" } }, "outputs": [], "source": [ "mode1.input.solver_parameters = {\n", " 'linear_solver': 'gmres',\n", " 'preconditioner': 'ilu'\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After this, we'll construct our geometry. We saw above how the constructor works, but we also want to be able to find nodes that belong to the boundaries so we can apply boundary conditions." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:54:58.587479Z", "start_time": "2020-12-26T13:54:58.579909Z" } }, "outputs": [], "source": [ "mode1.input.mesh_resolution = 30\n", "set_domain_to_cracked_box(mode1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "FEniCS, which pyiron is using under the hood, does this with a special python function that takes two arguments: the position of the node, `x`, which is three-dimensional for us, and a boolean array `on_boundary` which FEniCS uses internally to keep track of which nodes are on *any* boundary.\n", "\n", "For our boundary conditions, we'll hold the end of the sample opposite the crack fixed, and then displace the face of the sample where the crack starts according to which deformation mode we want. For this, we'll need functions to find the face where the crack *starts* and *ends*, as well as functions to see if we have the *top* half of the face above the crack or the *bottom* half of the face below it.\n", "\n", "Using `set_domain_to_cracked_box`, our crack penetrates along the z-axis runs the length of the x-axis, i.e. whether we're above or below it is determined by our y-position. Thus, our logical conditions will use `x[2]` (z) and `x[1]` (y)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:54:58.600894Z", "start_time": "2020-12-26T13:54:58.592612Z" } }, "outputs": [], "source": [ "def near_start(x, on_boundary):\n", " return job.fenics.near(x[2], 0.)\n", "\n", "def near_end(x, on_boundary):\n", " return job.fenics.near(x[2], 1.)\n", "\n", "def top_half(x, on_boundary):\n", " return on_boundary and x[1] > 0.5 and near_start(x, on_boundary)\n", "\n", "def bottom_half(x, on_boundary): \n", " return on_boundary and x[1] < 0.5 and near_start(x, on_boundary)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our actual boundary conditions will be \"Dirichlet\" (i.e. displacement) conditions on these faces. To communicate these to FEniCS we use special data types called `Constant` (for things that don't change, obviously) and `Expression` (in case we have something with variables). Here let's just apply a static strain using `Constant`.\n", "\n", "Activating mode 1, we'll move the top half up a bit and the bottom half down a bit in the y-direction:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:00.010269Z", "start_time": "2020-12-26T13:54:58.610148Z" } }, "outputs": [], "source": [ "strain = 0.01\n", "top_bc = mode1.create.bc.dirichlet(mode1.Constant((0, strain, 0)), top_half)\n", "bottom_bc = mode1.create.bc.dirichlet(mode1.Constant((0, -strain, 0)), bottom_half)\n", "rigid_bc = mode1.create.bc.dirichlet(mode1.Constant((0, 0, 0)), near_end)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These boundary conditions (BCs) get applied to the job as a list." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:00.037732Z", "start_time": "2020-12-26T13:55:00.014450Z" } }, "outputs": [], "source": [ "mode1.BC = [rigid_bc, top_bc, bottom_bc]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now all that's left to do is run our job.\n", "\n", "Note: FEniCS funcitonality in pyiron is still experimental, so unlike other jobs which get automatically saved and can be re-loaded later, these jobs exist only in the notebook. Thankfully for us, this problem does not use very much CPU time, so we can always simply re-run the calculations without too much headache." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:01.354113Z", "start_time": "2020-12-26T13:55:00.043735Z" } }, "outputs": [], "source": [ "mode1.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Analyzing run\n", "\n", "The main output from this run is the displacement field (i.e. the `solution`) and the von Mises stress, which is a scalar field value we post-process from the solution and indicates where plastic activity is most likely to start: $\\sigma_M = \\sqrt{\\frac{3}{2}s : s}$, where $s = \\sigma - \\frac{1}{3}\\mathrm{tr}(\\sigma)\\mathrm{I}$ is the deviatoric strain.\n", "\n", "Let's start by looking at a 3D plot of this stress:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:01.516970Z", "start_time": "2020-12-26T13:55:01.357723Z" } }, "outputs": [], "source": [ "mode1.plot.stress3d();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reassuringly, there is a line of high-stress dots along the very tip of the crack, as we would expect. Otherwise most of the stress seems to be focused in the half of the sample containing the crack.\n", "\n", "3D plots are useful, but sometimes it is more effective to look at some 2D projection. We also have helper functions to do this, which project all the results onto a single plane. Since we'd like to isolate the crack tip, it's helpful to project on the x-axis onto the yz-plane so that we're looking down the length of the crack." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:02.297549Z", "start_time": "2020-12-26T13:55:01.549081Z" } }, "outputs": [], "source": [ "mode1.plot.stress2d(projection_axis=0);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Oh no! Where did our crack go?\n", "\n", "Actually, it is still there -- this is just a simple interpolating colourmap, so along the crack where we have no data the plot is simply interpolating. This is actually visible: if you look carfully you will see that there is a triangle with its base on the horizontal axis that has horizontal striated colour bands where the interpolation is happening from one side of the open crack to the other.\n", "\n", "As with the 3D plot, but a bit easier to see now, the stress is indeed focused at the tip of the crack. It's so focused that this bright spot actually washes out information about the rest of the stress distribution. We can get around this by using a logarithmic scaling for our colours:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:02.999338Z", "start_time": "2020-12-26T13:55:02.303327Z" } }, "outputs": [], "source": [ "mode1.plot.stress2d(projection_axis=0, lognorm=True);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we might also be interested in the raw numeric values -- like what is the maximum nodal stress value? These nodal values are stored in the output:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:03.011435Z", "start_time": "2020-12-26T13:55:03.002546Z" } }, "outputs": [], "source": [ "print(\"Max nodal stress = {}\".format(mode1.output.von_Mises[-1].max()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Thinking Deeper\n", "\n", "Which mode do you expect will have the largest maximum stress, which the smallest, and why?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your thoughts here*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise\n", "\n", "Let's get some data to test your hypothesis! Run the same calculation for the other two modes; examine the stress plots and directly compare the maximum stresses. Does the data support your hypothesis, or illuminate something new?\n", "\n", "Hint: We only need to change the displacements in our Dirichlet boundary conditions." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:03.033168Z", "start_time": "2020-12-26T13:55:03.027896Z" } }, "outputs": [], "source": [ "# Your code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Quasistatic strain profiles\n", "\n", "You may have noticed that just like atomistic output we looked at the output with a `[-1]`, i.e. `job.output.von_Mises[-1]`. That's because just like the atomistic data we can look at a time series of results.\n", "\n", "In this case, it's fairly easy to set up an experiment studying the stress as a function of increasing strain in the quasistatic limit (i.e. ignoring any sort of momentum effects, etc.). The key difference will be that our boundary conditions are now a time (`t`) dependent `Expression` instead of a `Constant`. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:09.841205Z", "start_time": "2020-12-26T13:55:09.793735Z" } }, "outputs": [], "source": [ "mode1t = pr.create.job.FenicsLinearElastic('mode1t', delete_existing_job=True)\n", "mode1t.input.bulk_modulus = K_Al\n", "mode1t.input.shear_modulus = G_Al\n", "mode1t.input.solver_parameters = {\n", " 'linear_solver': 'gmres',\n", " 'preconditioner': 'ilu'\n", "}\n", "mode1t.input.mesh_resolution = 30\n", "set_domain_to_cracked_box(mode1t)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:10.942435Z", "start_time": "2020-12-26T13:55:09.845919Z" } }, "outputs": [], "source": [ "strain_step = 0.005\n", "dirichlet_top = mode1t.Expression(('0', 'a * t', '0'), degree=2, a=strain_step, t=0)\n", "dirichlet_bot = mode1t.Expression(('0', '-a * t', '0'), degree=2, a=strain_step, t=0)\n", "\n", "top_bc = mode1t.create.bc.dirichlet(dirichlet_top, top_half)\n", "bottom_bc = mode1t.create.bc.dirichlet(dirichlet_bot, bottom_half)\n", "rigid_bc = mode1t.create.bc.dirichlet(mode1t.Constant((0, 0, 0)), near_end)\n", "mode1t.BC = [rigid_bc, top_bc, bottom_bc]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lastly, we need to let the job know which expressions are time-dependent so it can update their `t` parameter, and tell it how many steps to run for:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:10.954165Z", "start_time": "2020-12-26T13:55:10.945759Z" } }, "outputs": [], "source": [ "mode1t.time_dependent_expressions.append(dirichlet_top)\n", "mode1t.time_dependent_expressions.append(dirichlet_bot)\n", "mode1t.input.n_steps = 10" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:19.837939Z", "start_time": "2020-12-26T13:55:10.957485Z" } }, "outputs": [], "source": [ "mode1t.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As before we can plot the output, but we can also choose which timestep to look at by setting the `frame` argument. Note how the top end of the stress plot is narrower and narrower compared to the bottom part as we apply more and more strain!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:20.581880Z", "start_time": "2020-12-26T13:55:19.949216Z" } }, "outputs": [], "source": [ "mode1t.plot.stress2d(frame=-1, projection_axis=0, lognorm=True);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also examine the peak stress as a function of strain:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:20.666418Z", "start_time": "2020-12-26T13:55:20.592951Z" } }, "outputs": [], "source": [ "fig, ax = plt.subplots() # A little bit of overhead because we're using interactive plots in this notebook\n", "peak_stress_mode1t = np.array(mode1t.output.von_Mises).max(axis=1)\n", "ax.scatter(strain_step * np.arange(len(peak_stress_mode1t)), peak_stress_mode1t)\n", "ax.set_xlabel('Strain')\n", "ax.set_ylabel('Peak stress');" ] }, { "cell_type": "markdown", "metadata": { "ExecuteTime": { "end_time": "2020-12-26T09:47:02.728725Z", "start_time": "2020-12-26T09:47:02.721920Z" } }, "source": [ "## Exercise\n", "\n", "What do the stress strain profiles look like for the other two deformation modes?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:20.678144Z", "start_time": "2020-12-26T13:55:20.669547Z" } }, "outputs": [], "source": [ "# Your code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Numerics\n", "\n", "At the heart of FEM is a spatial discretization of our PDE, thus the mesh density for this discretization is a critical feature. Let's examine how our solution changes as a function of our numeric paramters.\n", "\n", "Since changing the mesh density means the number and position of the nodes, we can't simply compare the output nodal positions. However, the stored solution at the end of the run can be evaluated at points other than just the mesh nodes by using the function defined below like `evaluate_solution(job.solution)`" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:41.393370Z", "start_time": "2020-12-26T13:55:41.373284Z" } }, "outputs": [], "source": [ "spacing = np.linspace(0, 1, num=50)\n", "x, y, z = np.meshgrid(spacing, spacing, spacing)\n", "grid_points = np.vstack((x.flatten(), y.flatten(), z.flatten())).T\n", "\n", "def evaluate_solution(solution, points=grid_points):\n", " values = []\n", " for p in points:\n", " try:\n", " values.append(solution(p))\n", " except RuntimeError:\n", " # There is no solution where we cut out the crack\n", " pass\n", " return np.array(values)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise\n", "\n", "For a variety of mesh values, e.g. `[20, 30, 40, 50, 60]`, look at the root mean square difference of the solution on our uniform grid -- how does this converge as a function of the `input.mesh_resolution`? Use only your favourite deformation mode for this exercise, and a strain of 1%.\n", "\n", "Note: we can't exploit the time-dependent expression here to do multiple mesh densities in a single calculation. Just make a regular for-loop and run multiple jobs.\n", "\n", "Hint: You can get the magnitude of the difference between two arrays (of the same shape) very quickly using `np.linalg.norm(array1 - array2)`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2020-12-26T13:55:41.417584Z", "start_time": "2020-12-26T13:55:41.404574Z" } }, "outputs": [], "source": [ "# Your code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise\n", "\n", "So far we have been using the default value for `input.element_order` -- 1. That means that values are linearly interpolated throughout each individual finite element formed by our mesh.\n", "\n", "Increase this value from 1 to 2, so that values are quadratically interpolated, and repeat the above exercise." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Your code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Thinking deeper\n", "\n", "Consider the relationship between mesh resolution and element order, their relative accuracy, and expense? What guidelines can we think of for when to use each setting?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Physics\n", "\n", "Finally, let's experiment with a couple of the actual material system parameters." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise\n", "\n", "Increase the crack width smoothly from 0.1 of the sample width to 0.25. What is the effect on the maximum stress value? Does this agree with our textbook knowledge?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Your code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise\n", "\n", "Suppose that Al had the same bulk modulus, but double the shear modulus. What effect do you expect this to have on the maximum stresses across the different deformation modes? Repeat the quasistatic strain exercise for all three modes with `job.input.shear_modulus = 2 * G_Al` and compare these results to the original results. We'll keep $\\lambda$ the same, as given below by modifying the bulk modulus." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Your code here" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:pyiron_38]", "language": "python", "name": "conda-env-pyiron_38-py" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.6" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }