how can I do the MD of twisting balloon

Hi, I am Yu-Chuan Cheng from NTHU,
I have a problem about twisting balloon.
I have a long cylinder balloon membrane, and I want to twist it.

How can I incorporate air conservation in simulation ?

the following link is the video of twisting balloon,

Best regards

the lammps code is following:

clear

Initialization ---------------------------------------------

units lj
boundary p p p # p = periodic
atom_style angle # only bonds and angle
neighbor 0.3 bin # 0.3 = extra distance beyond force cutoff (distance units)

bin = binning an operation that scales linearly with the number of atoms per processor

neigh_modify delay 5

Parameters ---------------------------------------------------

variable step equal step
variable k_bond equal 2000000.0
variable k_angle equal 50000.0
variable randomseed equal 239710 # 6位數字
variable epsilon equal 1.0
variable sigma1 equal 1.0
variable sigma2 equal 0.01
variable sigma3 equal 1.0
variable cutoff1 equal \${sigma1} # lj potential cutoff
variable lowerb equal -100.0
variable upperb equal 100.0
variable dumpstep equal 500 # per dumpstep output a dump
variable thermostep equal 500 # per thermostep output a thermo
variable totalstep1 equal 150000
variable timestep equal 0.00001
variable PEStep equal 23000

variable force equal 5.0
variable VELOCITY equal 1.0
variable R0 equal 11

Configuration ------------------------------------------------

region 1 block {upperb} EDGE EDGE EDGE EDGE EDGE region 2 block EDGE {lowerb} EDGE EDGE EDGE EDGE

region 3 cylinder x 0 0 v_R0 {lowerb}-10 {lowerb}+10 open 2
region 4 cylinder x 0 0 v_R0 {upperb}-10 {upperb}+10 open 1

group 1 region 1
group 2 region 2
group 3 id <= 16758
group 4 id > 16758

Force field ---------------------------------------------------

pair_style lj/cut {cutoff1} # lj potantial / 1.122 = result of differential pair_coeff 1 1 {epsilon} {sigma1} {cutoff1} # atom type 2 term [epsilon] [sigma] [cutoff]
pair_coeff 1 2 {epsilon} {sigma2} {cutoff1} # atom type 2 term [epsilon] [sigma] [cutoff] pair_coeff 2 2 {epsilon} {sigma3} {cutoff1} # atom type 2 term [epsilon] [sigma] [cutoff]
bond_style harmonic # potential = E = K(r-r0)^2
bond_coeff 1 {k_bond} 1.0 # [bond type] [K] [r0] #angle_style harmonic # potential = E = K(theta-theta0)^2 #angle_coeff 1 {k_angle} 180 # [angle type] [K] [theta0(degree)]

angle_style area_volume
angle_coeff 1 1 1 1 1 1 # type ka a0 kv v0 kl

Action ---------------------------------------------------------

fix a 1 nvt temp 0.0001 0.0001 1.0
fix b 2 nvt temp 100 100 1.0
#fix w all wall/region 3 lj93 {epsilon} {sigma1} {cutoff1} # [regionID] [epsilon] [sigma] [cutoff] #fix w all wall/region 4 lj93 {epsilon} {sigma1} {cutoff1} # [regionID] [epsilon] [sigma] [cutoff]

fix left 1 setforce NULL 0.0 0.0

fix right 2 setforce NULL 0.0 0.0

angular velocity type------------------------

#variable Vyu atom \${VELOCITY}*z/(y^2+z^2)^0.5

#variable Vzu atom -\${VELOCITY}*y/(y^2+z^2)^0.5

#variable Vyl atom -\${VELOCITY}*z/(y^2+z^2)^0.5

#variable Vzl atom \${VELOCITY}*y/(y^2+z^2)^0.5

#velocity 1 set NULL v_Vyu v_Vzu
#velocity 2 set NULL v_Vyl v_Vzl

torque type----------------------------------

compute a 1 property/atom y
compute b 1 property/atom z
compute c 2 property/atom y
compute d 2 property/atom z

variable ycoordu atom -{force}*c_a variable zcoordu atom {force}*c_b
variable ycoordl atom {force}*c_c variable zcoordl atom -{force}*c_d

fix aleft 1 addforce 0.0 v_zcoordu v_ycoordu
fix aright 2 addforce 0.0 v_zcoordl v_ycoordl

Compute and Dump -----------------------------------------------

compute 4 all pe/atom angle # per-atom energy is calculated by the various angle
compute 3 all pe/atom bond # per-atom energy is calculated by the various bond

#dump 1 all custom {totalstep1} pe%.R41L82_Kb{k_bond}_Ka\${k_angle}_s\${totalstep1} id c_3 c_4
dump 2 all custom \${dumpstep} atomlocal.* id x y z

Run -------------------------------------------------------------

timestep {timestep} thermo {thermostep}
thermo_style custom step temp eangle ebond

restart \${totalstep1} twistcylinder.example

run \${totalstep1}

I know people have modeled flexible membranes with LAMMPS,
e.g. for closed volumes like a biological cell. These are coarse-grained
models where the membrane is simply beads bonded together, like
a 2d hexagonal lattice of harmonic bonds. If you make the bead separation small enough then
nothing inside will escape. However you don’t need anything inside
to “inflate” the membrane. The fact that it wants to be a 2d plane
is enough to make it balloon like.

Whether you can twist it to create
balloon animal shapes in an MD simulation is a separate Q.
Just apply a spatial/time dependent fold-into-an-elephant external force field