Hydrostatic pressure in PDE toolbox

1 次查看（过去 30 天）
Jigar
Jigar 2024-5-17
回答： Anurag Ojha 2024-6-13
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
  3 个评论
显示 1更早的评论
Jigar
Jigar 2024-5-17
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
Jigar 2024-5-26
Hi guys, I am still looking for anything that might help here. Has someone tried before to provide uniformly varyling load as well?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Anurag Ojha
Anurag Ojha 2024-6-13
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Geometry and Mesh 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by