.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/03-advanced/composite_DCB_test_3D.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_03-advanced_composite_DCB_test_3D.py: Composite double-cantilever beam with a cohesive interface ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The model consists of two composite arms separated by a zero-thickness cohesive interface. Cohesive elements are inserted only after the prescribed starter-crack length. The initially unconnected part of the two arms therefore represents the pre-crack. The arms are made from a unidirectional composite whose fibres are aligned with the beam axis. A mixed-mode bilinear :class:`fedoo.constitutivelaw.CohesiveLaw` governs the interface, and the opening displacement is applied incrementally with :class:`fedoo.problem.NonLinear`. The example uses N, mm, and MPa. Its compact mesh demonstrates the cohesive-law workflow; it is not intended as a mesh-converged fracture benchmark. .. GENERATED FROM PYTHON SOURCE LINES 19-26 .. code-block:: Python from __future__ import annotations import fedoo as fd import numpy as np .. GENERATED FROM PYTHON SOURCE LINES 27-32 Geometry and discretization ~~~~~~~~~~~~~~~~~~~~~~~~~~~ We first define the DCB dimensions, loading, and compact mesh resolution. The mesh is deliberately small so that the complete nonlinear example remains fast enough for the documentation gallery. .. GENERATED FROM PYTHON SOURCE LINES 32-47 .. code-block:: Python # Number of nodes along the beam length, width, and through each arm. nx = 19 ny = 5 nz_per_arm = 3 length = 60.0 width = 10.0 arm_thickness = 1.0 crack_length = 10.0 opening = 10 fd.ModelingSpace("3D") .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 48-53 Mesh the two arms ~~~~~~~~~~~~~~~~ The lower and upper meshes deliberately retain separate nodes at z=0. Coincident-but-distinct nodes are necessary because the cohesive element interpolates the displacement jump between the two interface faces. .. GENERATED FROM PYTHON SOURCE LINES 53-99 .. code-block:: Python lower = fd.mesh.box_mesh( nx, ny, nz_per_arm, x_min=0.0, x_max=length, y_min=0.0, y_max=width, z_min=-arm_thickness, z_max=0.0, elm_type="hex8", ) upper = fd.mesh.box_mesh( nx, ny, nz_per_arm, x_min=0.0, x_max=length, y_min=0.0, y_max=width, z_min=0.0, z_max=arm_thickness, elm_type="hex8", ) # Mesh.stack concatenates the two node lists without merging coincident nodes. domain = fd.Mesh.stack(lower, upper, name="DCB_domain") upper_offset = lower.n_nodes # Define one volume mesh per arm using the common global node list. Separate # meshes make it possible to assign distinct materials or ply orientations to # the arms without altering the interface node numbering. lower_volume = fd.Mesh( domain.nodes, lower.elements, "hex8", ndim=3, ) upper_volume = fd.Mesh( domain.nodes, upper.elements + upper_offset, "hex8", ndim=3, ) .. GENERATED FROM PYTHON SOURCE LINES 100-108 Construct the cohesive interface ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ For hex8 elements, nodes 4:8 form the upper face and nodes 0:4 form the lower face. Pair the top face of the lower arm with the bottom face of the upper arm. A quad4interface element therefore contains eight nodes: [four nodes on the lower arm, four nodes on the upper arm] .. GENERATED FROM PYTHON SOURCE LINES 108-128 .. code-block:: Python n_interface_cells = (nx - 1) * (ny - 1) lower_interface_faces = lower.elements[-n_interface_cells:, 4:8] upper_interface_faces = upper.elements[:n_interface_cells, 0:4] + upper_offset interface_elements = np.hstack((lower_interface_faces, upper_interface_faces)) # Do not insert cohesive elements over the starter-crack region. The two arms # are consequently disconnected for x < crack_length and bonded by cohesive # elements for x >= crack_length. interface_centers_x = domain.nodes[interface_elements].mean(axis=1)[:, 0] interface_elements = interface_elements[interface_centers_x >= crack_length - 1.0e-12] interface = fd.Mesh( domain.nodes, interface_elements, "quad4interface", ndim=3, name="DCB_cohesive_interface", ) .. GENERATED FROM PYTHON SOURCE LINES 129-133 Constitutive laws and assemblies ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Both unidirectional composite arms have fibres along the beam axis (X). Separate objects are retained so either arm can easily be changed later. .. GENERATED FROM PYTHON SOURCE LINES 133-183 .. code-block:: Python lower_material = fd.constitutivelaw.CompositeUD( Vf=0.6, E_f=250000.0, E_m=3500.0, nu_f=0.33, nu_m=0.3, angle=0.0, ) upper_material = fd.constitutivelaw.CompositeUD( Vf=0.6, E_f=250000.0, E_m=3500.0, nu_f=0.33, nu_m=0.3, angle=0.0, ) # Bilinear mixed-mode cohesive law. Axis 2 is the local interface-normal # direction; the quad4interface element supplies the corresponding local frame. cohesive_law = fd.constitutivelaw.CohesiveLaw( GIc=0.3, SImax=60.0, KI=1.0e4, GIIc=1.6, SIImax=None, KII=1.0e4, tangent_mode="secant", ) lower_assembly = fd.Assembly.create( fd.weakform.StressEquilibrium(lower_material), lower_volume, ) upper_assembly = fd.Assembly.create( fd.weakform.StressEquilibrium(upper_material), upper_volume, ) cohesive_assembly = fd.Assembly.create( fd.weakform.InterfaceForce(cohesive_law), interface, ) # Keep a separate sum of the two volume assemblies for post-processing. It # contains only standard hex8 meshes and can therefore be converted directly # to a PyVista ``MultiMesh``. The cohesive assembly is then added only to the # global mechanical assembly. volume_assembly = lower_assembly + upper_assembly assembly = volume_assembly + cohesive_assembly .. GENERATED FROM PYTHON SOURCE LINES 184-188 Problem and boundary conditions ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Material softening makes the problem nonlinear even though large-deformation kinematics are disabled. .. GENERATED FROM PYTHON SOURCE LINES 188-206 .. code-block:: Python problem = fd.problem.NonLinear(assembly, nlgeom=False) # Load the end faces at x=0. The lower end is fully fixed; the upper end is # prevented from moving in X and Y and receives the prescribed opening in Z. lower_left = np.nonzero( np.isclose(domain.nodes[:, 0], 0.0) & (domain.nodes[:, 2] < -arm_thickness + 1.0e-12) )[0] upper_left = np.nonzero( np.isclose(domain.nodes[:, 0], 0.0) & (domain.nodes[:, 2] > arm_thickness - 1.0e-12) )[0] problem.bc.add("Dirichlet", lower_left, "Disp", 0.0) problem.bc.add("Dirichlet", upper_left, "DispX", 0.0) problem.bc.add("Dirichlet", upper_left, "DispY", 0.0) problem.bc.add("Dirichlet", upper_left, "DispZ", opening) .. rst-class:: sphx-glr-script-out .. code-block:: none Dirichlet boundary condition: var = 'DispZ' n_nodes = 5 value = 10 .. GENERATED FROM PYTHON SOURCE LINES 207-214 Incremental nonlinear solution ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The imposed opening is scaled by the problem time from zero to its final value. Adaptive increments help the Newton solver cross the damage events. Convergence may be slow for this kind of problem. During damage propagation, the residual may temporarily increase over several iterations before convergence, so the early-divergence check is disabled. .. GENERATED FROM PYTHON SOURCE LINES 214-225 .. code-block:: Python problem.set_nr_criterion(check_early_divergence=False) problem.nlsolve( dt=0.05, tmax=1.0, update_dt=True, dt_min=1.0e-5, tol_nr=5.0e-3, print_info=0, ) .. GENERATED FROM PYTHON SOURCE LINES 226-230 Cohesive-zone response ~~~~~~~~~~~~~~~~~~~~~~ Damage and relative displacement are stored at the interface Gauss points. A damage value of one denotes a fully failed cohesive point. .. GENERATED FROM PYTHON SOURCE LINES 230-244 .. code-block:: Python damage = np.asarray(cohesive_assembly.sv["DamageVariable"]) relative_disp = np.asarray(cohesive_assembly.sv["RelativeDisp"]) print("\nDCB cohesive-zone summary") print("-------------------------") print(f"n_nodes: {domain.n_nodes}") print("n_volume_elements: " f"{lower_volume.n_elements + upper_volume.n_elements}") print(f"n_cohesive_elements: {interface.n_elements}") print(f"max_damage: {damage.max()}") print(f"damaged_gauss_points: {np.count_nonzero(damage > 0.0)}") print("failed_gauss_points: " f"{np.count_nonzero(damage >= 1.0 - 1.0e-10)}") print(f"max_opening: {relative_disp[2].max()}") .. rst-class:: sphx-glr-script-out .. code-block:: none DCB cohesive-zone summary ------------------------- n_nodes: 570 n_volume_elements: 288 n_cohesive_elements: 60 max_damage: 1.0 damaged_gauss_points: 112 failed_gauss_points: 112 max_opening: 5.589760930506991 .. GENERATED FROM PYTHON SOURCE LINES 245-252 Plot the deformed DCB ~~~~~~~~~~~~~~~~~~~~~ Request results from the volume sum explicitly. Calling problem.get_results(...) without an assembly would select the summed assembly, which also contains the eight-node ``quad4interface`` elements. These elements are valid for the mechanical calculation but have no direct VTK/PyVista cell equivalent. .. GENERATED FROM PYTHON SOURCE LINES 252-275 .. code-block:: Python results = problem.get_results(volume_assembly, ["Stress", "Strain", "Disp"]) # The AssemblySum result is a single MultiMesh dataset containing both arms. # Element fields are retained independently on each submesh, whereas a # recovered nodal stress field would be treated as one shared global field and # would keep only the first arm's recovery. Displacement remains a shared nodal # field and is used to display the final deformed geometry. plotter = results.plot( "Stress", component="vm", data_type="Node", show_edges=True, show=False, title="DCB - von Mises stress", scalar_bar_args={"interactive": False}, ) # Use an isometric camera from the negative-X side. This reverses the visual # direction of the beam axis without changing the mesh, crack position, or # boundary conditions. plotter.view_vector((-1, 1, 1), viewup=(0, 0, 1)) plotter.reset_camera() plotter.show() .. image-sg:: /examples/03-advanced/images/sphx_glr_composite_DCB_test_3D_001.png :alt: composite DCB test 3D :srcset: /examples/03-advanced/images/sphx_glr_composite_DCB_test_3D_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.228 seconds) .. _sphx_glr_download_examples_03-advanced_composite_DCB_test_3D.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: composite_DCB_test_3D.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: composite_DCB_test_3D.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: composite_DCB_test_3D.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_