Skip to content

Von Karaman Beam Solver

Class to solve the von Karman beam bending problem using numerical methods.

Source code in lib/BeamModel.py 
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
class VonKarmanBeamSolver:
    """
    Class to solve the von Karman beam bending problem using numerical methods.
    """

    def __init__(self, P, L, E, b, h, point_offset, time_steps=100, kappa=5/6, poisson_ratio=0.33, yield_strength=250, fatigue_strength=160):
        """
        Initializes the Von Karman Beam Solver.

        Args:
            P (float): Maximum load in Newtons.
            L (float): Length of the beam in mm.
            E (float): Young's Modulus in MPa (N/mm^2).
            b (float): Width of the beam in mm.
            h (float): Thickness of the beam in mm.
            point_offset (float): Offset from the center of the beam in mm.
            time_steps (int, optional): Number of time steps for force increase. Defaults to 100.
            kappa (float, optional): Shear correction factor. Defaults to 5/6.
            poisson_ratio (float, optional): Poisson's ratio of the material. Defaults to 0.33.
            yield_strength (float, optional): Yield strength of the material in MPa. Defaults to 250 MPa.
            fatigue_strength (float, optional): Fatigue strength of the material at 1e7 cycles in MPa. Defaults to 160 MPa.
        """
        self.P_max = P
        self.L = L
        self.E = E
        self.G = E / (2 * (1 + poisson_ratio))  # Shear modulus
        self.b = b
        self.h = h
        self.point_offset = point_offset
        self.kappa = kappa
        self.time_steps = time_steps
        self.yield_strength = yield_strength
        self.fatigue_strength = fatigue_strength

        # Cross-sectional properties
        self.I = (b * h ** 3) / 12  # Second moment of area for rectangular cross-section
        self.A = b * h  # Cross-sectional area

        # Initialize safety factor lists
        self.safety_factors_yield = []  # List to store safety factors for yield strength
        self.safety_factors_fatigue = []  # List to store safety factors for fatigue strength

        # Variables to store forces at which safety factors fall below 2
        self.force_below_yield_sf = None
        self.force_below_fatigue_sf = None

    @classmethod
    def from_yaml(cls, yaml_file):
        """
        Load parameters from a YAML file and initialize the class.

        Args:
            yaml_file (str): Path to the YAML file containing beam configuration parameters.

        Returns:
            VonKarmanBeamSolver: Instance of the VonKarmanBeamSolver class.
        """
        with open(yaml_file, 'r') as file:
            data = yaml.safe_load(file)
        return cls(
            P=data['P'],
            L=data['L'],
            E=data['E'],
            b=data['b'],
            h=data['h'],
            point_offset=data['point_offset'],
            time_steps=data.get('time_steps', 100),
            kappa=data.get('kappa', 5 / 6),
            poisson_ratio=data.get('poisson_ratio', 0.33)
        )

    def von_karman_beam(self, t, y, P):
        """
        Differential equation for Von Karman beam bending.

        Args:
            t (float): Position along the beam.
            y (list of float): State vector [deflection, slope, horizontal displacement].
            P (float): Applied load at the current time step.

        Returns:
            list of float: List of derivatives [dw/dx, d²w/dx², du/dx].
        """
        w, dw_dx, u = y  # Deflection, slope, and horizontal displacement
        d2w_dx2 = (P * (self.L - t)) / (self.E * self.I)  # Curvature
        du_dx = 0.5 * (dw_dx ** 2)  # Horizontal stretch
        return [dw_dx, d2w_dx2, du_dx]

    def solve(self, save_path=None):
        """
        Solve the beam equations using numerical integration.

        Args:
            save_path (str, optional): Path to save the results as a CSV file. Defaults to None.

        Returns:
            None
        """
        y0 = [0, 0, 0]  # Initial conditions: w(0) = 0, slope = 0, u(0) = 0
        t_eval = np.linspace(0, self.L, 100)  # Evaluation points along the beam

        # Store solutions for each time step where the force P increases linearly
        self.solutions = []
        self.applied_force = []

        # Solve for t=0 with no force
        sol_initial = solve_ivp(self.von_karman_beam, [0, self.L], y0, t_eval=t_eval, args=(0,))
        self.solutions.append(sol_initial)

        # Solve for subsequent timesteps where force increases linearly
        for t_step in range(1, self.time_steps + 1):
            P_current = self.P_max * (t_step / self.time_steps)  # Linearly increase force
            self.applied_force.append(P_current)
            sol = solve_ivp(self.von_karman_beam, [0, self.L], y0, t_eval=t_eval, args=(P_current,))
            self.solutions.append(sol)

            # Calculate yield and fatigue safety factors
            yield_sf = self.calculate_safety_factor(P_current)
            fatigue_sf = self.calculate_fatigue_safety_factor(P_current)
            self.safety_factors_yield.append(yield_sf)
            self.safety_factors_fatigue.append(fatigue_sf)

            # Record the force when safety factors fall below 2
            if yield_sf < 2 and self.force_below_yield_sf is None:
                self.force_below_yield_sf = P_current

            if fatigue_sf < 2 and self.force_below_fatigue_sf is None:
                self.force_below_fatigue_sf = P_current

        # Extract the final solution for plotting and saving
        final_solution = self.solutions[-1]
        self.x = final_solution.t  # Position along the beam
        self.w = final_solution.y[0]  # Vertical deflection
        self.u = final_solution.y[2]  # Horizontal displacement
        self.slope = final_solution.y[1]  # Slope (dw/dx)
        self.angle = np.degrees(self.slope)  # Slope in degrees

        # Calculate total displacements
        self.u_total = self.u + self.point_offset * np.sin(self.slope)
        self.w_total = self.w + self.point_offset * (1 - np.cos(self.slope))

        if save_path:
            # Save data to a CSV file for the final time step
            data = {
                'Position (x)': self.x,
                'Vertical Deflection (w)': self.w_total,
                'Horizontal Displacement (u)': self.u_total,
                'Angle of Slope (degrees)': self.angle,
                'Safety Factor (Yield)': self.safety_factors_yield[-1],
                'Safety Factor (Fatigue)': self.safety_factors_fatigue[-1]
            }
            df = pd.DataFrame(data)
            df.to_csv(save_path, index=False)
            print(f"Data saved to {save_path}")

    def calculate_safety_factor(self, P):
        """
        Calculate the safety factor for yielding in bending.

        Args:
            P (float): Current applied load in Newtons.

        Returns:
            float: Safety factor for bending.
        """
        M_max = P * self.L / 4  # Maximum bending moment
        c = self.h / 2  # Distance from neutral axis to outermost fiber
        sigma_max = (M_max * c) / self.I  # Maximum stress in bending
        safety_factor = self.yield_strength / sigma_max
        return safety_factor

    def calculate_fatigue_safety_factor(self, P):
        """
        Calculate the safety factor for fatigue at 1e7 cycles.

        Args:
            P (float): Current applied load in Newtons.

        Returns:
            float: Safety factor for fatigue.
        """
        M_max = P * self.L / 4  # Maximum bending moment
        c = self.h / 2  # Distance from neutral axis to outermost fiber
        sigma_max = (M_max * c) / self.I  # Maximum bending stress
        safety_factor_fatigue = self.fatigue_strength / sigma_max
        return safety_factor_fatigue

    def plot(self, save_path=None):
        """
        Generate an interactive plot using Plotly to visualize beam deflection and slope.

        Args:
            save_path (str, optional): Path to save the plot as an HTML file. Defaults to None.

        Returns:
            None
        """
        fig = go.Figure()

        # Horizontal movement trace
        fig.add_trace(go.Scatter(x=self.x, y=self.u_total, mode='lines', name='Horizontal Movement (u)',
                                 line=dict(color='blue')))

        # Vertical movement trace
        fig.add_trace(go.Scatter(x=self.x, y=self.w_total, mode='lines', name='Vertical Movement (w)',
                                 line=dict(color='red', dash='dash')))

        # Angle trace
        fig.add_trace(go.Scatter(x=self.x, y=self.angle, mode='lines', name='Angle of Slope (degrees)',
                                 line=dict(color='green', dash='dot')))

        # Layout settings
        fig.update_layout(title='Movement of the Point at the Bending Edge',
                          xaxis_title='Position along the beam (mm)',
                          yaxis_title='Displacement (mm) / Angle (degrees)',
                          legend=dict(x=0.05, y=0.95),
                          template='plotly_white')

        if save_path:
            # Save plot as an HTML file
            fig.write_html(save_path)
            print(f"Plot saved to {save_path}")
        # Show the interactive plot
        fig.show()

from_yaml(yaml_file) classmethod

Load parameters from a YAML file and initialize the class.

Parameters:

Name Type Description Default
yaml_file str

Path to the YAML file containing beam configuration parameters.

 required

Returns:

Name Type Description
VonKarmanBeamSolver

Instance of the VonKarmanBeamSolver class.
Source code in lib/BeamModel.py 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
@classmethod
def from_yaml(cls, yaml_file):
    """
    Load parameters from a YAML file and initialize the class.

    Args:
        yaml_file (str): Path to the YAML file containing beam configuration parameters.

    Returns:
        VonKarmanBeamSolver: Instance of the VonKarmanBeamSolver class.
    """
    with open(yaml_file, 'r') as file:
        data = yaml.safe_load(file)
    return cls(
        P=data['P'],
        L=data['L'],
        E=data['E'],
        b=data['b'],
        h=data['h'],
        point_offset=data['point_offset'],
        time_steps=data.get('time_steps', 100),
        kappa=data.get('kappa', 5 / 6),
        poisson_ratio=data.get('poisson_ratio', 0.33)
    )

von_karman_beam(t, y, P)

Differential equation for Von Karman beam bending.

Parameters:

Name Type Description Default
t float

Position along the beam.

 required
y list of float

State vector [deflection, slope, horizontal displacement].

 required
P float

Applied load at the current time step.

 required

Returns:

Type Description

list of float: List of derivatives [dw/dx, d²w/dx², du/dx].
Source code in lib/BeamModel.py 
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
def von_karman_beam(self, t, y, P):
    """
    Differential equation for Von Karman beam bending.

    Args:
        t (float): Position along the beam.
        y (list of float): State vector [deflection, slope, horizontal displacement].
        P (float): Applied load at the current time step.

    Returns:
        list of float: List of derivatives [dw/dx, d²w/dx², du/dx].
    """
    w, dw_dx, u = y  # Deflection, slope, and horizontal displacement
    d2w_dx2 = (P * (self.L - t)) / (self.E * self.I)  # Curvature
    du_dx = 0.5 * (dw_dx ** 2)  # Horizontal stretch
    return [dw_dx, d2w_dx2, du_dx]

solve(save_path=None)

Solve the beam equations using numerical integration.

Parameters:

Name Type Description Default
save_path str

Path to save the results as a CSV file. Defaults to None.

 None

Returns:

Type Description

None
Source code in lib/BeamModel.py 
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
def solve(self, save_path=None):
    """
    Solve the beam equations using numerical integration.

    Args:
        save_path (str, optional): Path to save the results as a CSV file. Defaults to None.

    Returns:
        None
    """
    y0 = [0, 0, 0]  # Initial conditions: w(0) = 0, slope = 0, u(0) = 0
    t_eval = np.linspace(0, self.L, 100)  # Evaluation points along the beam

    # Store solutions for each time step where the force P increases linearly
    self.solutions = []
    self.applied_force = []

    # Solve for t=0 with no force
    sol_initial = solve_ivp(self.von_karman_beam, [0, self.L], y0, t_eval=t_eval, args=(0,))
    self.solutions.append(sol_initial)

    # Solve for subsequent timesteps where force increases linearly
    for t_step in range(1, self.time_steps + 1):
        P_current = self.P_max * (t_step / self.time_steps)  # Linearly increase force
        self.applied_force.append(P_current)
        sol = solve_ivp(self.von_karman_beam, [0, self.L], y0, t_eval=t_eval, args=(P_current,))
        self.solutions.append(sol)

        # Calculate yield and fatigue safety factors
        yield_sf = self.calculate_safety_factor(P_current)
        fatigue_sf = self.calculate_fatigue_safety_factor(P_current)
        self.safety_factors_yield.append(yield_sf)
        self.safety_factors_fatigue.append(fatigue_sf)

        # Record the force when safety factors fall below 2
        if yield_sf < 2 and self.force_below_yield_sf is None:
            self.force_below_yield_sf = P_current

        if fatigue_sf < 2 and self.force_below_fatigue_sf is None:
            self.force_below_fatigue_sf = P_current

    # Extract the final solution for plotting and saving
    final_solution = self.solutions[-1]
    self.x = final_solution.t  # Position along the beam
    self.w = final_solution.y[0]  # Vertical deflection
    self.u = final_solution.y[2]  # Horizontal displacement
    self.slope = final_solution.y[1]  # Slope (dw/dx)
    self.angle = np.degrees(self.slope)  # Slope in degrees

    # Calculate total displacements
    self.u_total = self.u + self.point_offset * np.sin(self.slope)
    self.w_total = self.w + self.point_offset * (1 - np.cos(self.slope))

    if save_path:
        # Save data to a CSV file for the final time step
        data = {
            'Position (x)': self.x,
            'Vertical Deflection (w)': self.w_total,
            'Horizontal Displacement (u)': self.u_total,
            'Angle of Slope (degrees)': self.angle,
            'Safety Factor (Yield)': self.safety_factors_yield[-1],
            'Safety Factor (Fatigue)': self.safety_factors_fatigue[-1]
        }
        df = pd.DataFrame(data)
        df.to_csv(save_path, index=False)
        print(f"Data saved to {save_path}")

calculate_safety_factor(P)

Calculate the safety factor for yielding in bending.

Parameters:

Name Type Description Default
P float

Current applied load in Newtons.

 required

Returns:

Name Type Description
float

Safety factor for bending.
Source code in lib/BeamModel.py 
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
def calculate_safety_factor(self, P):
    """
    Calculate the safety factor for yielding in bending.

    Args:
        P (float): Current applied load in Newtons.

    Returns:
        float: Safety factor for bending.
    """
    M_max = P * self.L / 4  # Maximum bending moment
    c = self.h / 2  # Distance from neutral axis to outermost fiber
    sigma_max = (M_max * c) / self.I  # Maximum stress in bending
    safety_factor = self.yield_strength / sigma_max
    return safety_factor

calculate_fatigue_safety_factor(P)

Calculate the safety factor for fatigue at 1e7 cycles.

Parameters:

Name Type Description Default
P float

Current applied load in Newtons.

 required

Returns:

Name Type Description
float

Safety factor for fatigue.
Source code in lib/BeamModel.py 
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
def calculate_fatigue_safety_factor(self, P):
    """
    Calculate the safety factor for fatigue at 1e7 cycles.

    Args:
        P (float): Current applied load in Newtons.

    Returns:
        float: Safety factor for fatigue.
    """
    M_max = P * self.L / 4  # Maximum bending moment
    c = self.h / 2  # Distance from neutral axis to outermost fiber
    sigma_max = (M_max * c) / self.I  # Maximum bending stress
    safety_factor_fatigue = self.fatigue_strength / sigma_max
    return safety_factor_fatigue

plot(save_path=None)

Generate an interactive plot using Plotly to visualize beam deflection and slope.

Parameters:

Name Type Description Default
save_path str

Path to save the plot as an HTML file. Defaults to None.

 None

Returns:

Type Description

None
Source code in lib/BeamModel.py 
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
def plot(self, save_path=None):
    """
    Generate an interactive plot using Plotly to visualize beam deflection and slope.

    Args:
        save_path (str, optional): Path to save the plot as an HTML file. Defaults to None.

    Returns:
        None
    """
    fig = go.Figure()

    # Horizontal movement trace
    fig.add_trace(go.Scatter(x=self.x, y=self.u_total, mode='lines', name='Horizontal Movement (u)',
                             line=dict(color='blue')))

    # Vertical movement trace
    fig.add_trace(go.Scatter(x=self.x, y=self.w_total, mode='lines', name='Vertical Movement (w)',
                             line=dict(color='red', dash='dash')))

    # Angle trace
    fig.add_trace(go.Scatter(x=self.x, y=self.angle, mode='lines', name='Angle of Slope (degrees)',
                             line=dict(color='green', dash='dot')))

    # Layout settings
    fig.update_layout(title='Movement of the Point at the Bending Edge',
                      xaxis_title='Position along the beam (mm)',
                      yaxis_title='Displacement (mm) / Angle (degrees)',
                      legend=dict(x=0.05, y=0.95),
                      template='plotly_white')

    if save_path:
        # Save plot as an HTML file
        fig.write_html(save_path)
        print(f"Plot saved to {save_path}")
    # Show the interactive plot
    fig.show()