Element Length
Contents
A test MATLAB program for computing the "element length" as defined by Equation 4.11 of Tezduyar (1992, p. 16):
\[ h = 2 \left( \sum_{a=1}^{n_{en}} | s \cdot \nabla N_a | \right)^{-1} \].
function h = elem_length
Element definition
Defines the nodal locations, arbitrary velocity vectors at the nodes, and the element shape functions for a 4-node quadralateral element.
% Define node locations (local) px = [-1,1,1,-1]; py = [-1,-1,1,1]; % Define nodal velocies v(1,:) = [1,2]; v(2,:) = [1,0]; v(3,:) = [1,1]; v(4,:) = [0,0]; % Define element shape functions as MATLAB symbolic expressions syms x y; N_sym{1}(x,y) = 1 / 4 * (1 - x)*(1 - y); N_sym{2}(x,y) = 1 / 4 * (1 - x)*(1 + y); N_sym{3}(x,y) = 1 / 4 * (1 + x)*(1 + y); N_sym{4}(x,y) = 1 / 4 * (1 + x)*(1 - y);
Compute summation components of the element length
Creates a vector of length \(n_{en}\) that contains the value to be summed: \(\left|\vec{s} \cdot \nabla N_a\right|\). If the velocity is zero then so is the associated value.
for i = 1:length(N_sym); % Compute the norm of the velocity vector, if it is zero set the value % to zero and progress to the next node n = norm(v(i,:)); if n == 0; h(i) = 0; continue; end % Compute the normal vector from the velocity s = v(i,:)/ norm(v(i,:)'); % Gradient of shape function at current node f{1} = matlabFunction(diff(N_sym{i},x)); f{2} = matlabFunction(diff(N_sym{i},y)); % Evaluate gradient at current node for j = 1:length(f); B(j) = f{j}(px(i),py(i)); end % Compute the summation component hvec(i) = abs(dot(s,B)); end
Compute the scalar element length, this is the value returned
h = 2*sum(hvec)^(-1);
References
