How to create beam bending moment (BMD) and shear force diagram (SFD) using PyMAPDL and PyDPF

Vishnu

shear force and bending moment diagram example below uses BEAM188 elements .

example uses pyMAPDL to run the simulation
Pydpf to extract result
matplotlib to plot the SF and BM diagrams

You can alter the beam SECTYPE,SECDATA, fix_location and force location as input arguments of the function

Vishnu

Inputs and import statements

#Import statements
from ansys.mapdl.core import launch_mapdl
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable
from matplotlib.colors import Normalize
from matplotlib import ticker
from ansys.dpf import core as dpf
import os
run_location=r"D:\apdl\test\4"
nproc=4
override=True
mapdl = launch_mapdl(run_location=run_location, nproc=nproc, override=override)

Creation of model using pyMAPDL

def create_model_setup(SECDATA, 
                       BEAM_LENGTH, 
                       LINE_DIV,
                       LENGTH_DIV, 
                       YM, 
                       MU, 
                       SECTYPE, 
                       SECSUBTYPE):
    mapdl.clear("All")
    mapdl.prep7()
    mapdl.et('', 'BEAM188')
    mapdl.mp(lab='EX', mat='1', c0='%s'%(YM))
    mapdl.mp(lab='PRXY', mat='1', c0='%s'%(MU))
    mapdl.sectype(secid='1', type_='BEAM', subtype='%s'%(SECTYPE))
    mapdl.secdata(tuple(SECDATA))
    mapdl.seccontrol(1E15,1E15)
    mapdl.k(npt='1', x='0', y='0', z='0')
    mapdl.k(npt='2', x='%s'%(BEAM_LENGTH), y='0', z='0')
    mapdl.l(p1='1', p2='2', ndiv='%s'%(LINE_DIV))
    mapdl.lesize(nl1='ALL', size='', angsiz='', ndiv='%s'%(LENGTH_DIV))
    mapdl.lmesh(nl1='ALL')

Adding Loads and Bc and Solving them

def run_simulation(fix_loc,
                   pin_loc, 
                   concentrated_loc, 
                   concentrated_force,
                   distributed_loc,
                   distributed_force,
                   elem_plot_flag):
    mapdl.solution()
    mapdl.ddele(node='ALL', lab='ALL')
    mapdl.fdele(node='ALL', lab='ALL')
    mapdl.sfedele(elem='ALL', lkey='ALL', lab='ALL')
    for fix_lo in fix_loc:
        fix_node = mapdl.queries.node(fix_lo[0], fix_lo[1], fix_lo[2])
        mapdl.d(node='%s'%(fix_node), lab='ALL', value='0')
    for pin_lo in pin_loc:
        pin_node = mapdl.queries.node(pin_lo[0], pin_lo[1], pin_lo[2])
        mapdl.d(node='%s'%(pin_node), lab='UX', value='0')
        mapdl.d(node='%s'%(pin_node), lab='UY', value='0')
        mapdl.d(node='%s'%(pin_node), lab='UZ', value='0')
        mapdl.d(node='%s'%(pin_node), lab='ROTY', value='0')
        mapdl.d(node='%s'%(pin_node), lab='ROTZ', value='0')
    for count, force_lo in enumerate(concentrated_loc):
        force_node = mapdl.queries.node(force_lo[0], force_lo[1], force_lo[2])
        mapdl.f(node='%s'%(force_node), lab='FY', value=concentrated_force[count])
    for count, force_lo in enumerate(distributed_loc):
        mapdl.esel(type_='S', item='CENT', comp='X', vmin=force_lo[0], vmax=force_lo[1])
        mapdl.sfbeam(elem='ALL', lkey='2', lab='PRES', vali=distributed_force[count])   
        mapdl.allsel(labt='ALL', entity='ALL')
    mapdl.allsel(labt='ALL', entity='ALL')
    mapdl.outres(item='ALL', freq='ALL')
    if elem_plot_flag == True:
        mapdl.eplot(plot_bc=True, bc_labels=["UX", "UY", "UZ", "FX", "FY", "FZ"], plot_bc_labels=True)
    mapdl.solve()
    mapdl.save(fname='test', ext='db')

Extracting shear force and bending moments using DPF and plot using matplotlib

def plot_bending_moment_shear_force(rstpath, 
                                    smiscI, 
                                    smiscJ, 
                                    plotname):
    # data source
    my_data_sources = dpf.DataSources(rstpath)
    # model
    my_model = dpf.Model(rstpath)
    # mesh
    my_mesh = my_model.metadata.meshed_region
    # Get result operator
    result_operator = dpf.operators.result.__getattribute__("smisc")
    result_field_container = result_operator(
        time_scoping=[1], data_sources=my_data_sources, mesh=my_mesh, item_index=smiscI
    ).outputs.fields_container.get_data()
    # First Field usually field corresponding to first time scoping
    resultfield = result_field_container[0]
    element_ids = resultfield.scoping.ids
    i_nodes_coordinates = [my_mesh.elements.element_by_id(
        x).nodes[1].coordinates[0] for x in element_ids]
    BMY_T_I = resultfield.data_as_list
    result_field_container = result_operator(
        time_scoping=[1], data_sources=my_data_sources, mesh=my_mesh, item_index=smiscJ
    ).outputs.fields_container.get_data()
    # First Field usually field corresponding to first time scoping
    resultfield = result_field_container[0]
    element_ids = resultfield.scoping.ids
    j_nodes_coordinates = [my_mesh.elements.element_by_id(
        x).nodes[1].coordinates[0] for x in element_ids]
    BMY_T_J = resultfield.data_as_list
    # Assuming i_nodes_coordinates, j_nodes_coordinates, BMY_T_I, BMY_T_J are defined
    x = np.concatenate((i_nodes_coordinates, j_nodes_coordinates))
    y = np.concatenate((BMY_T_I, BMY_T_J))
    # Create a colormap
    cmap = plt.get_cmap('viridis')
    norm = Normalize(vmin=np.min(y), vmax=np.max(y))
    # Create a ScalarMappable
    sm = ScalarMappable(cmap=cmap, norm=norm)
    sm.set_array([])
    # Create the bar chart
    fig, ax = plt.subplots(figsize=(15, 5))  # Adjust the figure size
    ax.bar(x, y, color=cmap(norm(y)))
    ax.set_xlabel('X_distance')
    ax.set_title(plotname+' Diagram')
    # Set y-axis limits
    ax.set_ylim([np.min(y), np.max(y)])
    # Create an array of 10 evenly spaced values between min(y) and max(y)
    ticks = np.linspace(np.min(y), np.max(y), 10)
    # Add a colorbar with specific ticks
    cbar = fig.colorbar(sm, ax=ax, ticks=ticks)
    cbar.set_label(plotname)
    # Format the colorbar tick labels to have two decimal places
    cbar.formatter = ticker.FuncFormatter(lambda x, pos: '{:.4f}'.format(x))
    cbar.update_ticks()
    plt.show()

Function calls Below

### Call the defined functions####
# create model
create_model_setup(
    SECDATA=[5, 8, 2],  # SECDATA = [BEAM_RADIUS,CIRCUM_DIV,RADIAL_DIV]
    BEAM_LENGTH=100,
    LINE_DIV=50,
    LENGTH_DIV=1,
    YM=2E8,
    MU=0.3,
    SECTYPE="CSOLID",
    SECSUBTYPE="CIRBEAM",
)
#Run simulation
run_simulation(
    fix_loc=[[10, 0, 0], [90, 0, 0]],
    pin_loc=[[50, 0, 0]],
    concentrated_loc = [[30, 0, 0], [80, 0, 0],[100, 0, 0]],
    concentrated_force = [1000, 500,500],
    distributed_loc=[[0,100]],
    distributed_force = [1000],
    elem_plot_flag = True
)
# Bending Moment diagram
plot_bending_moment_shear_force(
    rstpath=os.path.join(run_location, "file.rst"),
    smiscI=32,
    smiscJ=37,
    plotname="Bending Moment"
)
# Shear Force diagram
plot_bending_moment_shear_force(
    rstpath=r"D:\apdl\test\4\file.rst",
    smiscI=6,
    smiscJ=19,
    plotname="Shear Force"
)
mapdl.exit()

Vishnu

Vishnu

Also for further reference please refer BEAM188 documentation below,
https://ansyshelp.ansys.com/account/secured?returnurl=/Views/Secured/corp/v231/en/ans_elem/Hlp_E_BEAM188.html

deepak

Hi,
i am getting error while running above code.

The operator SMISC doesn't exist in the registry. Check its spelling in the documentation or verify its availability in your loaded plugins. The current available operator names can be accessed using 'available_operator_names' method."
how to rectify it?
thanks

Vishnu

@deepak
In recent releases its
dpf.operators.result.smisc
in earlier releases its
dpf.operators.result.mapdl.smisc

Please check these