Time-dependent dynamic problem with nonlinear stiffness using ode45
44 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Andres Bellei
el 8 de Jul. de 2024 a las 19:40
I am solving a time-dependent finite element problem with nonlinearities in the stiffness matrix which takes the following form:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1730356/image.png)
where the dot symbolizes derivative with respect to time, U is the displacement vector, M is the mass matrix, C is the damping matrix, K(U) is the stiffness matrix which depends on U, and F is an external force. I would like to use ode45 to solve the dynamic problem. The general flow for solving this problem would be the following:
- Calculate the stiffness matrix K using the latest displacement U (Initial condition if it is the first time step).
- Formulate the state-space equations using K matrix, M matrix, C matrix, and F vector.
- Integrate the state-space equations with ode45 to caclulate the new displacement U.
- Recalculate the K matrix with the new displacement U. Check if the residual error meets the preselected tolerance.
- If the error is met, update velocities and accelerations, and go to the next time step. If the error is not met, then go back to step 1 using this latest K matrix.
I am familiar with solving this problem without the nonlinearity in K. However, I am unsure how to update the K matrix at every time step within the ode45 framework.
Thank you for the guidance!
0 comentarios
Respuestas (1)
Torsten
el 8 de Jul. de 2024 a las 19:50
Movida: Torsten
el 8 de Jul. de 2024 a las 19:50
You solve for U and Udot - thus you get values for U and Udot from ode45. Simply compute the matrix K(U) with these values and proceed as if K was constant.
4 comentarios
Torsten
el 9 de Jul. de 2024 a las 17:04
Editada: Torsten
el 9 de Jul. de 2024 a las 17:28
Example:
M = [1 0;0 4];
invM = inv(M);
C = [4 -2;6 3];
K = @(u)[u(1) sin(u(2));exp(-u(2)) 1/u(1)];
F = [0 ;0];
tspan = [0 1];
y0 = [1; 2 ;3; 4]; % vector of initial values for (u(1),u(2),udot(1),udot(2))
[T,Y] = ode45(@(t,y)fun(t,y,invM,C,K,F),tspan,y0);
plot(T,[Y(:,1),Y(:,2)])
grid on
function dy = fun(t,y,invM,C,K,F)
n = numel(y)/2;
u = y(1:n);
udot = y(n+1:2*n);
dy = zeros(2*n,1);
dy(1:n) = udot;
dy(n+1:2*n) = invM*(F-C*udot-K(u)*u);
end
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!