From 17942cdbe5b8f23f5f37e461f3c0e3ce24c0f3e6 Mon Sep 17 00:00:00 2001 From: Akbar Rahman Date: Tue, 26 Dec 2023 23:37:47 +0000 Subject: [PATCH] calculate 2d beam element stiffness matrix --- src/bin/two_d.rs | 70 +++++++++++++++++ src/lib.rs | 7 ++ src/two_d.rs | 14 ++++ src/two_d/beam_elements.rs | 25 ++++++ src/two_d/beam_elements/beam.rs | 68 +++++++++++++++++ src/two_d/beam_elements/cross_section.rs | 21 +++++ src/two_d/beam_elements/world.rs | 97 ++++++++++++++++++++++++ 7 files changed, 302 insertions(+) create mode 100644 src/bin/two_d.rs create mode 100644 src/lib.rs create mode 100644 src/two_d.rs create mode 100644 src/two_d/beam_elements.rs create mode 100644 src/two_d/beam_elements/beam.rs create mode 100644 src/two_d/beam_elements/cross_section.rs create mode 100644 src/two_d/beam_elements/world.rs diff --git a/src/bin/two_d.rs b/src/bin/two_d.rs new file mode 100644 index 0000000..03c18ac --- /dev/null +++ b/src/bin/two_d.rs @@ -0,0 +1,70 @@ +use fea::two_d::beam_elements::*; +use fea::two_d::*; + +fn main() { + let student_id_last_digit = 9 as f32; + let d = student_id_last_digit / 100.0; + let r = d/2.0; + let e = 210e9; + let p = (25.0 - student_id_last_digit)*1000.0; + let theta = 30.0_f32.to_radians(); + let l = 2.0/theta.cos(); + + let p1 = Point{ + id: 1, + pos: Vector(0.0, 0.0), + bc: BoundaryCondition::Fixed, + beams: vec![], + }; + + let p2 = Point { + id: 2, + pos: Vector(0.0, l * theta.sin()), + bc: BoundaryCondition::Force(Vector(0.0, 1000.0)), + beams: vec![], + }; + + let l_e3 = 2.0/theta.sin(); + println!("l_e3 {l_e3}"); + + let p3 = Point { + id: 3, + pos: Vector(0.0, (l * theta.sin()) + (l_e3 * theta.cos())), + bc: BoundaryCondition::Fixed, + beams: vec![], + }; + + let p4 = Point { + id: 4, + pos: Vector(l * theta.cos(), l * theta.sin()), + bc: BoundaryCondition::Force(Vector::from_mag_angle(p, 45_f32.to_radians())), + beams: vec![] + }; + + let mut w = World::from(vec![p1, p2, p3, p4]); + + w.link(1,4, NewBeam { + cross_section: CrossSection::Circular(r), + material: fea::Material { + yield_stress: 95e6, + youngs_modulus: e, + } + }).unwrap(); + w.link(2,4, NewBeam { + cross_section: CrossSection::Circular(r), + material: fea::Material { + yield_stress: 95e6, + youngs_modulus: e, + } + }).unwrap(); + w.link(3,4, NewBeam { + cross_section: CrossSection::Circular(r), + material: fea::Material { + yield_stress: 95e6, + youngs_modulus: e, + } + }).unwrap(); + + println!("{:?}", w); + println!("{:?}", w.stiffness_matrix().unwrap()); +} diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..3f2da4b --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,7 @@ +pub mod two_d; + +#[derive(Debug)] +pub struct Material { + pub youngs_modulus: f32, + pub yield_stress: f32, +} diff --git a/src/two_d.rs b/src/two_d.rs new file mode 100644 index 0000000..e8b803c --- /dev/null +++ b/src/two_d.rs @@ -0,0 +1,14 @@ +pub mod beam_elements; + +#[derive(Debug)] +pub struct Vector(pub f32, pub f32); + +impl Vector { + pub fn distance(&self, other: &Self) -> f32 { + ((self.0 - other.0).powi(2) + (self.1 - other.1).powi(2)).sqrt() + } + + pub fn from_mag_angle(magnitude: f32, angle: f32) -> Self { + Self (magnitude * angle.cos(), magnitude * angle.sin()) + } +} diff --git a/src/two_d/beam_elements.rs b/src/two_d/beam_elements.rs new file mode 100644 index 0000000..7532466 --- /dev/null +++ b/src/two_d/beam_elements.rs @@ -0,0 +1,25 @@ +mod beam; +mod cross_section; +mod world; + +pub use beam::*; +pub use cross_section::*; +pub use world::*; + +use super::Vector; +use crate::Material; + +#[derive(Debug)] +pub struct Point { + pub id: usize, + pub pos: Vector, + pub bc: BoundaryCondition, + pub beams: Vec, +} + +#[derive(Debug)] +pub enum BoundaryCondition { + Free, + Fixed, + Force(Vector), +} diff --git a/src/two_d/beam_elements/beam.rs b/src/two_d/beam_elements/beam.rs new file mode 100644 index 0000000..f88829c --- /dev/null +++ b/src/two_d/beam_elements/beam.rs @@ -0,0 +1,68 @@ +use super::*; + +#[derive(Debug)] +pub struct Beam{ + pub id: usize, + pub points: (usize, usize), + pub material: Material, + pub cross_section: CrossSection, +} + +pub struct NewBeam { + pub material: Material, + pub cross_section: CrossSection, +} + +impl Beam { + pub fn new(p1: usize, p2: usize, props: NewBeam) -> Beam { + let b = Beam { + id: 0, + points: (p1, p2), + material: props.material, + cross_section: props.cross_section, + }; + + return b; + } + + fn get_points<'a>(&self, world: &'a World) -> (&'a Point, &'a Point) { + let p1 = world.points.get(&self.points.0).unwrap(); + let p2 = world.points.get(&self.points.1).unwrap(); + + (p1, p2) + } + + pub fn stiffness(&self, world: &World) -> f32 { + let s = self.cross_section.area() * self.material.youngs_modulus / self.length(world); + + println!("{:?} {}", self, s); + + s + + } + + pub fn length(&self, world: &World) -> f32 { + let (p1, p2) = self.get_points(world); + let l = p1.pos.distance(&p2.pos); + + println!("{:?} {}", self, l); + + l + } + + pub fn angle(&self, world: &World) -> f32 { + let (p1, p2) = self.get_points(world); + let dx = p1.pos.0 - p2.pos.0; + let dy = p1.pos.1 - p2.pos.1; + + (dy/dx).atan() + } + + pub fn other_point(&self, p: usize) -> usize { + if p == self.points.0{ + return self.points.1; + } else { + return self.points.0; + } + } +} diff --git a/src/two_d/beam_elements/cross_section.rs b/src/two_d/beam_elements/cross_section.rs new file mode 100644 index 0000000..ad77ad1 --- /dev/null +++ b/src/two_d/beam_elements/cross_section.rs @@ -0,0 +1,21 @@ +use std::f32::consts::PI; + +#[derive(Debug)] +pub enum CrossSection { + Area(f32), + Circular(f32), + Rectangular(f32, f32), +} + +impl CrossSection { + pub fn area(&self) -> f32 { + let a = match self { + Self::Area(a) => *a, + Self::Circular(r) => PI * r.powi(2), + Self::Rectangular(l1, l2) => l1 * l2, + }; + println!("{:?} {}", self, a); + + a + } +} diff --git a/src/two_d/beam_elements/world.rs b/src/two_d/beam_elements/world.rs new file mode 100644 index 0000000..7549eda --- /dev/null +++ b/src/two_d/beam_elements/world.rs @@ -0,0 +1,97 @@ +use super::*; +use ndarray::prelude::*; +use std::collections::HashMap; + +#[derive(Debug)] +pub struct World { + pub points: HashMap, + pub beams: HashMap, +} + +impl World { + pub fn stiffness_matrix(&self) -> Result, String> { + let dof = self.points.len() * 2; + let mut a: Vec = vec![0.0; dof.pow(2)]; + + // set up forces exerted by displacements of other beams + for (_point_id, point) in &self.points { + for beam_id in &point.beams { + let beam = self.beams.get(beam_id).unwrap(); + let stiffness = beam.stiffness(self); + let other_point = self.points.get(&beam.other_point(point.id)).unwrap(); + let ax = 2*(point.id-1); + let ay = ax + 1; + let bx = 2*(other_point.id-1); + let by = bx + 1; + + let angle = beam.angle(self); + println!("point: {} to point {} beam: {} beam_angle: {}", _point_id, other_point.id, beam.id, beam.angle(self).to_degrees()); + let c2 = angle.cos().powi(2); + let s2 = angle.sin().powi(2); + let cs = angle.cos() * angle.sin(); + + // F_AX += K_AB * c2 * (U_AX) + a[ax*dof + ax] += stiffness * c2; + // F_AX += K_AB * cs * U_AY + a[ax*dof + ay] += stiffness * cs; + // F_AY += K_AB * cs * U_AX + a[ay*dof + ax] += stiffness * cs; + // F_AY += K_AB * s^2 * U_AY + a[ay*dof + ay] += stiffness * s2; + + // F_AX += K_AB -c^2 U_BX + a[ax*dof + bx] += stiffness * -c2; + // F_AX += K_AB -cs U_BY + a[ax*dof + by] += stiffness * -cs; + // F_AY += K_AB -cs U_BX + a[ay*dof + bx] += stiffness * -cs; + // F_AY += K_AB -s^2 U_BY + a[ay*dof + by] += stiffness * -s2; + } + } + + match Array2::from_shape_vec((dof, dof), a) { + Ok(arr) => Ok(arr), + Err(e) => Err(e.to_string()), + } + } + + pub fn link(&mut self, id1: usize, id2: usize, new_beam: NewBeam) -> Result<(), &str>{ + let b_id = self.beams.len(); + let p1 = match self.points.get_mut(&id1) { + Some(point) => point, + None => return Err("Point with id1 not found"), + }; + p1.beams.push(b_id); + let p1_id = p1.id; + + let p2 = match self.points.get_mut(&id2) { + Some(point) => point, + None => return Err("Point with id1 not found"), + }; + let p2_id = p2.id; + p2.beams.push(b_id); + + let mut b = Beam::new(p1_id, p2_id, new_beam); + b.id = b_id; + self.beams.insert(b_id, b); + + Ok(()) + } +} + +impl From> for World { + /// Create a world from a list of points. + /// Expects each point to have a unique ID. + /// If this is not the case, then some points will be overriden. + fn from(points: Vec) -> World { + let mut points = points; + let mut w = World { points: HashMap::new(), beams: HashMap::new() }; + + while let Some(point) = points.pop() { + w.points.insert(point.id, point); + } + + w + } +}