Dear Simone Di Gregorio,
You are of course right saying that the approximation of transport
equation with a finite element method may induce some instabilities.
This is true except in certain cases where the coercivity is preserved
(for instance with a Dirichlet condition at the right end point and a
with monotone velocity field). Adding the viscous term solves the
problem of course for beta large enough as you explain.
Your formulation seems to be perfectly correct to me. The only thing is
that it seems to me that it is transposed. I would instead write
asm_real_or_complex_1_param(M,mim,mf,mf_data,V,rg, "a=data$1(#2);"
"M$1(#1,#1)+=comp(Grad(#1).Base(#1).Base(#2))(:,i,:,i,p).a(p);")
because there is a test function for each line of the linear system (not
for each column).
Of course, this is the expression for a 1D (straight) domain.
A suggestion : you should pass to the high level generic assembly which
is more easy to work with (and is now faster).
Best regards,
Yves
Le 16/03/2017 à 10:11, Simone Di Gregorio a écrit :
>
> Dear GetFem users,
>
>
> I'm dealing with a pure stationary advection problem in 1D
> straight single branch domain. So the equation is:
>
>
> /d/( v(s)c(s))//ds = 0 /
>
> /
> /
>
> /with /v(s) velocity field (NOT constant, so it varies with s) and
> c(s) solute concentration, while s is the arc length of the vessel.
> The boundary conditions are Dirichlet B.C. in inlet of the vessel
>
> /
> /
>
> The weak form is:
>
>
> (\partial_s \phi, v c)_{\Lambda} + v(\Gamma _{out})c(\Gamma
> _{out})\phi(\Gamma _{out}) = 0
>
>
> where \phi is the piecewise linear continuous test function, \Lambda
> si the domain, \Gamma _{out} is the outlet boundary.
>
>
> So I build the advection matrix with:
>
>
> asm_real_or_complex_1_param(M,mim,mf,mf_data,V,rg, "a=data$1(#2);"
> "M$1(#1,#1)+=comp(Base(#1).vGrad(#1).Base(#2))(:,:,i,j,p).a(p);")
>
>
> Then I subtract to the last element of the matrix (in position N,N)
> the value of v(\Gamma _{out}).
>
> As expected the solution is unstable, so i added an artificial
> diffusion in the following way:
>
>
> +D(\partial_s \phi, \partial_s c)_{\Lambda}
>
>
> with D= beta*max(v(s))*h, beta large enough and h dimension of the
> element, and imposing dc/ds=0 over the boundary. I assembled this with
>
> asm_real_or_complex_1_param(A,mim,mf,mf_data,D,rg, "a=data$1(#2);"
> "M$1(#1,#1)+=comp(vGrad(#1).vGrad(#1).Base(#2))(:,i,j,:,i,j,p).a(p);")
>
>
> The solution becomes stable.
>
>
> The problem is that:
>
> -the solution is always the same for very different function of v(s):
>
> -if D is not very high in the first elements of the vessel, the
> solution is constant, whereas it varies in the last element(this
> happens because add v(\Gamma _{out})
>
> otherwise the solution would be costant everywhere)
>
> -if D is quite high the solution is a simple exponential, but is not
> sensitive to variation of v
>
>
> From an analytical solution I expect that v(s)c(s) remain constant so
> c(s) should vary along the vessel dependently on value of v(s).
>
>
> So it's like the problem doesn't feel the variation of velocity and I
> can't understand what i'm doing wrong.
>
>
> Thanking you in advance, i hope i was clear enough in my explanation
>
>
> Best Regards
>
>
> Simone Di Gregorio
>
>
>
>
> _______________________________________________
> Getfem-users mailing list
> [email protected]
> https://mail.gna.org/listinfo/getfem-users
--
Yves Renard ([email protected]) tel : (33) 04.72.43.87.08
Pole de Mathematiques, INSA-Lyon fax : (33) 04.72.43.85.29
20, rue Albert Einstein
69621 Villeurbanne Cedex, FRANCE
http://math.univ-lyon1.fr/~renard
---------
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
