CFD Online Discussion Forums - Computation of the approximate gradient

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- Main CFD Forum (https://www.cfd-online.com/Forums/main/)

- - Computation of the approximate gradient (https://www.cfd-online.com/Forums/main/139869-computation-approximate-gradient.html)

Computation of the approximate gradient

Hello,

I compute eikonal equation, which gives me a matrix of potential on the domain and I need to compute gradient of this potential. I use following approximation:

grad_i = 1/|C_i| * sum((phi_i + phi_j)/2 * normal_ij),

where C_i is finite volume, sum is over all neighbors of cell i, phi_i is potential in cell i and normal_ij is outward unit normal to surface between cell i and cell j. My implementation works fine on rectangular grid, but on triangular mesh it gives me bad results. Don't you know, where could be mistake on triangular grid?

My implementation in Matlab for triangular grid, where variable v are finite volumes, p are points and phi is potential:

Code:

for i=1:num_v

    neighbors = pdeent(v,i); %neighbors 

    neighbors(neighbours == i) = []; %remove cell i from neighbors 

    num_nbs = size(neighbours,2); %neighbors count

    if(num_nbs ~= 3) %on boundary zero

        continue;

    end

    

    for j=1:num_nbs

        nb_v = neighbours(j); %adjacent cell id

        edge_points = intersect(v(1:3,i),v(1:3,nb_v)); %edge points

        

        %edge size

        x = p(1,edge_points(1)) - p(1,edge_points(2));

        y = p(2,edge_points(1)) - p(2,edge_points(2));

        edges_h(1,j) = sqrt(x^2 + y^2);

        

        %computation of outward unit normal

        p3_i = v(1:3,i); p3_i(p3_i==edge_points(1)) = []; p3_i(p3_i==edge_points(2)) = [];

        normal = -flipud(p(:,edge_points(2)) - p(:,edge_points(1))).*[-1;1]; %normal

        e2_i =  p(:,p3_i) - p(:,edge_points(1)); %vector for control of outward

        if(dot(normal,e2_i) > 0) %if not outward

            normal = -normal; %switch

        end

        unit_normal(1:2,j) = normal/norm(normal); %unit normal



        average(j,1) = (phi(i) + phi(nb_v))/2;

    end

    

    gradient(1:2,i) = (unit_normal(1:2,:)*average)/max(edges_h);

end

Thank you.
Jakub

Your formulation is wrong for a Gauss gradient calculation. It's likely that the regular mesh that you used before hid the error.

The Gauss gradient formula turns the gradient averaged over the cell (a volume integral) into a surface integral of the function value over the boundary. Dividing by the cell volume is correct--that part comes from the volume integral. But the summation in the formula is approximating the surface integral. That means that the normal_ij in your formula is NOT a unit vector. It is the face area vector and should have a magnitude that is equal to the face's physical area.

I haven't debugged your code, but if you use the face area vector in place of the unit normal vectors, the formulation should be correct. Good luck.

PS. Have a look here: http://www.cfd-online.com/Wiki/Gradient_computation

Also, you are using a pretty sloppy interpolation for the face value in just taking the average of the two. Generally, it is better to at least use linear interpolation. In fact, even linear interpolation without correction may be inadequate if your mesh is skewed.