Transient in Mixed formulation


109
views
0
6 weeks ago by
Hello, I am not able to solve transient one way coupling, of 1*1 domain where i have applied dritchlet boundary condition in both the Weak forms ,  i want to project initial temperature of the domain as zero.   want variation of temperature at particular point answer coming as [nan nan nan]

from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np

T = 5.0 # final time
num_steps = 50 # number of time steps
dt = T / num_steps # time step size
mu = 0.001 # dynamic viscosity
rho = 1 # density
alpha = 0.001 #Cofficient of thermal Expansion
beta = 1.25
lambda_ = beta

# Create mesh
mesh = UnitSquareMesh(8, 8)

# Define function spaces
P1 = VectorElement("CG", mesh.ufl_cell(), 2) #displacement u
P2 = FiniteElement("CG", mesh.ufl_cell(), 1) #Temperature
a = MixedElement([P1,P2])
TH = FunctionSpace(mesh, a)

# Define boundaries
def clamped_boundary(x, on_boundary):
return on_boundary and near(x[0],0)

def right_boundary(x, on_boundary):
return on_boundary and near(x[0],1)

# Define boundary condition For Displacement
bc_leftu = DirichletBC(TH.sub(0),Constant((0,0)),clamped_boundary)

bc_rightu = DirichletBC(TH.sub(0).sub(0), Constant(0.5),right_boundary)

# Define boundary condition For Temperature

bc_leftt = DirichletBC(TH.sub(1),Constant((0.0)),clamped_boundary)

bc_rightt = DirichletBC(TH.sub(1),Constant((100.0)), right_boundary)


bc = [bc_leftu,bc_rightu, bc_leftt, bc_rightt]

# Define trial and test functions
(u,p) = TrialFunctions(TH)
(v,q) = TestFunctions(TH)

# Define functions for solutions at previous and current time steps
th =Function(TH)
th0 =Function(TH)

u , t = split(th)
u0, t0 = split(th0)

#t0.sub(1)=project(Constant((0.0),th0.sub(1))
#t0 = project(Constant(0.0), t0)

# Define expressions used in variational forms
alpha= Constant(alpha)
k = Constant(dt)



# Define symmetric gradient
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)

# Define stress tensor
def sigma(u,p):
return lambda_*nabla_div(u)*Identity(2) + 2*mu*epsilon(u) - p*alpha*Identity(2)

# Define variational problem
F1 = inner(sigma(u,p), epsilon(v))*dx
F2 = inner(grad(p),grad(q))*dx + dot((p - t0) / k, q)*dx
F= F1+F2
a = lhs(F)
L = rhs(F)

# Time-stepping
pta = Point(0.5,0.5)

file = File("transient.pvd")
for n in range(num_steps):
solve(a == L, th, bc)
th0.assign(th)
print (th(pta))
file << th
interactive()


SOLUTION-
Solving linear variational problem.
[ nan nan nan]
Solving linear variational problem.
[ nan nan nan]
Solving linear variational problem.
[ nan nan nan]
Community: FEniCS Project

1 Answer


3
6 weeks ago by
The problem is the order in which you specify different 'functions', it should go as follows:
# Define trial and test functions
(u, p) = TrialFunctions(TH)
(v, q) = TestFunctions(TH)

# Define functions for solutions at previous and current time steps
th0 = Function(TH)
(u0, t0) = split(th0)

# Define symmetric gradient
...

# Define stress tensor
...

# Define variational problem
...

th = Function(TH)

# Time-stepping
for n in range(num_steps):
    solve(a == L, th, bc)
    (u, t) = th.split()
    th0.assign(th)​
Thanks Adam for your valuable time,

how can i impose/project initial values of temperature, i have tried by 2 codes below, its not working

#t0.sub(1)=project(Constant((0.0),th0.sub(1))
#t0 = project(Constant(0.0), t0)

these above codes work in previous version of fenics , i m using 2017.1

please have a glance over it, thanks
written 5 weeks ago by Md Ayub Ansari  
1
You have to project onto function spaces not functions, i.e. something like
t0 = project(Constant(0.0),TH.sub(1).collapse()) should work.
written 5 weeks ago by Lukas O.  
1
Actually, there is no need to impose zero initial condition, th0 =Function(TH) sets all DOFs to zero. For non-zero conditions, just follow the tutorial, p. 80: https://fenicsproject.org/pub/tutorial/pdf/fenics-tutorial-vol1.pdf
written 5 weeks ago by Adam Janecka  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »