Fix/example00 block code

examples/00-mapdl-examples/3d_notch.py

@@ -12,8 +12,8 @@
 First, start MAPDL as a service and disable all but error messages.
 """
 # sphinx_gallery_thumbnail_number = 3
-
 import numpy as np
+from matplotlib import pyplot as plt
 from ansys.mapdl.core import launch_mapdl
 
 mapdl = launch_mapdl(loglevel="ERROR")
@@ -27,12 +27,7 @@
 length = 0.4
 width = 0.1
 
-# ratio = 0.3  # diameter/width
-# diameter = width*ratio
-# radius = diameter*0.5
-
 notch_depth = 0.04
-# notch_radius = 0.002
 notch_radius = 0.01
 
 # create the half arcs
@@ -61,7 +56,6 @@
 
 rect_anum = mapdl.blc4(width=length, height=width)
 
-
 # Note how pyansys parses the output and returns the area numbers
 # created by each command.  This can be used to execute a boolean
 # operation on these areas to cut the circle out of the rectangle.
@@ -73,17 +67,16 @@
 mapdl.lplot(vtk=True, show_keypoint_numbering=True)
 mapdl.lsel("all")
 
-
 # plot the area using vtk/pyvista
 mapdl.aplot(vtk=True, show_area_numbering=True, show_lines=True, cpos="xy")
 
-
-# ###############################################################################
 # Next, extrude the area to create volume
 thickness = 0.01
 mapdl.vext(cut_area, dz=thickness)
 
-mapdl.vplot(vtk=True, show_lines=True, show_axes=True, smooth_shading=True)
+# Checking volume plot
+_ = mapdl.vplot(vtk=True, show_lines=True, show_axes=True, smooth_shading=True)
+
 
 ###############################################################################
 # Meshing
@@ -100,7 +93,6 @@
 
 # define a PLANE183 element type with thickness
 
-
 # ensure there are at 25 elements around the hole
 notch_esize = np.pi * notch_radius * 2 / 50
 plate_esize = 0.01
@@ -181,7 +173,7 @@
 mapdl.f("ALL", "FX", 1000)
 
 # finally, be sure to select all nodes again to solve the entire solution
-_ = mapdl.allsel()
+_ = mapdl.allsel(mute=True)
 
 
 ###############################################################################
@@ -191,7 +183,8 @@
 mapdl.run("/SOLU")
 mapdl.antype("STATIC")
 mapdl.solve()
-mapdl.finish()
+_ = mapdl.finish()
+
 
 ###############################################################################
 # Post-Processing
@@ -219,6 +212,7 @@
 # Must use nanmax as stress is not computed at mid-side nodes
 max_stress = np.nanmax(von_mises)
 
+
 ###############################################################################
 # Compute the Stress Concentration
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -244,7 +238,6 @@
 print("Far field von mises stress: %e" % far_field_stress)
 # Which almost exactly equals the analytical value of 10000000.0 Pa
 
-# result.plot_element_result(0, 'ENS', 0)
 
 ###############################################################################
 # Since the expected nominal stress across the cross section of the
@@ -261,49 +254,218 @@
 print("Stress Concentration: %.2f" % stress_con)
 
 
-# ###############################################################################
-# # Batch Analysis
-# # ~~~~~~~~~~~~~~
-# # The above script can be placed within a function to compute the
-# # stress concentration for a variety of hole diameters.  For each
-# # batch, MAPDL is reset and the geometry is generated from scratch.
-
-# ###############################################################################
-# # Run the batch and record the stress concentration
-# k_t_exp = []
-# ratios = np.linspace(0.001, 0.5, 20)
-# print('    Ratio  : Stress Concentration (K_t)')
-# for ratio in ratios:
-#     stress_con = compute_stress_con(ratio)
-#     print('%10.4f : %10.4f' % (ratio, stress_con))
-#     k_t_exp.append(stress_con)
-
-
-# ###############################################################################
-# # Analytical Comparison
-# # ~~~~~~~~~~~~~~~~~~~~~
-# # Stress concentrations are often obtained by referencing tablular
-# # results or polynominal fits for a variety of geometries.  According
-# # to Peterson's Stress Concentration Factors (ISBN 0470048247), the analytical
-# # equation for a hole in a thin plate in uniaxial tension:
-# #
-# # :math:`k_t = 3 - 3.14\frac{d}{h} + 3.667\left(\frac{d}{h}\right)^2 - 1.527\left(\frac{d}{h}\right)^3`
-# #
-# # Where:
-# #
-# # - :math:`k_t` is the stress concentration
-# # - :math:`d` is the diameter of the circle
-# # - :math:`h` is the height of the plate
-# #
-# # As shown in the following plot, ANSYS matches the known tabular
-# # result for this geometry remarkably well using PLANE183 elements.
-# # The fit to the results may vary depending on the ratio between the
-# # height and width of the plate.
-
-# # where ratio is (d/h)
-# k_t_anl = 3 - 3.14*ratios + 3.667*ratios**2 - 1.527*ratios**3
-
-# plt.plot(ratios, k_t_anl, label=r'$K_t$ Analytical')
-# plt.plot(ratios, k_t_exp, label=r'$K_t$ ANSYS')
-# plt.legend()
-# plt.show()
+###############################################################################
+# Batch Analysis
+# ~~~~~~~~~~~~~~
+# The above script can be placed within a function to compute the
+# stress concentration for a variety of hole diameters.  For each
+# batch, MAPDL is reset and the geometry is generated from scratch.
+#
+# .. note::
+#    This section has been disabled to reduce the execution time of
+#    this example. Enable it by setting ``RUN_BATCH = TRUE``
+
+RUN_BATCH = False
+
+# The function to compute the batch analysis is the following:
+
+def compute_stress_con(ratio):
+
+    notch_depth = ratio*width/2
+
+    mapdl.clear()
+    mapdl.prep7()
+
+    # Notch circle.
+    circ0_kp = mapdl.k(x=length / 2, y=width + notch_radius)
+    circ_line_num = mapdl.circle(circ0_kp, notch_radius)
+    circ_line_num = circ_line_num[2:]  # only concerned with the bottom arcs
+
+    circ0_kp = mapdl.k(x=0, y=0)
+    k1 = mapdl.k(x=0, y=-notch_depth)
+    l0 = mapdl.l(circ0_kp, k1)
+    mapdl.adrag(*circ_line_num, nlp1=l0)
+
+    circ1_kp = mapdl.k(x=length / 2, y=-notch_radius)
+    circ_line_num = mapdl.circle(circ1_kp, notch_radius)
+    circ_line_num = circ_line_num[:2]  # only concerned with the top arcs
+
+    k0 = mapdl.k(x=0, y=0)
+    k1 = mapdl.k(x=0, y=notch_depth)
+    l0 = mapdl.l(k0, k1)
+    mapdl.adrag(*circ_line_num, nlp1=l0)
+
+    rect_anum = mapdl.blc4(width=length, height=width)
+    cut_area = mapdl.asba(rect_anum, "ALL")  # cut all areas except the plate
+
+    mapdl.allsel()
+    mapdl.vext(cut_area, dz=thickness)
+
+    notch_esize = np.pi * notch_radius * 2 / 50
+    plate_esize = 0.01
+
+    mapdl.asel("S", "AREA", vmin=1, vmax=1)
+
+    mapdl.lsel("NONE")
+    for line in [7, 8, 20, 21]:
+        mapdl.lsel("A", "LINE", vmin=line, vmax=line)
+
+    mapdl.ksel('NONE')
+    mapdl.ksel('S', 'LOC', 'X', length/2 - notch_radius * 1.1,
+               length/2 + notch_radius*1.1)
+    mapdl.lslk('S', 1)
+    mapdl.lesize("ALL", notch_esize, kforc=1)
+    mapdl.lsel("ALL")
+
+    mapdl.mopt("EXPND", 0.7)  # default 1
+
+    esize = notch_esize * 5
+    if esize > thickness / 2:
+        esize = thickness / 2  # minimum of two elements through
+
+    mapdl.esize()  # this is tough to automate
+    mapdl.et(1, "SOLID186")
+    mapdl.vsweep("all")
+
+    mapdl.allsel()
+
+    # Uncomment if you want to print geometry and mesh plots.
+    # mapdl.vplot(savefig=f'vplot-{ratio}.png', off_screen=True)
+    # mapdl.eplot(savefig=f'eplot-{ratio}.png', off_screen=True)
+
+    mapdl.units("SI")  # SI - International system (m, kg, s, K).
+
+    mapdl.mp("EX", 1, 210e9)  # Elastic moduli in Pa (kg/(m*s**2))
+    mapdl.mp("DENS", 1, 7800)  # Density in kg/m3
+    mapdl.mp("NUXY", 1, 0.3)  # Poisson's Ratio
+
+    mapdl.nsel("S", "LOC", "X", 0)
+    mapdl.d("ALL", "UX")
+
+    mapdl.nsel("R", "LOC", "Y", width / 2)
+    mapdl.d("ALL", "UY")
+    mapdl.d("ALL", "UZ")
+
+    mapdl.nsel("S", "LOC", "X", length)
+    mapdl.cp(5, "UX", "ALL")
+
+    mapdl.nsel("R", "LOC", "Y", width / 2)  # selects more than one
+    single_node = mapdl.mesh.nnum[0]
+    mapdl.nsel("S", "NODE", vmin=single_node, vmax=single_node)
+    mapdl.f("ALL", "FX", 1000)
+
+    mapdl.allsel(mute=True)
+
+    mapdl.run("/SOLU")
+    mapdl.antype("STATIC")
+    mapdl.solve()
+    mapdl.finish()
+
+    result = mapdl.result
+    _, stress = result.principal_nodal_stress(0)
+    von_mises = stress[:, -1]  # von-Mises stress is the right most column
+    max_stress = np.nanmax(von_mises)
+
+    mask = result.mesh.nodes[:, 0] == length
+    far_field_stress = np.nanmean(von_mises[mask])
+
+    adj = width / (width - notch_depth * 2)
+    stress_adj = far_field_stress * adj
+
+    return max_stress / stress_adj
+
+
+###############################################################################
+# Run the batch and record the stress concentration
+
+if RUN_BATCH:
+    k_t_exp = []
+    ratios = np.linspace(0.05, 0.75, 9)
+    print('    Ratio  : Stress Concentration (K_t)')
+    for ratio in ratios:
+        stress_con = compute_stress_con(ratio)
+        print('%10.4f : %10.4f' % (ratio, stress_con))
+        k_t_exp.append(stress_con)
+
+
+###############################################################################
+# Analytical Solution
+# ~~~~~~~~~~~~~~~~~~~
+# Stress concentrations are often obtained by referencing tabular
+# results or polynominal fits for a variety of geometries.  According
+# to *Roark’s Formulas for Stress and Strain* (Warren C. Young and
+# Richard G. Budynas, Seventh Edition, McGraw-Hill) the analytical
+# equation for two U notches in a thin plate in uniaxial tension:
+#
+# .. math::
+#
+#    K_t = C_1 + C_2 \left(\dfrac{2h}{D}\right) + C_3 \left(\dfrac{2h}{D}\right)^2 + C_4 \left(\dfrac{2h}{D}\right)^3
+#
+# where:
+#
+# .. math::
+#    \begin{array}{c|c|c}
+#        & 0.1 \leq h/r \leq 2.0                 & 2.0 \leq h/r \leq 50.0 \\ \hline
+#    C_1 & 0.85 + 2.628 \sqrt{h/r} - 0.413 h/r   & 0.833 + 2.069 \sqrt{h/r} - 0.009 h/r \\
+#    C_2 & -1.119 - 4.826 \sqrt{h/r} + 2.575 h/r & 2.732 - 4.157   \sqrt{h/r} + 0.176 h/r \\
+#    C_3 & 3.563 - 0.514 \sqrt{h/r} - 2.402 h/r  & -8.859 + 5.327 \sqrt{h/r} - 0.32 h/r \\
+#    C_4 & -2.294 + 2.713 \sqrt{h/r} + 0.240 h/r & 6.294 - 3.239 \sqrt{h/r} + 0.154 h/r
+#    \end{array}
+#
+# Where:
+#
+# - :math:`K_t` is the stress concentration
+# - :math:`r` is the radius of the notch
+# - :math:`h` is the notch depth
+# - :math:`D` is the width of the plate
+#
+# In this example the ratio is given as :math:`2h/D`.
+#
+# These formulas are converted in the following function:
+
+def calc_teor_notch(ratio):
+    notch_depth = ratio*width/2
+    h = notch_depth
+    r = notch_radius
+    D = width
+
+    if 0.1 <= h/r <= 2.0:
+        c1 = 0.85 + 2.628*(h/r)**0.5 - 0.413*h/r
+        c2 = -1.119 - 4.826*(h/r)**0.5 + 2.575*h/r
+        c3 = 3.563 - 0.514*(h/r)**0.5 - 2.402*h/r
+        c4 = -2.294 + 2.713*(h/r)**0.5 + 0.240*h/r
+    elif 2.0 <= h/r <= 50.0:
+        c1 = 0.833 + 2.069*(h/r)**0.5 - 0.009*h/r
+        c2 = 2.732 - 4.157 * (h/r)**0.5 + 0.176*h/r
+        c3 = -8.859 + 5.327*(h/r)**0.5 - 0.32*h/r
+        c4 = 6.294 - 3.239*(h/r)**0.5 + 0.154*h/r
+
+    return c1 + c2*(2*h/D) + c3*(2*h/D)**2 + c4*(2*h/D)**3
+
+
+###############################################################################
+# which is used later to calculate the concentration factor for the given ratios:
+
+if RUN_BATCH:
+    print('    Ratio  : Stress Concentration (K_t)')
+    k_t_anl = []
+    for ratio in ratios:
+        stress_con = calc_teor_notch(ratio)
+        print('%10.4f : %10.4f' % (ratio, stress_con))
+        k_t_anl.append(stress_con)
+
+
+###############################################################################
+# Analytical Comparison
+# ~~~~~~~~~~~~~~~~~~~~~
+#
+# As shown in the following plot, MAPDL matches the known tabular
+# result for this geometry remarkably well using PLANE183 elements.
+# The fit to the results may vary depending on the ratio between the
+# height and width of the plate.
+
+if RUN_BATCH:
+    plt.plot(ratios, k_t_anl, label=r'$K_t$ Analytical')
+    plt.plot(ratios, k_t_exp, label=r'$K_t$ ANSYS')
+    plt.legend()
+    plt.show()