Problem in boundary condition for Elasticity Equation (2D IVP) in a disk

[Arnab Roy]

Hi all,

I am new to Dedalus and trying to solve a linear elasticity problem in a disk (2D IVP). I am trying to give stress free boundary condition in terms of both "sigma(r=1)=h" (both sigma and h are tensors) or in terms of traction vector  "radial(sigma(r=1))=m" (m is vector) or separately I can do 'traction_vector(r=1)=m'. But I am always getting the error <dedalus.core.coords.AzimuthalCoordinate object at 0x13e1049d0>. Can anyone suggest or help me to figure out what is the problem here? I searched for this error and it feels like a problem coming from NCC along 'phi'(which is the Fourier domain here)..Is there a way to solve the problem? It will be very helpful to get any little help.

Another query (to understand things better) - 1. What are edge basis (edge = disk.edge) or S1_basis or boundary basis actually mean? 2. How to define tensors component wise like d3.radial(d3.radial(sigma)) is this one correct or d3.radial(d3.radial(sigma),index=1) is correct for representing sigma_rr ? And similarly what will be correct for sigma_rphi? Any response will be highly appreciated.

Thank you.

My code is attached here.

# Elasticity Vector eqn  (2D IVP)
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
import dedalus.public as d3
import logging
logger = logging.getLogger(__name__)
# Parameters
Kb = 10    #Bulk Modulus
mu = 1     #Shear modulus
Nphi, Nr = 128,128
dealias = 1
stop_sim_time = 0.01
timestepper = d3.SBDF2
timestep = 1e-6
dtype = np.float64

# Bases
coords = d3.PolarCoordinates('phi', 'r')
dist = d3.Distributor(coords, dtype=dtype)
disk = d3.DiskBasis(coords, shape=(Nphi, Nr), radius=1, dealias=dealias, dtype=dtype)
S1_basis = basis.S1_basis(radius=1)
edge = disk.edge

# Substitutions
phi, r = dist.local_grids(disk,scales=disk.dealias)
lift = lambda A: d3.Lift(A, disk, -1)
t = dist.Field()
der = lambda f: d3.grad(f)
#graddiv = lambda g: d3.grad(d3.div(g))

# Fields
u = dist.VectorField(coords,name='u',bases=disk)
tau_u = dist.VectorField(coords, name='tau_u', bases=edge)
strain = dist.TensorField((coords, coords), name= 'strain', bases=disk)
#Deviatoric strain tensor
dev_strain = dist.TensorField((coords, coords), name= 'dev_strain', bases=disk)
##Stress tensor
sigma = dist.TensorField((coords, coords), name= 'sigma', bases=disk)
#traction 
traction = dist.VectorField(coords, name='traction', bases= disk)

#Identity tensor in polar coordinates
I = dist.TensorField((coords, coords), name= 'I', bases=disk)
I['g'][0,0] = 1
I['g'][1,1] = r*r
I['g'][0,1] = 0
I['g'][1,0] = 0

# Radial unit vecor
er = dist.VectorField(coords, bases=disk.radial_basis)
er['g'][0] = 0
er['g'][1] = 1
#Tangent unit vector
ephi = dist.VectorField(coords, bases=disk.radial_basis)
ephi['g'][0] = 1
ephi['g'][1] = 0

#equation
strain = (d3.grad(u) + d3.trans(d3.grad(u)))/2
dev_strain  = strain - (d3.div(u)*I)/2
sigma = Kb*((d3.div(u))*I) + 2*mu*(dev_strain)
#shear_stress = d3.angular(d3.radial(strain_rate(r=1), index=1))
sigma_rr = d3.radial(d3.radial(strain(r=1),index=1))
sigma_rt = d3.azimuthal(d3.radial(strain(r=1),index=1))
strain_rr = d3.radial(d3.radial(strain(r=1)))
print(sigma_rr)
print(sigma_rt)

# h field
h = dist.TensorField((coords, coords), name= 'h', bases=disk)
h['g'][1,1] = 0
h['g'][0,0] = 0
h['g'][0,1] = 0
h['g'][1,0] = 0

traction = er@sigma

m['g'][1]=-1
m['g'][0]=0
# Problem, two variable
problem = d3.IVP([u,tau_u],
                  time=t, namespace=locals())
problem.add_equation("dt(u) - Kb*grad(div(u))- mu*lap(u)+ lift(tau_u) = 0")

# Add the B.C

problem.add_equation("sigma(r=1) = h")
#problem.add_equation("radial(sigma(r=1))=m")
#problem.add_equation("sigma_rr = 0")
#problem.add_equation("sigma_rt = 0")
#problem.add_equation("traction = dot(sigma,er)")
#problem.add_equation("traction(r=1)=m")

# Initial conditions
u.fill_random('g', seed=42, distribution='standard_normal') # Random noise
u.low_pass_filter(scales=0.25) # Keep only lower fourth of the modes

# Solver
solver = problem.build_solver(timestepper)   # Build Solver
solver.stop_sim_time = stop_sim_time         # Set integration limits
#print(grad(dot(u,er))(r=1))
# Analysis
snapshots = solver.evaluator.add_file_handler('snaps', sim_dt=0.001, max_writes=10)
snapshots.add_task(u, name='u')

# Main loop
try:
    logger.info('Starting main loop')
    while solver.proceed:
        # step forward
        solver.step(timestep)
        #perform some analysis
        #print(np.mean(u['g']),np.std(u['g']))
        if (solver.iteration-1) % 10 == 0:
            logger.info("Iteration=%i, Time=%e, dt=%e" %(solver.iteration, solver.sim_time, timestep))
except:
    logger.error('Exception raised, triggering end of main loop.')
    raise
finally:
    solver.log_stats()