Search for the crack path ================================== Generalities ----------- It is assumed that the crack belongs to a certain group of elements (GROUP_MA), on which field :math:`X` (damage, plastic deformation) is defined. The crack is identified by an area where the values in this field are especially high, and normally greater than a certain value: damage greater than :math:`0` for example. In the following, the procedure for finding the crack path is described. We exploit the fact that, in the plane orthogonal to the crack, field :math:`X` has a maximum on the crack itself. The procedure is step-by-step, with each step a new crack point is found. It can be described by the *not-type*, by its *initiation* and by the *stopping criteria*. Not typical -------- At each stage, research requires two points. The last point found :math:`{P}_{i}` is the *starting point* and the vector :math:`\overrightarrow{{P}_{i\mathrm{-}1}{P}_{i}}` (:math:`{P}_{i\mathrm{-}1}` being the penultimate point found) gives the *search direction*. .. image:: images/10000000000002EA000001A34D6CDC1A76A17F35.png :width: 4.0909in :height: 2.289in .. _RefImage_10000000000002EA000001A34D6CDC1A76A17F35.png: Figure 1: Diagram of the pass type At each typical step, the following actions are executed (see also the diagram of): 1. the position of the next point is first estimated (*prediction point* :math:`{P}_{i}^{\mathit{pr}}`) at the distance :math:`a` from the starting point in the search direction, :math:`a` being the *no progress;* 2. field :math:`X` is projected onto a line orthogonal to the search direction and passing through the prediction point; .. _RefNumPara__1649_1916982135: 3. the projected field is smoothed by convolution (**eq**) to partially get rid of the mesh (for example, if we rely on finite elements with linear interpolation, the projected field is linear in pieces); 4. the new point in crack :math:`{P}_{i+1}` is where the smoothed field reaches its maximum value. The smoothing mentioned in point :ref:`3 ` is necessary: for example, if the field projected onto the orthogonal profile is linear in pieces, its maximum is necessarily on the edge of an element. The smoothing is achieved by convolution: :math:`\overline{X}(s)\mathrm{=}\frac{{}_{{l}_{\mathit{orth}}}\mathrm{\int }X(\zeta )\Psi (\mathrm{\mid }\zeta \mathrm{-}s\mathrm{\mid })d\zeta }{{}_{{l}_{\mathit{orth}}}\mathrm{\int }\Psi (\mathrm{\mid }\zeta \mathrm{-}s\mathrm{\mid })d\zeta }` **eq** 2.2-1 where :math:`\overline{X}` is the smoothed field, :math:`{l}_{\mathit{orth}}` is the length of the line orthogonal to the search direction, :math:`s` is the curvilinear abscissa on this line, and :math:`\Psi` is a Gauss function: :math:`\Psi (\mathrm{\mid }\zeta \mathrm{-}s\mathrm{\mid })\mathrm{=}\mathrm{exp}(\mathrm{-}{(\frac{2(\zeta \mathrm{-}s)}{{l}_{\mathit{reg}}})}^{2})` **eq** 2.2-2 An example of fields :math:`X` and :math:`\overline{X}` based on the parameter :math:`{l}_{\mathit{reg}}` is given in. .. image: images/1000000000000242000001 BB53FDCD79BE30BF60 .tif :width: 3.789 in :height: 2.9035 in .. _refImage_1000000000000242000001 BB53FDCD79BE30BF60 .tif: Figure 2: Smoothing field :math:`X` projected onto the orthogonal line. An additional parameter is given by the number of points on the orthogonal line, :math:`{N}_{\mathit{orth}}`. This determines the precision of the discrete convolution product, as well as the :math:`{\delta }_{\mathit{orth}}\mathrm{=}{l}_{\mathit{orth}}\mathrm{/}{N}_{\mathit{orth}}` precision (distance between points). Discrete convolution is performed by the *convolve* function from the *numpy* python library. The functions s :math:`X(s)` and :math:`\Psi (s)` are sampled over :math:`{l}_{\mathit{orth}}` at the constant distance :math:`{\delta }_{\mathit{orth}}`. By calling :math:`\mathrm{X}`, :math:`\Psi` the two sampled vectors, the point :math:`\overline{{X}_{j}}` of the vector :math:`\overline{\mathrm{X}}` that approximates the smoothed function :math:`\overline{X}(s)` by points will be given by: :math:`\overline{{X}_{j}}\mathrm{=}\frac{\mathrm{\sum }_{i}{\Psi }_{i}{X}_{i+j}}{\mathrm{\sum }_{i}{\Psi }_{i}}` **eq** 2.2-3 Initialization -------------- The purpose of initialization is to get the first two points, to continue after as described in the pas-type. A diagram of the initialization procedure is given in: .. image:: images/1000000000000245000001CFE77C4B917749C01A.png :width: 3.2598in :height: 2.6346in .. _RefImage_1000000000000245000001CFE77C4B917749C01A.png: Figure 3: diagram of the initialization procedure. In detail: 1. field :math:`X` is projected onto a circle with radius :math:`a` and centers on node :math:`1` where the field at the nodes has its maximum on the chosen GROUP_MA; 2. field :math:`X` is smoothed over the circle via a discrete convolution product (**eq.**), the second point of the crack :math:`{P}_{2}` is the point in the circle where :math:`\overline{X}` is maximum; 3. the position of the point :math:`1` is corrected by projecting and smoothing the field :math:`\overline{X}` on a line orthogonal to :math:`\overrightarrow{{P}_{2}1}` and going through :math:`{P}_{2}`: the new point :math:`3={P}_{1}` is found; 4. two possible research directions are thus determined: :math:`\overrightarrow{{P}_{1}{P}_{2}}` and :math:`\overrightarrow{{P}_{2}{P}_{1}}`. Stop criteria ----------------- The search for the crack path is completed in both directions :math:`\overrightarrow{{P}_{1}{P}_{2}}` and :math:`\overrightarrow{{P}_{2}{P}_{1}}` until stopped, in one of the following cases: * the value of the field corresponding to the new point is less than a threshold value, defined a priory. This value must be entered among the command parameters under the BORNE_MIN operand of the RECHERCHE factor keyword. * the prediction point is outside the subject. Fields for the search for the crack path ------------------------------------------------- The method described can be applied to all scalar fields that represent the cracking of the material. Typically, damage :math:`d` can be used. Sometimes the damage field is too "flat", and it breaks a plateau with :math:`d\mathrm{=}1`. As there would then not be a well-defined maximum, the procedure described in paragraphs 2.2 - 2.3 may not give the desired results. It is then possible to choose other relevant fields: a main deformation component, an equivalent deformation... An example of choosing the relevant field is given by the test cases zzzz264b, c [:external:ref:`V1.01.264 `].