Optimization example using PyMAPDL and SciPy Optimize library

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

Below is a PyMAPDL example using FLUID218 elements (3D Hydrodynamic Bearing Element) on how you could find an optimized Inputs (Shaft center X and Y location) for a preferred FLUID forces in X and Y directions.

Please refer FLUID218 Documentation below
https://ansyshelp.ansys.com/account/secured?returnurl=/Views/Secured/corp/v241/en/ans_elem/Hlp_E_FLUID218.html?q=FLUID218

Tagged:

Comments

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

    The first step in optimization is formulating an objective function, in this case the function named extract_fluid_forces is evaluated by using PyMAPDL

    scipy.optimize ,minimize is used for optimization .please refer the documentation for it below:

    https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

    Example code below

    # Import necessary libraries
    import matplotlib.pyplot as plt
    import numpy as np
    from ansys.dpf import core as dpf
    from scipy.optimize import minimize
    
    from ansys.mapdl.core import launch_mapdl
    
    # Define the path to the result file
    rstpath = r"D:\apdl\test\file.rst"
    
    # Launch MAPDL with specified run location and number of processors
    mapdl = launch_mapdl(run_location=r"D:\apdl\test",nproc=1,override=True)
    
    mapdl.clear("All")
    mapdl.prep7()
    mapdl.et('', 'FLUID218')
    mapdl.mp(lab='VISC', mat='1', c0='0.135')
    mapdl.mp(lab='DENS', mat='1', c0='910')
    mapdl.r(nset='1', r1='0', r2='10')
    mapdl.cyl4(xcenter='', ycenter='', rad1='10', theta1='0', rad2='', theta2='90', depth='50')
    mapdl.cyl4(xcenter='', ycenter='', rad1='10', theta1='90', rad2='', theta2='180', depth='50')
    mapdl.cyl4(xcenter='', ycenter='', rad1='10', theta1='180', rad2='', theta2='270', depth='50')
    mapdl.cyl4(xcenter='', ycenter='', rad1='10', theta1='270', rad2='', theta2='360', depth='50')
    mapdl.vdele(nv1='ALL', nv2='', ninc='', kswp='0')
    mapdl.run('CSYS,1')
    mapdl.asel('U', 'LOC', 'X', 9, 10)
    mapdl.adele('ALL')
    mapdl.allsel()
    mapdl.mshape(0, '2D')
    mapdl.mshkey(1)
    mapdl.lesize('ALL', size='', angsiz='',  ndiv=20)
    mapdl.amesh('ALL')
    mapdl.eplot()
    mapdl.nummrg('ALL')
    clearance = 0.1
    mapdl.rmodif(1, 8, clearance)
    mapdl.rmodif(1, 9, clearance)
    mapdl.rmodif(1, 10, clearance)
    mapdl.rmodif(1, 11, clearance)
    mapdl.nsel('S', 'LOC', 'Z', '', 0)
    mapdl.d('ALL', 'PRES', 0)
    mapdl.allsel()
    mapdl.nsel('S', 'LOC', 'Z', 50, 50)
    mapdl.d('ALL', 'PRES', 100)
    mapdl.nsel('S', 'LOC', 'Y', -20, 20)
    mapdl.d('ALL', 'PRES', 100)
    mapdl.nsel('S', 'LOC', 'Y', 160, 200)
    mapdl.d('ALL', 'PRES', 100)
    mapdl.allsel()
    mapdl.omega('', '', 10)
    mapdl.finish()
    
    # Define a function to extract fluid forces
    def extract_fluid_forces(xoffset,yoffset):
            mapdl.prep7()
            mapdl.rmodif(1, 3,xoffset)
            mapdl.rmodif(1, 4,yoffset)
            mapdl.solution()
            mapdl.solve()       
            # data source
            my_data_sources = dpf.DataSources(rstpath)
    
            # model
            my_model = dpf.Model(rstpath)
    
            # mesh
            my_mesh = my_model.metadata.meshed_region
    
            # Field Containers at all mesh
            # Get result operator
            result_operator = dpf.operators.result.__getattribute__("nmisc")
    
            result_field_container = result_operator(
                time_scoping=[1], data_sources=my_data_sources, mesh=my_mesh,item_index=11
            ).outputs.fields_container.get_data()
    
            # First Field usually field corresponding to first time scoping
            resultfield = result_field_container[0]
    
            FXTOTAL=sum(resultfield.data_as_list)
    
            result_field_container = result_operator(
                time_scoping=[1], data_sources=my_data_sources, mesh=my_mesh,item_index=12
            ).outputs.fields_container.get_data()
    
            # First Field usually field corresponding to first time scoping
            resultfield = result_field_container[0]
    
            FYTOTAL=sum(resultfield.data_as_list)
    
            fluid_forces = np.array([FXTOTAL,FYTOTAL])
            #plt.bar(range(len(fluid_forces)), fluid_forces)
    
            return fluid_forces
    
    
    # preferred forces
    preferred_fx = 10000  # replace with your preferred FXtotal
    preferred_fy = 90000  # replace with your preferred FYtotal
    
    # Define objective function
    def objective(offsets):
        try:
            [fx, fy] = extract_fluid_forces(offsets[0],offsets[1])
            return (((fx - preferred_fx)**2 + (fy - preferred_fy)**2))
        except:  # catch all exceptions
            return 9E50  # return a large value
    
    
    # initial guess for x and y offsets
    x0 = [0.0,0.0]
    
    # Minimize the objective function
    result = minimize(objective, x0,method='nelder-mead',options={'maxfev':1000,'maxiter':1000,'disp': True,'fatol' : 1e-10, 'xatol' : 1e-4})
    
    # Get the optimized x and y offsets
    x, y = result.x
    fx,fy = extract_fluid_forces(x, y)
    
    # Exit MAPDL
    mapdl.exit()
    
    # Print the optimized values
    print("Optimized X and Y Offsets are %f and %f"%(x,y))
    print("The Fx,Fy forces at these X and Y offsets are %f and %f "%(fx,fy))
    
    
    # Define the labels for the bars
    labels = ['FX', 'FY']
    
    # Create a figure and a set of subplots
    fig, ax = plt.subplots()
    
    # Define the width of the bars
    width = 0.35
    
    # Define the positions for the groups of bars
    x = np.arange(len(labels))
    
    # Plot the preferred and obtained values as bars next to each other
    ax.bar(x - width/2, [preferred_fx,preferred_fy], width, label='Preferred')
    ax.bar(x + width/2, [fx,fy], width, label='optimized')
    
    # Add some text for labels, title and custom x-axis tick labels, etc.
    ax.set_ylabel('Values')
    ax.set_title('Preferred and Obtained Values of FX and FY')
    ax.set_xticks(x)
    ax.set_xticklabels(labels)
    ax.legend()
    
    
    # Display the plot
    plt.show()
    
    
  • Vishnu
    Vishnu Member, Employee Posts: 198
    Name Dropper First Anniversary Solution Developer Community of Practice Member First Comment
    edited February 6

    Output --->>>