Abaqus - Automatic free surface generation

I am trying to find a way (Python code) such that the exterior element surfaces can automatically be detected. However, the following states that this is not possible in Abaqus/CAE (https://classes.engineering.wustl.edu/2009/spring/mase5513/abaqus/docs/v6.6/books/usb/default.htm?startat=pt01ch02s03aus16.html). Is there another way in which I can approach this problem?

My goal is to automatically put (pressure) forces on only the outside element surfaces.


For now this is my code:

from part import *

from material import *

from section import *

from assembly import *

from step import *

from interaction import *

from load import *

from mesh import *

from optimization import *

from job import *

from sketch import *

from visualization import *

from connectorBehavior import *


import subprocess #Python script in abaqus


ss = 0.5 #Seed size


session.journalOptions.setValues(replayGeometry=COORDINATE, recoverGeometry=COORDINATE)


Mdb()

#: A new model database has been created.

#: The model "Model-1" has been created.


argument = 'Circle.cae'


openMdb(

    pathName='X:/temp/30_May/'+argument)

p = mdb.models['Model-1'].parts['Part-1']

session.viewports['Viewport: 1'].setValues(displayedObject=p)

p.seedPart(size=ss, deviationFactor=0.1, minSizeFactor=0.1)

p.generateMesh()


#Create Sin Amplitude

mdb.models['Model-1'].PeriodicAmplitude(name='Sin', timeSpan=STEP,

    frequency=6.28319, start=0.0, a_0=0.0, data=((0.0, 1.0), ))


#Material

mdb.models['Model-1'].Material(name='AluminumAlloy')

mdb.models['Model-1'].materials['AluminumAlloy'].Density(table=((2770.0, ), ))

mdb.models['Model-1'].materials['AluminumAlloy'].Elastic(table=((71000000000.0,

    0.33), ))

mdb.models['Model-1'].materials['AluminumAlloy'].Damping(alpha=0.01)


#Section

c = p.cells

cells = c.findAt(((0.0, 9.854317, 0.0), ))

region = p.Set(cells=cells, name='Set-1')


mdb.models['Model-1'].HomogeneousSolidSection(name='Section-1',

    material='AluminumAlloy', thickness=None)

region = p.sets['Set-1']

p.SectionAssignment(region=region, sectionName='Section-1', offset=0.0,

    offsetType=MIDDLE_SURFACE, offsetField='',

    thicknessAssignment=FROM_SECTION)


mdb.save()

#: The model database has been saved to "C:\Users\tshar\OneDrive\Documenten\Toloue\Research Project\Abaqus Files\26_May\Circle.cae".


#Assembly-Instance is created

a = mdb.models['Model-1'].rootAssembly

a.DatumCsysByDefault(CARTESIAN)

a.Instance(name='Part-1-1', part=p, dependent=ON)


#Step is created

mdb.models['Model-1'].ImplicitDynamicsStep(name='LoadingStep',

    previous='Initial', description='SineLoadingStep', timePeriod=50.0,

    maxNumInc=10000, initialInc=0.1, minInc=0.0005, maxInc=0.1)


import random

maxelements = len(mdb.models['Model-1'].parts['Part-1'].elements)

x = random.randint(0,20)


for i in range(1,15,1):

   q = random.randint(0,maxelements)


   #Surface is created

   f = p.elements

   side1Elements1 = f[q-1:q]

   p.Surface(side1Elements=side1Elements1, name='Surf-'+str(i))

   

   #Load is created

   region = a.instances['Part-1-1'].surfaces['Surf-'+str(i)]

   mdb.models['Model-1'].Pressure(name='Load-'+str(i), createStepName='LoadingStep',

    region=region, distributionType=UNIFORM, field='', magnitude=1.0,

    amplitude='Sin')