Structural-Mechanical
/
FEM-Course-Matlab
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
FEM-Course-Matlab/16.几何非线性有限元matlab编程/几何非线性有限元-Williams frame/solve.m

213 lines
9.3 KiB
Mathematica
Raw Permalink Normal View History

init project
2024-01-28 16:46:36 +00:00
%% Function to solve nonlinear system of equations
%
% INPUT ARGUMENTS:
% * Model -> Structure with model parameters:
% - nnp: Number of nodes
% - nel: Number of elements
% - neq: Number of equation
% - neqf: Number of free DOF equations
% - neqc: Number of constrained DOF equations
% * Anl -> Structure with analysis parameters:
% - form: Nonlinear beam formulation
% - alg: Solution algorithm
% - incr_type: Increment size update strategy in predictor step
% - iter_type: Iteration strategy
% - geom_mtx: Geometric stiffness matrix
% - init_incr: Initial load ratio increment
% - max_ratio: Maximum load ratio
% - max_steps: Maximum number of steps
% - max_iter: Maximum number of iterations
% - trg_iter: Target number of iterations
% - tol: Tolerance for iteration convergence
% * Node -> Vector of structures with nodes information:
% - x: X-axis coordinate
% - y: Y-axis coordinate
% - px: Nodal load in global X-axis direction
% - py: Nodal load in global Y-axis direction
% - mz: Concentrated moment about global Z-axis
% - dof: Vector with global DOF numbers
% * Elem -> Vector of structures with elements information:
% - n1: Handle to initial node
% - n2: Handle to final node
% - E: Elasticity modulus
% - A: Cross-section area
% - I: Cross-section moment of inertia
% - L: Length
% - L_1: Length from beggining of step
% - fi: Vector of internal forces in local system [fx1,fy1,mz1,fx2,fy2,mz2]
% - fi_1: Vector of internal forces in local system from beggining of step
% - fn: Vector of internal forces in natural system [N2,M1,M2]
% - fn_1: Vector of internal forces in natural system from beggining of step
% - angle: Angle with horizontal axis
% - angle_1: Angle with horizontal axis from beggining of step
% - gle: Gather vector with global DOF numbers
% - rot: Rotation matrix from global to local coordinate system
% - ke: Elastic stiffness matrix
% - kt: Tangent stiffness matrix
% * Pref -> Reference load vector (total applied nodal load)
%
% OUTPUT:
% * Result -> Vector of structures with results to be plotted:
% - U_step: Matrix of nodal displacements (each column is the displacements vector in each step of analysis)
% - U_iter: Matrix of nodal displacements (each column is the displacements vector in each iteration of analysis)
% - lr_step: Vector of load ratio in all steps
% - lr_iter: Vector of load ratio in all iterations
%
% NOTATION:
% * lr -> load ratio
% * U -> Vector of nodal displacements
% * D_ -> Indicate accumulated step increment of a variable
% * d_ -> Indicate iterative increment of a variable
%
function Result = solve(Model,Anl,Elem,Pref)
% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD>λ<EFBFBD>ƺͺ<EFBFBD><EFBFBD>ر<EFBFBD>
U = zeros(Model.neq,1);
lr = 0;
% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD>result<EFBFBD>ṹ<EFBFBD>壬<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʼƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Result = struct('U_step',[],'U_iter',[],'lr_step',[],'lr_iter',[]);
Result.U_step(:,1) = U;
Result.U_iter(:,1) = U;
Result.lr_step(1) = lr;
Result.lr_iter(1) = lr;
%====================================================================================================================================
% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
step = 0;
while (step < Anl.max_steps)
step = step + 1;
% <EFBFBD><EFBFBD><EFBFBD>²ο<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>µ<EFBFBD>UL<EFBFBD><EFBFBD><EFBFBD>̣<EFBFBD>L_1<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>U<EFBFBD><EFBFBD><EFBFBD>и<EFBFBD><EFBFBD>£<EFBFBD>angle_1<EFBFBD><EFBFBD><EFBFBD>ŵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>²<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>angle<EFBFBD><EFBFBD>Ϊ<EFBFBD>ò<EFBFBD>angle_1;<EFBFBD><EFBFBD><EFBFBD><EFBFBD>fi_1ͬangle<EFBFBD><EFBFBD>
% <EFBFBD>˴<EFBFBD><EFBFBD><EFBFBD><EFBFBD>µ<EFBFBD>L_1,angle_1<EFBFBD><EFBFBD>phi_1<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>õ<EFBFBD>intForcesUL
for i = 1:Model.nel
Elem(i).L_1 = elemLength(Elem(i),U); %L_1 Length from beggining of step
Elem(i).angle_1 = Elem(i).angle; %angle_1 Angle with horizontal axis from beggining of step
Elem(i).fi_1 = Elem(i).fi; %fi_1 Vector of internal forces in local system from beggining of step
end
% <EFBFBD><EFBFBD><EFBFBD>߸նȾ<EFBFBD><EFBFBD>󣨸<EFBFBD><EFBFBD><EFBFBD>Ui-1<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ki0)
[Kt,Elem] = tangStiffMtxUL(Model,Anl,Elem,U) ; %%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>i-1<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>½ڵ<EFBFBD>λ<EFBFBD>ú<EFBFBD><EFBFBD>õ<EFBFBD><EFBFBD>ĵ<EFBFBD>Ԫ<EFBFBD>նȾ<EFBFBD><EFBFBD>󣨵<EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD>ȣ<EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD>Ƕȣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ص<EFBFBD>Elemֻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ke<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% Ԥ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD>Ƶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>delta_Ui1
d_Ut = solveLinearSystem(Model,Kt,Pref);%delta_Ui1<EFBFBD><EFBFBD>Ԥ<EFBFBD><EFBFBD><EFBFBD>׶ε<EFBFBD>λ<EFBFBD>Ƶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ǽ<EFBFBD><EFBFBD><EFBFBD>Ԥ<EFBFBD><EFBFBD><EFBFBD>׶εĺ<EFBFBD><EFBFBD>ر<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ƿ<EFBFBD>Ϊ<EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD>ĺ<EFBFBD><EFBFBD>ر<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD>ͬ<EFBFBD><EFBFBD>
if (step == 1)
s = 1; %<EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԥ<EFBFBD><EFBFBD><EFBFBD>׶ε<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD><EFBFBD><EFBFBD>Ĭ<EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD>
d_lr = Anl.init_incr; %delta_lamda_11 <EFBFBD><EFBFBD>Ϊ<EFBFBD>涨Ԥ<EFBFBD><EFBFBD><EFBFBD>׶<EFBFBD>d_lr
n2 = d_Ut'*d_Ut;%<EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD>ֵ,<EFBFBD><EFBFBD><EFBFBD>ں<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>GSP
else
% Generalized Stiffness Parameter
GSP = n2 / (d_Ut'*d_Ut_1);%<EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD>41<EFBFBD><EFBFBD> d_Ut_1 <EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> d_Ut <EFBFBD><EFBFBD>ǰ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
if (GSP < 0)
s = -s; %<EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD>42<EFBFBD><EFBFBD>
end
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>J
% J = sqrt((Anl.trg_iter + 1) / (iter + 1));%Ŀ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>/<EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD>40<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>뽫<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> 1 <EFBFBD>Ա<EFBFBD><EFBFBD>ⱻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ΪԤ<EFBFBD><EFBFBD><EFBFBD>Ľ<EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ε<EFBFBD><EFBFBD><EFBFBD>
J = 1;
%% Ԥ<EFBFBD><EFBFBD><EFBFBD>׶κ<EFBFBD><EFBFBD>ر<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ALG_ALCM_CYL <EFBFBD>̶<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>-Բ<EFBFBD><EFBFBD>
% Extract free DOF components
D_U_temp = D_U(1:Model.neqf);
d_Ut_temp = d_Ut(1:Model.neqf);
d_lr = J * sqrt((D_U_temp'*D_U_temp) /(d_Ut_temp'*d_Ut_temp));%<EFBFBD><EFBFBD>ʽ38<EFBFBD><EFBFBD>Բ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% d_lr = J * abs(d_lr_1);%% ALG_LCM (Load Increment)
% d_lr = J * sqrt(abs((D_lr*Pref'*D_U) /(Pref'*d_Ut))); %ALG_WCM (Work Increment)
% Apply correct sign
d_lr = s * d_lr;
end
% <EFBFBD><EFBFBD><EFBFBD>ر<EFBFBD>ҪС<EFBFBD>ڹ涨ֵ<EFBFBD><EFBFBD>==1<EFBFBD><EFBFBD>
if ((Anl.max_ratio > 0.0 && lr + d_lr > Anl.max_ratio) ||...
(Anl.max_ratio < 0.0 && lr + d_lr < Anl.max_ratio))
d_lr = Anl.max_ratio - lr;%<EFBFBD><EFBFBD>֤ǡ<EFBFBD>ôﵽ<EFBFBD><EFBFBD><EFBFBD>յĺ<EFBFBD><EFBFBD>ر<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>1ʵ<EFBFBD>ֺ<EFBFBD><EFBFBD>ص<EFBFBD>ȫ<EFBFBD><EFBFBD>ʩ<EFBFBD><EFBFBD>
end
% Ԥ<EFBFBD><EFBFBD><EFBFBD>׶<EFBFBD>λ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
d_U = d_lr * d_Ut;%<EFBFBD><EFBFBD><EFBFBD>ݹ<EFBFBD>ʽ<EFBFBD><EFBFBD>36<EFBFBD><EFBFBD>Ԥ<EFBFBD><EFBFBD><EFBFBD>׶μٶ<EFBFBD><EFBFBD>в<EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ˣ<EFBFBD>ֻ<EFBFBD><EFBFBD>ǰһ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% % <EFBFBD>洢Ԥ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
d_lr_1 = d_lr;
d_Ut_1 = d_Ut;%<EFBFBD>ж<EFBFBD>GSP<EFBFBD><EFBFBD><EFBFBD>ţ<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD>41<EFBFBD><EFBFBD>
d_U_1 = d_U;
% ǰ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>غɱȺ<EFBFBD>λ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ŵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϸ<EFBFBD><EFBFBD>µ<EFBFBD><EFBFBD>ӣ<EFBFBD><EFBFBD><EFBFBD>дD<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD><EFBFBD>λ<EFBFBD>ƺͺ<EFBFBD><EFBFBD>ر<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Сдd<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD>ƺͺ<EFBFBD><EFBFBD>ر<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
D_lr = d_lr;
D_U = d_U;
% <EFBFBD>ܺ<EFBFBD><EFBFBD>رȺ<EFBFBD>λ<EFBFBD><EFBFBD>
lr = lr + d_lr;
U = U + d_U;
%-------------------------------------------------------------------------------------------------------------------------------
% Start iterative process
iter = 0;
conv = 0;
while (conv == 0 && iter < Anl.max_iter)
% External and internal forces
P = lr * Pref;
[F,Elem] = intForcesUL(Model,Elem,U,D_U,1);%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD>ε<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>õ<EFBFBD><EFBFBD><EFBFBD>Ke,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>D_U<EFBFBD><EFBFBD>Ele.angle<EFBFBD><EFBFBD>Ele.fi<EFBFBD><EFBFBD><EFBFBD>и<EFBFBD><EFBFBD><EFBFBD>
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ǻ<EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD>Ƽ<EFBFBD><EFBFBD>㣬<EFBFBD><EFBFBD>λ<EFBFBD>ƴ<EFBFBD><EFBFBD>ڷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ԣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>в<EFBFBD>
% Unbalanced forces
R = P - F;
% Store iteration results
Result.U_iter(:,size(Result.U_iter,2)+1) = U;
Result.lr_iter(size(Result.lr_iter,2)+1) = lr;
% Check for convergence
conv = (norm(R(1:Model.neqf))/norm(Pref(1:Model.neqf)) < Anl.tol);
if (conv == 1)
break;
end
% Start/keep corrective cycle of iterations
iter = iter + 1;
[Kt,Elem] = tangStiffMtxUL(Model,Anl,Elem,U); %ÿ<EFBFBD>ε<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>£<EFBFBD><EFBFBD><EFBFBD>׼ţ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɭ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Elem.ke<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% Tangent and residual increments of displacements
d_Ut = solveLinearSystem(Model,Kt,Pref);
d_Ur = solveLinearSystem(Model,Kt,R);
% <EFBFBD><EFBFBD><EFBFBD>ر<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD>45<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>㶨<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>-Բ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
a = d_Ut'*d_Ut;
b = d_Ut'*(d_Ur + D_U);
c = d_Ur'*(d_Ur + 2*D_U);%D_UΪ<EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD>d_UrΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>º<EFBFBD><EFBFBD><EFBFBD>ֵ
s_iter = sign(D_U'*d_Ut);
d_lr = -b/a + s_iter*sqrt((b/a)^2 - c/a);
%Ϊ<EFBFBD>λ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ָ<EFBFBD><EFBFBD><EFBFBD>
if (~isreal(d_lr))
conv = -1;
break;
end
% <EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD>36<EFBFBD><EFBFBD>
d_U = d_lr * d_Ut + d_Ur;
% Increments of load ratio and displacements for current step
D_lr = D_lr + d_lr;
D_U = D_U + d_U;
% Total values of load ratio and displacements
lr = lr + d_lr;
U = U + d_U;
end
%--------------------------------------------------------------------------------------------------------------------------------
% Check for convergence fail reached or complex value of increment
if (conv == 0 && Anl.method == METH_INCR_ITER)
fprintf(1,'Convergence not reached!\nStep = %d\nIter = %d\n',step,iter);
return;
elseif (conv == -1)
fprintf(1,'Unable to compute load increment!\nStep = %d\nIter = %d\n',step,iter);
return;
end
% Store equilibrium configuration
Result.U_step(:,step+1) = U;
Result.lr_step(step+1) = lr;
fprintf(1,'Step:%d | ratio:%.2f | Iter:%d\n',step,lr,iter);
% Check if maximum load ratio was reached
if (lr >= 0.999*Anl.max_ratio)
break;
end
end
end