### How to adapt undocumented tutorial dg_advection_diffusion from 2D to 3D simulation

203
views
0
5 months ago by

Tested: ubuntu 16.04, official repo, 2017.1 and 2017.2

The tutorial is based on a 2D case, it is fine to get some result on my PC.

but after adapted into 3D heat transfer with convection PDE, the result is NAN.
Actually, I also try CG version of advection tutorial using SPUG stablization, it works perfect for 2d, but failed in 3D with Nan result.

I wonder how to adapt the tutorial code into 3D application.

Thanks very much!

The minimium working code are attached
from __future__ import print_function, division
import math
import numpy as np

from dolfin import *
from mshr import *

""" pin on disc tribology heat transfer, 3D quasi-static simulation
# see paper:  Qingfeng Xia, D.G., Andrew Owen, Gervas Franceschini.
Quasi-static thermal modelling of multi-scale sliding contact for unlubricated brush seal materials.
in Proceedings of ASME Turbo Expo 2018: Power for Land, Sea and Air. 2018. Lillestrom (Oslo), Norway.
"""

using_3D = True
r_p = 1.5e-3
thickness_p = 0.002
r_d = 1.5e-2
thickness_d = 1e-3
eccentricity = 1e-2  # track position

z_original = 0
disc_zmin, disc_zmax = z_original-thickness_d, z_original
pin_zmin, pin_zmax = z_original, z_original + thickness_p
#heat_source_zmin, heat_source_zmax = z_original, z_original+heat_source_depth

pin_pos = Point(0,eccentricity,pin_zmin)  # theta = 0
pin_pos_top = Point(0, eccentricity, pin_zmax)

if using_3D:
geom_disc = Cylinder(Point(0,0, disc_zmax), Point(0,0, disc_zmin), r_d, r_d) # disc height in Z- direction, z<0
geom_pin = Cylinder(pin_pos, pin_pos_top, r_p, r_p)
bpin = AutoSubDomain(lambda x, on_boundary:  x[2]>z_original and on_boundary)
bdisc = AutoSubDomain(lambda x, on_boundary: (near(x[2], z_original) or x[2]<z_original) and on_boundary)
else:
geom_disc = Circle(Point(0,0, disc_zmax), r_d) # disc height in Z- direction, z<0
geom_pin = Circle(pin_pos, r_p)
bpin = AutoSubDomain(lambda x, on_boundary:  x[1]>eccentricity-0.5*r_p and x[1]<eccentricity+0.5*r_p and on_boundary)
#bdisc = AutoSubDomain(lambda x, on_boundary: (near(x[2], z_original) or x[2]<z_original) and on_boundary)
geom = CSGUnion(geom_disc, geom_pin)
mesh = generate_mesh(geom, 60)

parameters["ghost_mode"] = "shared_facet"

#################################################

fe_degree = 1
# Defining the function spaces
V_dg = FunctionSpace(mesh, "DG", fe_degree)
V_cg = FunctionSpace(mesh, "CG", 1)
V_u  = VectorFunctionSpace(mesh, "CG", 2)

T_hot = 360
T_cold = 300
T_ambient = 300
htc =200
omega = 0.1  # angular velocity
velocity_code = '''

class Velocity : public Expression
{
public:

// Create expression with any components
Velocity() : Expression(3) {}

// Function for evaluating expression on each cell
void eval(Array<double>& values, const Array<double>& x, const ufc::cell& cell) const
{
const double x0 = %f;
const double y0 = %f;
const double omega = %f;

double r = sqrt((x[0]-x0)*(x[0]-x0) + (x[1]-y0)*(x[1]-y0));
double v = omega * r;
double a = atan2((x[1]-y0), (x[0]-x0));
values[0] = -v * sin(a);
values[1] = v * cos(a);
values[2] = 0.0;
}
};
'''%(0, 0, omega)
vector_degree = 2
if using_3D:
v_e = Expression(cppcode=velocity_code, degree=vector_degree)
v_e = Expression(("1.0", "0", "0"), degree=vector_degree)
else:
v_e = Expression(("1.0", "0"), degree=vector_degree)
u = interpolate(v_e, V_u)

# Create velocity Function
#u = Function(V_u, "../unitsquare_64_64_velocity.xml.gz")

# Test and trial functions
v   = TestFunction(V_dg)
phi = TrialFunction(V_dg)

# Diffusivity
kappa = Constant(1e-3)

# Source term
f = Constant(0.0)

# Penalty term
alpha = Constant(5.0)

# Mesh-related functions
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2

# ( dot(v, n) + |dot(v, n)| )/2.0
un = (dot(u, n) + abs(dot(u, n)))/2.0

# Bilinear form

a_fac = kappa*(alpha/h('+'))*dot(jump(v, n), jump(phi, n))*dS \

a_vel = dot(jump(v), un('+')*phi('+') - un('-')*phi('-') )*dS  + dot(v, un*phi)*ds

a = a_int + a_fac + a_vel

# Linear form
L = v*f*dx

# Set up boundary condition (apply strong BCs)
g = Expression("350", degree=2)
bc = DirichletBC(V_dg, g, bpin, "geometric")

# Solution function
phi_h = Function(V_dg)

# Assemble and apply boundary conditions
A = assemble(a)
b = assemble(L)
bc.apply(A, b)

# Solve system
solve(A, phi_h.vector(), b)

# Project solution to a continuous function space
up = project(phi_h, V=V_cg)

file = File("temperature.pvd")
file << up

# Plot solution
plot(up)
import matplotlib.pyplot as plt
plt.show()​
Community: FEniCS Project

0
5 months ago by
I had given other trials on this piece of code. It is found for very simple 3D case, and proper selected diffusion coefficient, and small velocity magnitude, it is possible to get non-NaN solution for 3D cases.

I wonder if there is some magic adjustment to deal with Peclet number, and how to select the mesh size and penalty number?

And is there a way to transform stiffness matrix into numpy array, to check if there is any NaN value generated during assebling?

Thanks

from __future__ import print_function, division
import math,sys
import numpy as np

from dolfin import *
from mshr import *

""" pin on disc tribology heat transfer, 3D quasi-static simulation
# see paper:  Qingfeng Xia, D.G., Andrew Owen, Gervas Franceschini.
Quasi-static thermal modelling of multi-scale sliding contact for unlubricated brush seal materials.
in Proceedings of ASME Turbo Expo 2018: Power for Land, Sea and Air. 2018. Lillestrom (Oslo), Norway.
"""

using_3D = True
using_simple_geomtry= True and using_3D
fe_degree = 1 # 2 is possible but use more memory,
parameters["ghost_mode"] = "shared_facet"

r_p = 2.5e-2
thickness_p = 0.02
r_d = 1.5e-1
thickness_d = 1e-2
eccentricity = 1e-1  # track position

z_original = 0
disc_zmin, disc_zmax = z_original-thickness_d, z_original
pin_zmin, pin_zmax = z_original, z_original + thickness_p
#heat_source_zmin, heat_source_zmax = z_original, z_original+heat_source_depth

pin_pos = Point(0,eccentricity,pin_zmin)  # theta = 0
pin_pos_top = Point(0, eccentricity, pin_zmax)

if using_simple_geomtry:
res = int(10 / fe_degree)
mesh = UnitCubeMesh(res*2,res*2,res)
bpin = AutoSubDomain(lambda x, on_boundary:  near(x[2], 1) and x[0]>0.5 and x[1]>0.5 and on_boundary)
bdisc = AutoSubDomain(lambda x, on_boundary: near(x[2], 0) and on_boundary)
else:
if using_3D:
geom_disc = Cylinder(Point(0,0, disc_zmax), Point(0,0, disc_zmin), r_d, r_d) # disc height in Z- direction, z<0
geom_pin = Cylinder(pin_pos, pin_pos_top, r_p, r_p)
bpin = AutoSubDomain(lambda x, on_boundary:  x[2]>z_original and on_boundary)
bdisc = AutoSubDomain(lambda x, on_boundary: (near(x[2], disc_zmin)) and on_boundary)
else:
geom_disc = Circle(Point(0,0, disc_zmax), r_d) # disc height in Z- direction, z<0
geom_pin = Circle(pin_pos, r_p)
bpin = AutoSubDomain(lambda x, on_boundary:  x[1]>eccentricity-0.5*r_p and x[1]<eccentricity+0.5*r_p \
and x[0]>eccentricity-0.5*r_p and x[0]<eccentricity+0.5*r_p and on_boundary)
#bdisc = AutoSubDomain(lambda x, on_boundary: (near(x[2], z_original) or x[2]<z_original) and on_boundary)
geom = CSGUnion(geom_disc, geom_pin)
mesh = generate_mesh(geom, 60)

#################################################

# Defining the function spaces
V_dg = FunctionSpace(mesh, "DG", fe_degree)
V_cg = FunctionSpace(mesh, "CG", fe_degree)
vector_degree = fe_degree
V_u  = VectorFunctionSpace(mesh, "CG", vector_degree)
print("functionspace.dim()", V_dg.dim())
print("element().signature()", V_dg.element().signature())
print("value rank", V_dg.element().value_rank())
#V_dg.element().value_dimension()

#sys.exit(0)

T_hot = 360
T_cold = 300
T_ambient = 300
htc =200
omega = 100  # angular velocity,  stable for 2D,  only a small omege can generate some result for 3D
if using_3D:
omega = 0.1
velocity_code_3D = '''

class Velocity : public Expression
{
public:

// Create expression with any components
Velocity() : Expression(3) {}

// Function for evaluating expression on each cell
void eval(Array<double>& values, const Array<double>& x, const ufc::cell& cell) const
{
const double x0 = %f;
const double y0 = %f;
const double omega = %f;

double r = sqrt((x[0]-x0)*(x[0]-x0) + (x[1]-y0)*(x[1]-y0));
double v = omega * r;
double a = atan2((x[1]-y0), (x[0]-x0));
values[0] = -v * sin(a);
values[1] = v * cos(a);
values[2] = 0.0;
}
};
'''%(0, 0, omega)

velocity_code_2D = '''

class Velocity : public Expression
{
public:

// Create expression with any components
Velocity() : Expression(2) {}

// Function for evaluating expression on each cell
void eval(Array<double>& values, const Array<double>& x, const ufc::cell& cell) const
{
const double x0 = %f;
const double y0 = %f;
const double omega = %f;

double r = sqrt((x[0]-x0)*(x[0]-x0) + (x[1]-y0)*(x[1]-y0));
double v = omega * r;
double a = atan2((x[1]-y0), (x[0]-x0));
values[0] = -v * sin(a);
values[1] = v * cos(a);
}
};
'''%(0, 0, omega)

if using_3D:
v_e = Expression(cppcode=velocity_code_3D, degree=vector_degree)
if using_simple_geomtry:
v_e = Expression(("1.0", "-1", "-1"), degree=vector_degree)
else:
v_e = Expression(cppcode=velocity_code_2D, degree=vector_degree)
#v_e = Expression(("-1.0", "-1"), degree=vector_degree)
u = interpolate(v_e, V_u)
plot(u)

# Create velocity Function
#u = Function(V_u, "../unitsquare_64_64_velocity.xml.gz")

# Test and trial functions
v   = TestFunction(V_dg)
phi = TrialFunction(V_dg)

# Diffusivity
kappa = Constant(1e-2)

# Source term
f = Constant(1.0)

# Penalty term
alpha = Constant(5.0)

# Mesh-related functions
n = FacetNormal(mesh)
h = CellSize(mesh)  # CellDiameter() in later version?
h_avg = (h('+') + h('-'))/2
#plot(h)

# ( dot(v, n) + |dot(v, n)| )/2.0
un = (dot(u, n) + abs(dot(u, n)))/2.0
#plot(un)

# Bilinear form

a_fac = kappa*(alpha/h('+'))*dot(jump(v, n), jump(phi, n))*dS \

a_vel = dot(jump(v), un('+')*phi('+') - un('-')*phi('-') )*dS  + dot(v, un*phi)*ds
a = a_int  + a_fac + a_vel

#paper 2008,automated code generation for DG methods
#a += a_extra

# Linear form
L = v*f*dx

# Set up boundary condition (apply strong BCs)
g = Expression("350", degree=fe_degree)
bc = DirichletBC(V_dg, g, bpin, "geometric")

# Solution function
phi_h = Function(V_dg)

# Assemble and apply boundary conditions
A = assemble(a)
b = assemble(L)
bc.apply(A, b)
if using_3D:
bc_d = DirichletBC(V_dg, 300, bdisc, "geometric")
bc_d.apply(A, b)

print(type(A))
#M = A.as_backend_type()
print(A)
print(b.array())

# Solve system
solve(A, phi_h.vector(), b)
print(phi_h.vector().array())

# Project solution to a continuous function space
up = project(phi_h, V=V_cg)

file = File("temperature.pvd")
file << up

# Plot solution
plot(up)
interactive()
#import matplotlib.pyplot as plt
#plt.show()​