PDE Coefficients in the pde toolbox
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Hi,
I'm going to solve the following elliptic equation in a domain which consists of two materials:
According to pde toolbox syntax the coefficients of the above (nonlinear) elliptic equation are:
c = K a = 0 f = d(K)/dz = d(K)/dh * dh/dz
but when I plot the flux terms (inside square brackets terms) using:
uintrp = pdeintrp(p,t,u);
[ux,uz] = pdegrad(p,t,u) ;
for j=1:2
idx = Sidx==j ; % 'Sidx' contains the material index (1 or 2) for each triangle
qx(idx) = -F(j).K(uintrp(idx)) .* ux(idx) ;
qz(idx) = -F(j).K(uintrp(idx)) .* (uz(idx)+1) ;
end
figure ; pdeplot(p,e,t,'flowdata',[qx; qz])
which should be conserved in the entire domain, I get this:
By setting "f=0" which is equivalent for:
the result is completely acceptable (the fluxes are equal everywhere). So it seems something is wrong with "f". For the reference "f" is calculated by:
function f=fCoef(p,t,u,time)
numElems = size(t,2) ;
f = zeros(1,numElems) ;
uintrp = pdeintrp(p,t,u); % Interpolated values at centroids
[~,uz] = pdegrad(p,t,u); % Approximate derivatives
for j=1:numElems
z1 = p(2,t(1,j)) ;
z2 = p(2,t(2,j)) ;
z3 = p(2,t(3,j)) ;
zm = (z1+z2+z3)/3 ;
if zm<=4
f(j) = F(2).dKdh(uintrp(j)).* uz(j) ;
else
f(j) = F(1).dKdh(uintrp(j)).* uz(j) ;
end
end
Where am I missing something?
Thanks
0 comentarios
Respuestas (0)
Ver también
Categorías
Más información sobre Boundary Conditions en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!