### Transient in Mixed formulation

383

views

0

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]

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

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)
```

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.

t0 = project(Constant(0.0),TH.sub(1).collapse()) should work.

written
9 months 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
9 months ago by
Adam Janecka

Please login to add an answer/comment or follow this question.

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