Thank you
I resumed the calculation with finite element calculations. For mathematical calculation, there is a double flow at the coordinate point (0,0). This is what the code does according to mathematical theory. For what I’m doing is adapting this code to the experimental, i.e. passing the same electric wire on both edges as well as in the corner, so that there is no the double flow at point (0,0)?
How to rewrite the boundary conditions of flow to insert it into the variational formulation?
from __future__ import print_function
from dolfin import *
from fenics import *
#from ufl import *
from dolfin import cpp
import numpy as np
from numpy import cos, sin, pi, exp
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
#"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
#Constants (atous en SI unités)
sf = 1.0e-3
l = 1
numElems =1
numNoeudx= numElems+1
numNoeudy= numElems+1
nt=(numNoeudx)*(numNoeudy)-1
mesh = RectangleMesh(Point(0,0),Point(l,l), numElems, numElems,"left")
V = FunctionSpace(mesh, "Lagrange", 1)
class Left_bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and ((abs(x[0] - 0.) < DOLFIN_EPS)) or (abs(x[1] - 0) < DOLFIN_EPS)
vertical = Left_bottom()
tdim = mesh.topology().dim()
#print(tdim)
boundaries = MeshFunction('size_t', mesh, tdim-1)
boundaries.set_all(0)
vertical.mark(boundaries, 1)
File("bord/marker.pvd") << boundaries
u_0 = Expression("20",degree=1)
f = Expression("0",degree=1)
# Equation coefficients
File("marker.pvd") << boundaries
K = as_matrix([[ 0.3,0.],[0.,0.2]])
Q_heat = 0.4e3
g = Expression('Q_heat', degree=1, Q_heat=Constant(Q_heat)) #
#g= Constant(0.4e3)
rhoC = Constant(1.46555e6) # constant J/m³/°C
dsB = Measure("ds", subdomain_id=1, subdomain_data=boundaries) #
#print(coor.shape[0] )
def operator(u, v):
return ( inner(grad(v), K*grad(u)) ) #- f*v)*dx # - g*v*dsB #g*v*dsC #
u = TrialFunction(V)
v = TestFunction(V)
u0 = Function(V)
t_end =36
#num_steps =10 # nombres de pas de temps( number of time steps)
#dt = t_end/ num_steps
dt = 4
theta = Constant(2.0/3.0)
#F = rhoC*inner((u-u0), v)*dx +dt*theta*operator(u, v) + dt*(1.0-theta)*operator(u0, v) -dt*theta*fe -dt*(1.0-theta)*fe #
File("bord/marker.pvd") << boundaries
ke = inner(grad(v),K*grad(u))*dx
ce = rhoC*inner(u,v)*dx
fe = dt*g*v*dsB
C = assemble(ce)
KE = assemble(ke)
a = ce + dt*theta*ke
L = rhoC*inner(u0, v)*dx - dt*(1.0-theta)*inner(grad(v),K*grad(u0))*dx +fe
u = Function(V)
problem = LinearVariationalProblem(a, L, u)
solver = LinearVariationalSolver(problem) #
solver.parameters['linear_solver'] = 'cg'
solver.parameters['preconditioner'] = 'ilu'
A = assemble(a) #
#print (A.array())
b = None
u0.interpolate(u_0) #
Flux = assemble(fe)
Flux[nt]/=2
print(Flux.get_local())
vtkfile = File('janvier01/solution.pvd')
u.rename("Temerature", "temperature")
u.interpolate(u_0) #
t = 0.0
vtkfile << (u, t)
td = 0.0 # le compteur
while t <= t_end: #tant que ( boucle while)
b = (assemble(rhoC*inner(u0, v)*dx - dt*(1.0-theta)*inner(grad(v),K*grad(u0))*dx) + Flux )
#b = assemble(L, tensor=b)
A= assemble(a)
solver.solve()
vtkfile << (u, t)
u0.assign(u)
t += dt
td += 1
#L = rhs(F)
print(u.compute_vertex_values(mesh) ,t)
plt.figure(1)
P = plot(u,title='temperature à t = %g' % t)
plt.colorbar(P)
#plt.show()
plt.figure(2)
p = plot(u,mode ="color", title='temperature à t = %g' % t)
plt.colorbar(p)
plot(mesh)
plt.show()
j’ai assemblé le flux seul comme vous l’avez dit et l’intégrant dans la formulation variationnelle mais ça ne fonctionne pas, ça n’agit pas dans le calcul.
flux = dt*g*v*dsB
Flux = assemble(fe)
Flux[nt]/=2
je l’ai mis dans la forme ci dessous
b = (assemble(rhoC*inner(u0, v)*dx - dt*(1.0-theta)*inner(grad(v),K*grad(u0))*dx) + Flux )
How to write it so that the variational formulation for the solution?
