Hydrostatic pressure in PDE toolbox
5 次查看(过去 30 天)
显示 更早的评论
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
回答(1 个)
Anurag Ojha
2024-6-13
Hello Jigar
To apply varying hydrostatic pressure on the side face of the meshed rectangular box, you can follow these steps:
- Load the STL file using the stlread function. This function reads the vertices and faces of the mesh from the STL file.
- Calculate centeroid of each face
- Define the hydrostatic pressure as a function of the centeroid coordinates.
- 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
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!