4

When dealing with conservation laws like your case, you can often make use of the divergence theorem (as you did). You can then express the fact that the total mass within your integration region is preserved by the following surface integral: $$\oint_{\partial \Omega} k \nabla T \cdot \mathbf{n} ~\partial S = 0$$ Now, as it stands, it is irrelevant which ...


3

The correct way is to average the source term which you can do easily in this case as you have a polynomial. In general you can do the average with a quadrature. For second order accuracy if the source term is smooth as in your case, you can also just evaluate it at the cell center, which is like mid-point quadrature. It should not matter much in your ...


2

Normal direction depends on the cell that you are writing equation for. the word outward is relative to the cell under study. In order to write equation for each of cells, i.e. $\Sigma \nabla T.n S_f=0$, stick to this : $\nabla T_{face}=\frac{T_c-T_i}{r_c-r_i}$ and assume $n$ as outward pointing normal vector for that face. I think your problem is that you ...


2

An extended answer. For more arbitrary meshes you have to consider that generally CFD/FEM solvers rely on generic data-structures with element and side lists: Element list Side list Consider the following pictures, which is the standard case for simple Cartesian meshes. Since there is a single plus and a single minus side on each face, the definition is ...


1

The smooth solution turned out to have BC's applied in the following way: Walls and inlet: $\frac{\partial p}{\partial n}=0$ Outlet: $p=0$


1

This is speculation, as I do not know the Lipschitz constant or the derivative scales of your simulation. Also, there might be some resonance effect in the interplay of the discrete and continuous parts. But what I would first draw attention to is that the error of RK4 (and any other method) has a V shape in a loglog plot vs. the step size. This is the ...


1

I believe that the issue you are facing emanates from the type of triangular mesh you are using. This particular discretisation has in-built anisotropy; note the alignment of all of the longest edges is parallel to one of the diagonals of the square. You will observe a different behaviour in the results if you choose the alignment parallel to the other ...


1

Here is a fast implementation using sparse matrices and sparse Jacobian estimation: % define square domain [-1,1] x [-1,1] n = 51; x=linspace(-1,1,n); y=x; [X,Y]=meshgrid(x,x); % build finite differences operators dx=x(2)-x(1); e=ones(n,1); d0x=ones(n,1); grad = spdiags([-e 0*e e],-1:1,n,n)/2/dx; % use Kronecker product to build matrix of d/dx and d/dy ...


Only top voted, non community-wiki answers of a minimum length are eligible