Commit 5563d5d3 authored by matthieu's avatar matthieu
Browse files

Energie cinétique du fluide

parent 55de2d13
......@@ -31,7 +31,7 @@ class Brinkman :
nu = self.mu/self.rhof
if self.c is not None :
val[..., _U_] = - solgrad[..., _p_, :]/self.rhof \
+ nu * (self.dot(solgrad[...,_U_,:], self.dca)-sol[...,_U_]/self.ca**2*self.dot(self.dca, self.dca))
- nu * (self.dot(solgrad[...,_U_,:], self.dca)-sol[...,_U_]/self.ca**2*self.dot(self.dca, self.dca))
val[..., _p_] = - self.c * (solgrad[..., _u_, 0] + solgrad[..., _v_, 1] + (0 if self.D == 2 else solgrad[..., _w_, 2]))
val[..., _p_] += - self.c * (self.dvs[..., 0, 0] + self.dvs[..., 1, 1] + (0 if self.D == 2 else self.dvs[..., 2, 2]))
else :
......@@ -82,10 +82,10 @@ class Brinkman :
F = solgrad[..., _p_, :]
#drag
GoU = f * Cd_u * self.rhof / 2 * Ap
reaction = F * self.pf.volume + (GoU*U) / (1 + self.stab/self.pf.mass * GoU)
reaction = (F * self.pf.volume + (GoU*U) / (1 + self.stab/self.pf.mass * GoU))*2
if fonpart is not None:
fonpart[...] = -reaction - self.g * self.pf.volume * self.rhof
val[..., _U_] = reaction / self.rhof
val[..., _U_] = reaction / self.rhof
val[..., _p_] = 0
return val
......@@ -100,6 +100,7 @@ class Brinkman :
self.pf = pf
self.g = np.array(g)
self.epsilon = epsilon
self.Ecf = 0
def particleInteractionTerm(self, meshJ, solution, res, jac) :
_p_ = self.D
......@@ -125,6 +126,7 @@ class Brinkman :
jac.add_to_field(f, g, pp[:, None, :] * ff[:, :, None], self.pf.eid)
def fluidTerm(self, meshJ, solution, res, jac) :
print("FLUIDTERM")
dofm = solution.dof_manager()
qp, qw = quadrature.get_triangle(3) if self.D == 2 else quadrature.get_tetrahedron(3)
self.vs = self.pf.momentum.at_qp(qp)
......@@ -148,6 +150,14 @@ class Brinkman :
dPsiW = basis.df(qp) * qw[:, None, None]
r1 = np.tensordot(dPsiW, x1[:, :, i, :], [[0, 2], [1, 2]])
np.add.at(res, dofm._fields[i][0], r0 + r1)
#test pour l'energie
_U_ = range(self.D)
solu = solution.at_qp(qp)
vL = (solu[...,_u_]**2 +solu[...,_v_]**2)* 0.5
Vg = meshJ.multJ(vL)
Ecf1 = np.tensordot(qw,Vg[:,:],[0,1]) * self.rhof
self.Ecf = sum(Ecf1)
if jac is None:
return
......@@ -171,6 +181,13 @@ class Brinkman :
if (np.max(np.abs(d11[:,:,f,:,g,:])) > 1e-8) :
dPsiDPsiW = np.einsum("qmx, qny, q -> qxymn", basisf.df(qp), basisg.df(qp), qw)
jac.add_to_field(f, g, np.tensordot(meshJ.multDXiDXLeft(meshJ.multDXiDXLeft(meshJ.multJ(d11[:,:,f,:,g,:]), 3), 2), dPsiDPsiW, 3))
#test pour l'energie
sol = solution.at_qp(qp)
vL = sol[...,_U_]**2 * 0.5
Vg = meshJ.multJ(vL)
self.Ecf = np.tensordot(qw,Vg[:,:],[0,1]) * self.rhof
print("ECFFFF",self.Ecf)
def massTerm(self, alpha, solution, res, jac) :
pass
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment