Hydrostatic pressure in PDE toolbox

Hi,

I have this rectangular box in stl file which i have imported. I did mesh it and now i want to apply varying hydrostatic pressure on the side face, how do i do that.

Thanks,

Jigar

[R]

The function structureboundaryload is useful to specify boundary loads for structural models. This example might help:

https://www.mathworks.com/help/pde/ug/pde.structuralmodel.structuralboundaryload.html#mw_971d0107-79cf-453d-8811-240e2667116f:~:text=Pressure%22%2CpLoad)%3B-,Specify%20Nonconstant%20Pressure%20For%20Transient%20Model%20by%20Using%20Function%20Handle,-Apply%20Sinusoidal%20Pressure

[Jigar]

Hi @R, thank you so much for the help. It helped. This are for transient and for varying load/pressure on full face or a point. I am trying to understand how i can give it to a face where the pressure will be higher at the bottom and almost non at the top of the side wall.

I am trying to understand how i can provide H (height) in an equation to wall pressure to vary pressure on the surface as the hydrostatic pressure will be rho times g times H.

[Jigar]

Hi guys, I am still looking for anything that might help here. Has someone tried before to provide uniformly varyling load as well?

Answer 1

[Anurag Ojha]

Hello Jigar

To apply varying hydrostatic pressure on the side face of the meshed rectangular box, you can follow these steps:

1. Load the STL file using the stlread function. This function reads the vertices and faces of the mesh from the STL file.
2. Calculate centeroid of each face
3. Define the hydrostatic pressure as a function of the centeroid coordinates.
4. Iterate over each face and calculate the pressure at its centroid. Apply the pressure as a force vector to the face.

I am adding code for your reference . I have taken mock data and certain assumptions to write the below code. Kindly modify it as per your use case:

% Mock data for a simple rectangular box
% Vertices of the box
vertices = [0 0 0; 1 0 0; 1 1 0; 0 1 0; 0 0 1; 1 0 1; 1 1 1; 0 1 1];

% Faces of the box (triangular faces)
faces = [
    1 2 3; 1 3 4; % Bottom face
    5 6 7; 5 7 8; % Top face
    1 2 6; 1 6 5; % Front face
    2 3 7; 2 7 6; % Right face
    3 4 8; 3 8 7; % Back face
    4 1 5; 4 5 8  % Left face
];

% Create triangulation object
TR = triangulation(faces, vertices);

% Save the triangulation to an STL file
stlwrite(TR, 'box.stl');

% Load the STL file
[vertices, faces] = customStlRead('box.stl');

% Calculate centroids of each face
centroids = (vertices(faces(:, 1), :) + vertices(faces(:, 2), :) + vertices(faces(:, 3), :)) / 3;

% Define the hydrostatic pressure as a function of the centroid coordinates
pressure = @(z) 1000 * z; % Example linear pressure function

% Initialize forces matrix
forces = zeros(size(faces, 1), 3);

% Iterate over each face and calculate the pressure at its centroid
for i = 1:size(faces, 1)
    centroid = centroids(i, :);
    z = centroid(3); % Height coordinate of the centroid
    p = pressure(z); % Calculate pressure at the centroid
    
    % Calculate face normal vector
    v1 = vertices(faces(i, 2), :) - vertices(faces(i, 1), :);
    v2 = vertices(faces(i, 3), :) - vertices(faces(i, 1), :);
    normal = cross(v1, v2);
    normal = normal / norm(normal); % Normalize the vector
    
    % Calculate force vector (pressure times area times normal vector)
    area = norm(cross(v1, v2)) / 2; % Area of the triangle
    force = p * area * normal; % Force vector
    forces(i, :) = force;
end

% Visualize the forces on the mesh
figure;
patch('Faces', faces, 'Vertices', vertices, 'FaceColor', 'cyan', 'EdgeColor', 'black');
hold on;
quiver3(centroids(:, 1), centroids(:, 2), centroids(:, 3), forces(:, 1), forces(:, 2), forces(:, 3), 'r');
hold off;
xlabel('X');
ylabel('Y');
zlabel('Z');
title('Hydrostatic Pressure on Side Face of Rectangular Box');

% Custom STL read function
function [vertices, faces] = customStlRead(filename)
    % Open the STL file in binary mode
    fid = fopen(filename, 'r');
    if fid == -1
        error('File could not be opened, check name or path.')
    end

    % Read file header
    fread(fid, 80, 'uchar=>schar');
    numFaces = fread(fid, 1, 'int32');

    % Initialize arrays to hold the data
    faces = zeros(numFaces, 3);
    vertices = zeros(3 * numFaces, 3);

    % Read the data for each face
    for i = 1:numFaces
        fread(fid, 3, 'float32'); % Normal vector
        vertices(3 * i - 2, :) = fread(fid, 3, 'float32'); % Vertex 1
        vertices(3 * i - 1, :) = fread(fid, 3, 'float32'); % Vertex 2
        vertices(3 * i, :) = fread(fid, 3, 'float32'); % Vertex 3
        fread(fid, 2, 'uchar'); % Attribute byte count
        faces(i, :) = 3 * i - 2:3 * i;
    end

    % Close the file
    fclose(fid);

    % Remove duplicate vertices
    [vertices, ~, J] = unique(vertices, 'rows');
    faces = J(faces);
end