213 lines
9.3 KiB
Matlab
213 lines
9.3 KiB
Matlab
%% 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)
|
||
|
||
|
||
% <20><>ʼ<EFBFBD><CABC><EFBFBD>ڵ<EFBFBD>λ<EFBFBD>ƺͺ<C6BA><CDBA>ر<EFBFBD>
|
||
U = zeros(Model.neq,1);
|
||
lr = 0;
|
||
% <20><>ʼ<EFBFBD><CABC>result<6C>ṹ<EFBFBD>壬<EFBFBD><E5A3AC><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʼƽ<CABC><C6BD><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;
|
||
|
||
%====================================================================================================================================
|
||
% <20><>ʼ<EFBFBD><CABC><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
||
step = 0;
|
||
while (step < Anl.max_steps)
|
||
step = step + 1;
|
||
% <20><><EFBFBD>²ο<C2B2><CEBF><EFBFBD><EFBFBD><EFBFBD><EFBFBD>µ<EFBFBD>UL<55><4C><EFBFBD>̣<EFBFBD>L_1<5F><31><EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1>λ<EFBFBD><CEBB>U<EFBFBD><55><EFBFBD>и<EFBFBD><D0B8>£<EFBFBD>angle_1<5F><31><EFBFBD>ŵ<EFBFBD><C5B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>²<EFBFBD><C2B2><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>angle<6C><65>Ϊ<EFBFBD>ò<EFBFBD>angle_1;<3B><><EFBFBD><EFBFBD>fi_1ͬangle<6C><65>
|
||
% <20>˴<EFBFBD><CBB4><EFBFBD><EFBFBD>µ<EFBFBD>L_1,angle_1<5F><31>phi_1<5F><31><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĺ<EFBFBD><C4BA><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
|
||
% <20><><EFBFBD>߸նȾ<D5B6><C8BE><EFBFBD><F3A3A8B8><EFBFBD>Ui-1<><31><EFBFBD><EFBFBD>Ki0)
|
||
[Kt,Elem] = tangStiffMtxUL(Model,Anl,Elem,U) ; %%<25><><EFBFBD><EFBFBD>i-1<><31><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><CEBB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>½ڵ<C2BD>λ<EFBFBD>ú<EFBFBD><C3BA>õ<EFBFBD><C3B5>ĵ<EFBFBD>Ԫ<EFBFBD>նȾ<D5B6><C8BE><EFBFBD>Ԫ<EFBFBD><D4AA><EFBFBD>ȣ<EFBFBD><C8A3><EFBFBD>Ԫ<EFBFBD>Ƕȣ<C7B6><C8A3><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ص<EFBFBD>Elemֻ<6D><D6BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ke<4B><65><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
||
% Ԥ<><D4A4><EFBFBD><EFBFBD>λ<EFBFBD>Ƶ<EFBFBD><C6B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>delta_Ui1
|
||
d_Ut = solveLinearSystem(Model,Kt,Pref);%delta_Ui1<69><31>Ԥ<EFBFBD><D4A4><EFBFBD>ε<D7B6>λ<EFBFBD>Ƶ<EFBFBD><C6B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ǽ<EFBFBD><C7BC><EFBFBD>Ԥ<EFBFBD><D4A4><EFBFBD>εĺ<CEB5><C4BA>ر<EFBFBD><D8B1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ƿ<EFBFBD>Ϊ<EFBFBD><CEAA>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD>ĺ<EFBFBD><C4BA>ر<EFBFBD><D8B1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><CABD>ͬ<EFBFBD><CDAC>
|
||
if (step == 1)
|
||
s = 1; %<25><>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԥ<EFBFBD><D4A4><EFBFBD>ε<D7B6><CEB5><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD><C4B7><EFBFBD>Ĭ<EFBFBD><C4AC>Ϊ<EFBFBD><CEAA>
|
||
d_lr = Anl.init_incr; %delta_lamda_11 <20><>Ϊ<EFBFBD>涨Ԥ<E6B6A8><D4A4><EFBFBD><EFBFBD>d_lr
|
||
n2 = d_Ut'*d_Ut;%<25><>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><CEBB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ķ<EFBFBD><C4B7><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><C6BD>ֵ,<2C><><EFBFBD>ں<EFBFBD><DABA><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>GSP
|
||
else
|
||
% Generalized Stiffness Parameter
|
||
GSP = n2 / (d_Ut'*d_Ut_1);%<25><>ʽ<EFBFBD><CABD>41<34><31> d_Ut_1 <20><>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> d_Ut <20><>ǰ<EFBFBD><C7B0><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
||
% <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
||
if (GSP < 0)
|
||
s = -s; %<25><>ʽ<EFBFBD><CABD>42<34><32>
|
||
end
|
||
% <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><CFB5>J
|
||
% J = sqrt((Anl.trg_iter + 1) / (iter + 1));%Ŀ<><C4BF><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>/<2F><>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĵ<EFBFBD><C4B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><CABD>40<34><30><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>뽫<EFBFBD><EBBDAB><EFBFBD><EFBFBD><EFBFBD><EFBFBD> 1 <20>Ա<EFBFBD><D4B1>ⱻ<EFBFBD><E2B1BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ΪԤ<CEAA><D4A4><EFBFBD>Ľ<EFBFBD><C4BD><EFBFBD>Ӧ<EFBFBD><D3A6><EFBFBD><EFBFBD><EFBFBD>ε<EFBFBD><CEB5><EFBFBD>
|
||
J = 1;
|
||
%% Ԥ<><D4A4><EFBFBD>κ<D7B6><CEBA>ر<EFBFBD><D8B1><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD>ALG_ALCM_CYL <20>̶<EFBFBD><CCB6><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>-Բ<><D4B2>
|
||
% 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));%<25><>ʽ38<33><38>Բ<EFBFBD><D4B2><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
|
||
% <20><><EFBFBD>ر<EFBFBD>ҪС<D2AA>ڹ涨ֵ<E6B6A8><D6B5>==1<><31>
|
||
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;%<25><>֤ǡ<D6A4>ôﵽ<C3B4><EFB5BD><EFBFBD>յĺ<D5B5><C4BA>ر<EFBFBD><D8B1><EFBFBD><EFBFBD><EFBFBD>1ʵ<31>ֺ<EFBFBD><D6BA>ص<EFBFBD>ȫ<EFBFBD><C8AB>ʩ<EFBFBD><CAA9>
|
||
end
|
||
% Ԥ<><D4A4><EFBFBD><EFBFBD>λ<EFBFBD><CEBB><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
||
d_U = d_lr * d_Ut;%<25><><EFBFBD>ݹ<EFBFBD>ʽ<EFBFBD><CABD>36<33><36>Ԥ<EFBFBD><D4A4><EFBFBD>μٶ<CEBC><D9B6>в<EFBFBD>Ϊ<EFBFBD><CEAA><EFBFBD><EFBFBD><EFBFBD>ˣ<EFBFBD>ֻ<EFBFBD><D6BB>ǰһ<C7B0><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
||
% % <20>洢Ԥ<E6B4A2><D4A4><EFBFBD><EFBFBD>
|
||
d_lr_1 = d_lr;
|
||
d_Ut_1 = d_Ut;%<25>ж<EFBFBD>GSP<53><50><EFBFBD>ţ<EFBFBD><C5A3><EFBFBD>ʽ<EFBFBD><CABD>41<34><31>
|
||
d_U_1 = d_U;
|
||
% ǰ<><C7B0><EFBFBD><EFBFBD><EFBFBD>غɱȺ<C9B1>λ<EFBFBD><CEBB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ŵ<EFBFBD><C5B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϸ<EFBFBD><CFB8>µ<EFBFBD><C2B5>ӣ<EFBFBD><D3A3><EFBFBD>дD<D0B4><44><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD><D3A6>λ<EFBFBD>ƺͺ<C6BA><CDBA>ر<EFBFBD><D8B1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Сдd<D0B4><64><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD>ƺͺ<C6BA><CDBA>ر<EFBFBD><D8B1><EFBFBD><EFBFBD><EFBFBD>
|
||
D_lr = d_lr;
|
||
D_U = d_U;
|
||
% <20>ܺ<EFBFBD><DCBA>رȺ<D8B1>λ<EFBFBD><CEBB>
|
||
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);%<25><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD>ε<EFBFBD><CEB5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>õ<EFBFBD><C3B5><EFBFBD>Ke,<2C><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>D_U<5F><55>Ele.angle<6C><65>Ele.fi<66><69><EFBFBD>и<EFBFBD><D0B8><EFBFBD>
|
||
%<25><><EFBFBD><EFBFBD><EFBFBD>ǻ<EFBFBD><C7BB><EFBFBD>λ<EFBFBD>Ƽ<EFBFBD><C6BC>㣬<EFBFBD><E3A3AC>λ<EFBFBD>ƴ<EFBFBD><C6B4>ڷ<EFBFBD><DAB7><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ԣ<EFBFBD><D4A3><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><CEB5><EFBFBD><EFBFBD><EFBFBD><EFBFBD>£<EFBFBD><C2A3><EFBFBD>ţ<D7BC><C5A3><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɭ<EFBFBD><C9AD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Elem.ke<6B><65><EFBFBD><EFBFBD>
|
||
|
||
% Tangent and residual increments of displacements
|
||
d_Ut = solveLinearSystem(Model,Kt,Pref);
|
||
d_Ur = solveLinearSystem(Model,Kt,R);
|
||
|
||
% <20><><EFBFBD>ر<EFBFBD><D8B1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><CABD>45<34><35><EFBFBD><EFBFBD><EFBFBD>㶨<EFBFBD><E3B6A8><EFBFBD><EFBFBD><EFBFBD><EFBFBD>-Բ<><D4B2><EFBFBD><EFBFBD>
|
||
a = d_Ut'*d_Ut;
|
||
b = d_Ut'*(d_Ur + D_U);
|
||
c = d_Ur'*(d_Ur + 2*D_U);%D_UΪ<55><CEAA>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><D6B5>d_UrΪ<72><CEAA><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>º<EFBFBD><C2BA><EFBFBD>ֵ
|
||
s_iter = sign(D_U'*d_Ut);
|
||
d_lr = -b/a + s_iter*sqrt((b/a)^2 - c/a);
|
||
|
||
%Ϊ<>λ<EFBFBD><CEBB><EFBFBD><EFBFBD>ָ<EFBFBD><D6B8><EFBFBD>
|
||
if (~isreal(d_lr))
|
||
conv = -1;
|
||
break;
|
||
end
|
||
|
||
% <20><>ʽ<EFBFBD><CABD>36<33><36>
|
||
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
|
||
|
||
|
||
|