Skip to content
GitLab
Commit c4e1cc8a authored by p7k's avatar p7k
Browse files

Partial implementation of physics settings for Synchronous_Reluctance_Machine test

parent 13f94c7d
Hide whitespace changes
Inline Side-by-side
build_gallery.py
View file @ c4e1cc8a
......@@ -74,14 +74,14 @@ def plot_vector(fname):
xc = (nodes[tris[:,0],0]+nodes[tris[:,1],0]+nodes[tris[:,2],0])/3.0
yc = (nodes[tris[:,0],1]+nodes[tris[:,1],1]+nodes[tris[:,2],1])/3.0
plt.figure()
plt.set_cmap('inferno')
plt.gca().set_aspect('equal')
plt.tripcolor(nodes[:,0],nodes[:,1],tris,facecolors = mag,edgecolors='k')
#plt.tripcolor(nodes[:,0],nodes[:,1],tris,facecolors = mag,edgecolors='k')
plt.tripcolor(nodes[:,0],nodes[:,1],tris,facecolors = mag)
plt.colorbar()
#plt.grid(True)
#plt.quiver(xc,yc,vals[:,0],vals[:,1],mag,pivot='mid',edgecolors='k',linewidths=0.1)
plt.grid(True)
plt.axes().set_aspect('equal','datalim')
for ext in imgtypes:
plt.savefig(fname+ext)
......@@ -93,13 +93,13 @@ root = os.getcwd()
for path, _, files in os.walk(root):
for name in files:
name, fext = os.path.splitext(name)
if len(sys.argv) == 1 or name.startswith(sys.argv[1]):
fpath = os.path.join(path, name)
if fext == '.oesk':
plot_sketch(fpath)
elif fext == '.oeme':
plot_mesh(fpath)
elif fext == '.oesc':
plot_scalar(fpath)
elif fext == '.oeve':
plot_vector(fpath)
if len(sys.argv) == 1 or name.startswith(sys.argv[1]):
fpath = os.path.join(path, name)
if fext == '.oesk':
plot_sketch(fpath)
elif fext == '.oeme':
plot_mesh(fpath)
elif fext == '.oesc':
plot_scalar(fpath)
elif fext == '.oeve':
plot_vector(fpath)
src/Mesh/src/Mesh.cpp
View file @ c4e1cc8a
......@@ -701,9 +701,10 @@ void Mesh::get_triangles() {
/*
* Parses the mesh for darts representing unique triangles
*
* // TODO: Instead of starting from scratch for each iteration, sweep and mark existing triangles, then parse the remainder
*/
// TODO: Instead of starting from scratch for each iteration, sweep and mark existing triangles, then parse the remainder
Triangles.resize(0);
Triangles.reserve(2 * num_points());
......@@ -944,24 +945,18 @@ void Mesh::save_as(std::string path, std::string file_name) const {
std::fstream fs;
fs.open(path + file_name + ".oeme", std::fstream::out);
for (size_t e = 0; e != Edges.size(); ++e) {
for (auto const &t : Triangles) {
size_t e = t.edge();
Point const v0 = base(e);
Point const v1 = base(next(e));
Point const v2 = base(next(next(e)));
fs << v0.X << ',' << v1.X << ',' << v2.X << ',' << v0.Y << ',' << v1.Y << ',' << v2.Y << '\n';
}
/*
for (auto t : Triangles) {
size_t e{t.edge()};
e = next(e);
Point const v1 = base(e);
Point const v0 = base(e);
Point const v1 = base(next(e));
Point const v2 = base(next(next(e)));
e = next(e);
Point const v2 = base(e);
fs << v0.X << ',' << v1.X << ',' << v2.X << ',' << v0.Y << ',' << v1.Y << ',' << v2.Y << '\n';
}
*/
fs.close();
}
......@@ -1787,10 +1782,14 @@ std::shared_ptr<BoundaryConstraint> Mesh::boundary_constraint(std::shared_ptr<Cu
}
std::shared_ptr<Contour const> Mesh::select_contour(Point const p) {
/*
* Selects boundary contour containing Point p
*/
size_t e{0};
LocateTriangleResult result = locate_triangle(p, e);
if (result != LocateTriangleResult::Interior) {
std::cerr << "Point was not found in the interior of the mesh. Returned Contour is null." << std::endl;
return std::make_shared<Contour const>();
}
......
src/Mesh/src/Mesh.h
View file @ c4e1cc8a
......@@ -164,9 +164,9 @@ public:
/*
* Constructs a MappedBoundaryPair mesh constraint from a ContinuousBoundaryPair
*
* //TODO: Present implementation strategy allows cyclical mapping which should converge for boundaries restricted to power-of-2 parametric discretizations, but is untested
*/
//TODO: Present implementation strategy allows cyclical mapping which should converge for boundaries restricted to power-of-2 parametric disc
MappedBoundaries.push_back(std::make_shared<MappedBoundaryPair>(*this, cbp));
return MappedBoundaries.back();
}
......
src/Model Templates/src/DistributedWindingStator.cpp
View file @ c4e1cc8a
#include "DistributedWindingStator.h"
void DistributedWindingStator::add_to_sketch(Sketch &s, std::shared_ptr<Vertex> origin) {
// Dimensions are approximate, sketch.solve() + constraints will impose real values
double_t at = M_PI / TotalTeeth;
double_t ri = InnerRadius;
......@@ -91,4 +90,20 @@ void DistributedWindingStator::add_to_sketch(Sketch &s, std::shared_ptr<Vertex>
if (AddBoundingBox) {
add_bounding_box(s, origin, V[1], ro, ModeledTeeth * (360.0 / TotalTeeth));
}
}
std::vector<std::shared_ptr<const Contour>> DistributedWindingStator::select_slot_contours(Mesh &mesh) const {
std::vector<std::shared_ptr<const Contour>> slot_contours;
double_t r{slot_interior_radius()};
double_t a{M_PI / TotalTeeth};
double_t da{2 * a};
for (size_t i = 0; i != ModeledTeeth; ++i) {
Point p{r * std::cos(a),r * std::sin(a)};
slot_contours.push_back(mesh.select_contour(p));
a += da;
}
return slot_contours;
}
\ No newline at end of file
src/Model Templates/src/DistributedWindingStator.h
View file @ c4e1cc8a
......@@ -2,9 +2,18 @@
#define OERSTED_DISTRIBUTEDWINDINGSTATOR_H
#include "PolarModelTemplate.h"
#include "Mesh.hpp"
class DistributedWindingStator : public PolarModelTemplate {
public:
double_t inner_radius() const override { return InnerRadius; };
double_t outer_radius() const override {return InnerRadius + ToothFaceThickness + SlotDepth + BackironThickness; };
double_t slot_interior_radius() const { return InnerRadius + ToothFaceThickness + SlotDepth - SlotWidth / 2.0; };
std::vector<std::shared_ptr<const Contour>> select_slot_contours(Mesh &m) const;
void add_to_sketch(Sketch &s, std::shared_ptr<Vertex> origin) override;
size_t TotalTeeth;
......
src/Model Templates/src/PolarModelTemplate.cpp
View file @ c4e1cc8a
......@@ -7,6 +7,8 @@ void PolarModelTemplate::convexify(Sketch &s, std::shared_ptr<Vertex> origin, st
}
void PolarModelTemplate::add_bounding_box(Sketch &s, std::shared_ptr<Vertex> origin, std::shared_ptr<Vertex> vp00, double_t radius, double_t angle) {
std::cerr << "Implement bounding box scaling option" << std::endl;
double_t rb = radius * std::min(2.0, std::abs(sqrt(8.0 / (1.0 + std::cos(angle * M_PI / 180.0))) - 1.0));
auto vp10 = s.select_periodic_vertex(vp00, origin, angle);
......@@ -14,12 +16,12 @@ void PolarModelTemplate::add_bounding_box(Sketch &s, std::shared_ptr<Vertex> ori
auto vp01 = s.new_element<Vertex>(rb, 0.0);
auto vp11 = s.new_element<Vertex>(rb * std::cos(angle * M_PI / 180.0), rb * std::sin(angle * 180.0 / M_PI));
auto l0 = s.new_element<LineSegment>(vp00, vp01);
auto l1 = s.new_element<LineSegment>(vp10, vp11);
auto co = s.new_element<CircularArc>(vp01, vp11, origin, rb);
BoundingBoxLineSegment[0] = s.new_element<LineSegment>(vp00, vp01);
BoundingBoxLineSegment[1] = s.new_element<LineSegment>(vp10, vp11);
BoundingBoxCircularArc = s.new_element<CircularArc>(vp01, vp11, origin, rb);
s.new_element<Radius>(co, rb);
s.new_element<Horizontal>(l0);
s.new_element<Angle>(l0, l1, angle);
s.new_element<Coincident<LineSegment>>(origin, l1);
s.new_element<Radius>(BoundingBoxCircularArc, rb);
s.new_element<Horizontal>(BoundingBoxLineSegment[0]);
s.new_element<Angle>(BoundingBoxLineSegment[0], BoundingBoxLineSegment[1], angle);
s.new_element<Coincident<LineSegment>>(origin, BoundingBoxLineSegment[1]);
}
\ No newline at end of file
src/Model Templates/src/PolarModelTemplate.h
View file @ c4e1cc8a
......@@ -5,11 +5,19 @@
class PolarModelTemplate : public ModelTemplate {
public:
virtual double_t inner_radius() const = 0;
virtual double_t outer_radius() const = 0;
void convexify(Sketch &s, std::shared_ptr<Vertex> origin, std::shared_ptr<Vertex> vp, double_t angle);
void add_bounding_box(Sketch &s, std::shared_ptr<Vertex> origin, std::shared_ptr<Vertex> vp, double_t radius, double_t angle);
std::shared_ptr<CircularArc> bounding_box_circular_arc() const { return BoundingBoxCircularArc; };
protected:
std::array<std::shared_ptr<LineSegment>,2> BoundingBoxLineSegment;
std::shared_ptr<CircularArc> BoundingBoxCircularArc;
};
#endif //OERSTED_POLARMODELTEMPLATE_H
\ No newline at end of file
src/Model Templates/src/SynchronousReluctanceRotor.h
View file @ c4e1cc8a
......@@ -5,6 +5,10 @@
class SynchronousReluctanceRotor : public PolarModelTemplate {
public:
double_t inner_radius() const override { return InnerRadius; };
double_t outer_radius() const override { return OuterRadius; };
void add_to_sketch(Sketch &s, std::shared_ptr<Vertex> origin) override;
size_t Poles;
......
test/LibraryIntegration/Mesh_To_FEM_Test.cpp
View file @ c4e1cc8a
......@@ -1010,16 +1010,7 @@ public:
me = Mesh{sk};
// Periodic boundaries
for (auto &pb : periodic_boundary) {
me.add_mapped_boundary_pair(pb);
/*
auto bc0 = me.boundary_constraint(pb.curve0());
bc0->uniform_discretization(true);
auto bc1 = me.boundary_constraint(pb.curve1());
bc1->uniform_discretization(true);
*/
}
me.add_mapped_boundary_pair(periodic_boundary);
// Airgap boundary
auto discrete_airgap = me.boundary_constraint(continuous_airgap);
......
test/Model Templates/test_ModelTemplates.hpp
View file @ c4e1cc8a
......@@ -4,6 +4,7 @@
#include "Sketch.hpp"
#include "Mesh.hpp"
#include "ModelTemplates.hpp"
#include "Physics.hpp"
#include "gtest.h"
......
test/Model Templates/test_SynchronousReluctanceMachine.cpp
View file @ c4e1cc8a
......@@ -6,9 +6,13 @@ public:
}
};
TEST_F(Synchronous_Reluctance_Machine, geometry_test) {
TEST_F(Synchronous_Reluctance_Machine, test) {
// Create sketch
size_t Np = 8;
size_t Nt = 48;
SynchronousReluctanceRotor srr;
srr.Poles = 8;
srr.Poles = Np;
srr.InnerRadius = 25.4e-3;
srr.OuterRadius = 89e-3-0.03*25.4e-3;
......@@ -30,8 +34,8 @@ TEST_F(Synchronous_Reluctance_Machine, geometry_test) {
dws.SlotWidth = 6.05e-3;
dws.SlotOpening = 1.74e-3;
dws.TotalTeeth = 48;
dws.ModeledTeeth = 6;
dws.TotalTeeth = Nt;
dws.ModeledTeeth = Nt / Np;
dws.Convexify = false;
dws.AddBoundingBox = true;
......@@ -44,15 +48,24 @@ TEST_F(Synchronous_Reluctance_Machine, geometry_test) {
srr.add_to_sketch(sketch,origin);
dws.add_to_sketch(sketch,origin);
CylindricalAirgap<2> cag{sketch,origin,srr.OuterRadius,0.0,dws.InnerRadius,0.0,360.0/8.0};
CylindricalAirgap<2> airgap{sketch,origin,srr.OuterRadius,0.0,dws.InnerRadius,0.0,360.0/Np};
EXPECT_LE(sketch.solve(), 100e-3 * FLT_EPSILON);
EXPECT_TRUE(sketch.build());
sketch.save_as<SaveMethod::Rasterize>(SAVE_DIR,"Synchronous_Reluctance_Machine");
// Sketch selections for boundary conditions
auto airgap_motion_interface = airgap.interface(1);
auto periodic_boundaries = sketch.select_periodic_boundary_pairs(origin, 360.0/Np);
std::shared_ptr<Curve const> magnetic_insulation_boundary = dws.bounding_box_circular_arc();
// Create mesh
Mesh mesh{sketch};
mesh.boundary_constraint(airgap_motion_interface)->uniform_discretization(true);
mesh.add_mapped_boundary_pair(periodic_boundaries);
mesh.create();
EXPECT_TRUE(mesh.edges_are_valid());
......@@ -67,6 +80,92 @@ TEST_F(Synchronous_Reluctance_Machine, geometry_test) {
mesh.save_as(SAVE_DIR,"Synchronous_Reluctance_Machine");
// Mesh selections for material properties and forcing functions
double_t r = srr.inner_radius() * (1 + FLT_EPSILON);
double_t x = r * std::cos(M_PI/Np);
double_t y = r * std::sin(M_PI/Np);
auto rotor_iron = mesh.select_contour(Point{x,y});
EXPECT_TRUE(rotor_iron);
r = dws.outer_radius() * (1 - FLT_EPSILON);
x = r * std::cos(M_PI/Np);
y = r * std::sin(M_PI/Np);
auto stator_iron = mesh.select_contour(Point{x,y});
EXPECT_TRUE(stator_iron);
auto slot_contours = dws.select_slot_contours(mesh);
for (auto sc : slot_contours) {
EXPECT_TRUE(sc);
}
// Setup FEA simulation
auto fem = std::make_shared<FiniteElementMesh<1, 1>>(mesh);
auto msp = std::make_shared<Magnetostatic>(fem);
// Set material properties
MaterialProperties electrical_steel = JFE_35JN210();
//MaterialProperties electrical_steel = MaterialProperties(std::make_shared<LinearIsotropicMagneticMaterial>(1000.0));
fem->region(rotor_iron)->material(electrical_steel);
fem->region(stator_iron)->material(electrical_steel);
// Conditions
FunctionArguments fargs;
fargs.add_argument("t",1.0/8.0);
fargs.add_argument("J",M_SQRT2 * 30e6 * 0.5);
// Current Density
msp->add_current_density([](FunctionArguments args) { return -args["J"] * sin(2 * M_PI * args["t"] - M_PI * 2 / 3); }, {slot_contours[0],slot_contours[1]});
msp->add_current_density([](FunctionArguments args) { return +args["J"] * sin(2 * M_PI * args["t"] + M_PI * 0 / 3); }, {slot_contours[2],slot_contours[3]});
msp->add_current_density([](FunctionArguments args) { return -args["J"] * sin(2 * M_PI * args["t"] + M_PI * 2./ 3); }, {slot_contours[4],slot_contours[5]});
// Boundary Conditions
msp->add_magnetic_insulation({magnetic_insulation_boundary});
// Rotating Interface
auto position = msp->add_sliding_interface(airgap_motion_interface, true);
//double_t t = 0.0;
//double_t dt = 1.0 / (4.0 * position->size());
msp->add_periodic_boundary(periodic_boundaries, true);
// Post Processor
//PostProcessorInterface pp;
//pp.add<AirgapTorque>(msph,"Torque",fem->region(airgap_contour), rro*2.0/3.0 + rsi*1.0/3.0, rro*1.0/3.0 + rsi*2.0/3.0, 1.0, 4.0);
// Solve
//Oe::VectorXd v_prev;
msp->assemble();
//size_t step = 32;
//for (size_t iter = 0; iter != position->size() / step; ++iter) {
auto solution = msp->new_solution(); // TODO: What about warm start using old initial condition?
// if (iter > 1)
// solution->v = v_prev;
msp->solve(solution, fargs);
// v_prev = solution->v;
//Verify equation is solved
EXPECT_LE(solution->r.norm(), FLT_EPSILON * solution->f.norm() * sqrt(solution->f.size()));
// Save image
Eigen::VectorXd vv = msp->recover_boundary(solution->v);
fem->write_scalar(vv, SAVE_DIR, std::string("Synchronous_Reluctance_Machine_A"));
Eigen::ArrayXd Bx = msp->derivatives().dy(0).transpose() * solution->v;
Eigen::ArrayXd By = -msp->derivatives().dx(0).transpose() * solution->v;
fem->write_vector(Bx, By, SAVE_DIR, std::string("Synchronous_Reluctance_Machine_B"));
// Increment Position
// *position += step;
//t += 32 * dt;
//fargs["t"] = t;
// msph->assemble();
// std::cout << iter << "," << pp("Torque",solution->v) << std::endl;
//}
// TODO: Periodic Boundaries
// TODO: Rotational Boundaries
}
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment