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

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

    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: 198
    Name Dropper First Anniversary Solution Developer Community of Practice Member First Comment
    edited February 12


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