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

Options
Vishnu
Vishnu Member, Employee Posts: 214
Name Dropper First Anniversary Solution Developer Community of Practice Member First Comment
edited February 7 in Structures

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

Tagged:

Comments

  • Vishnu
    Vishnu Member, Employee Posts: 214
    Name Dropper First Anniversary Solution Developer Community of Practice Member First Comment
    edited February 12
    Options

    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 Member, Employee Posts: 214
    Name Dropper First Anniversary Solution Developer Community of Practice Member First Comment
    edited February 12
    Options


  • Vishnu
    Vishnu Member, Employee Posts: 214
    Name Dropper First Anniversary Solution Developer Community of Practice Member First Comment
    edited February 12
    Options
  • deepak
    deepak Member Posts: 2
    First Comment First Anniversary
    Options

    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
    Vishnu Member, Employee Posts: 214
    Name Dropper First Anniversary Solution Developer Community of Practice Member First Comment
    edited May 16
    Options

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

    Please check these