Week - 2 - Explicit and Implicit Analysis
Explicit and Implicit Analysis:
Aim: To use implicit and explicit analysis methods numerically to solve a simple problem for better understanding and drawing conclusion for which method is better in what case.
Solution:
Both Implicit and explicit methods are widely used for numerical analysis such as in FEA. Both have their own pro’s and con’s.
Implicit method: Widely used in static problems and modal analysis and other simpler problems in FEA where the deformations are generally smaller. It involves inversion of mass and stiffness matrix and hence is avoided in complex problems where a lot of nodes and elements are involved as it will make things complicated and costly:
u(n+1) = u(n) + Δt*f(u(n+1)), u(n) = (I-A)* Δt*u(n+1)
In implicit method problems we are only interested in final answers and not the changes occurring till the final answer. As a result the time steps are larger- order of seconds, and there are no conditions for the time step for a stable analysis and hence this method is unconditionally stable.
Implicit method is used wherever geometric complexity or scale of problem is too much to be solved numerically by hand which is usually in case of linear static problems of beams with different c/s, heat transfer over a complex shape etc. It can be used for other types of problems which are dynamic, but is not suitable and efficient due to the reasons mentioned above.
Explicit Method: Finds application in crashworthiness, drop test, and other such dynamic analysis where the deformations are very large. Contrary to the implicit method, there is no inversion of stiffness matrix and hence is more suitable for complex problems with number of components and elements.
u(n+1) = u(n) + Δt*f(u(n))
Unlike implicit method, this method is heavily dependent on the time step for not only accuracy of the solutions but for the stability of the analysis too. The condition for the time step is it should be smaller than the time it takes for the shock-wave which travels at speeds equivalent to the speed of sound through the given medium (material) to propagate across an element and is given by Courant-Freidrichs-Levy Criterion :
where, ‘L’ is critical length, ‘ρ’ is density if the material,’ µ’ is the modulus if rigidity, and ‘v’ is the effective kinematic viscosity.
This makes the explicit method a ‘conditionally stable’ analysis method as the stability of the analysis depends crucially on time step and hence the time steps are very small – order of ‘ms’. Large element processing and small time steps endure the analysis is stable and thorough. In non-linear dynamic problem such as crash test, the final result and the accuracy of the output are not much important, but the process it follows till the end is. As a result the accuracy of results is measured by the distribution of results over time steps and as a result explicit method can be used effectively.
Problem:
F(u) = u3+9u2+4u
use this equation and solve using both Explicit and Implicit Methods ( have a tolerance of 10-2)
Explicit Method:
let u0 = 0 be the initial displacement of the rod.
f(u) =u3+9u2+ 4u –given
delta_F = 1 – External force step.
f'(u) = k(u) = 3u2+18u+4
k(u) being the stiffness of the rod
Step 1:
u1= u0 + delta_u1
delta_u1 being the change in ‘u’ due to the external force
delta_u = delta_F / k(u)
f(0) = 0, k(0) = 4
u(0)=0
delta_u1 = 1/4 = 0.25
u1 = 0.25
f_int = f(u1) = 1.578125
Residual = f_int - F_ext = 0.578125
Similarly we can calculate u2, u3.....u6 with an incremental external load ‘delta_F = 1’ at each step.
2.
f(u1) = 1.578125
k(u1) = 8.6875
u1= 0.25
delta_u1 = 1/8.6875 = 0.115107914
u2 = u1 + delta_u2, u2 = 0.365107914
f_int = f(u2) = 2.708836, F_ext = 1+1=2
Residual = f_int - F_ext = 0.70883602
3.
f(u2) = 2.708836
k(u2) = 10.97185
u2 = 0.365107914
delta_u3 = 1/10.97185 = 0.091142
u3 = u2 + delta_u3, u3 = 0.45625
f_int = f(u3) = 3.793454, F_ext = 1+1+1= 3
Residual = f_int - F_ext = 0. 793454
The procedure is repeated for further steps:
| f(u)_n-1 | k(u)_n-1 | u_n-1 | delta_u(n) | u_n | Residual f_int-F_ext |
| 0 | 4 | 0 | 0.25 | 0.25 | 0.578125 |
| 1.578125 | 8.6875 | 0.25 | 0.115108 | 0.365108 | 0.70883602 |
| 2.708836 | 10.97185 | 0.365108 | 0.091142 | 0.45625 | 0.793454169 |
| 3.793454 | 12.837 | 0.45625 | 0.0779 | 0.53415 | 0.856848469 |
| 4.856848 | 14.47065 | 0.53415 | 0.069105 | 0.603255 | 0.907811078 |
| 5.907811 | 15.95035 | 0.603255 | 0.062695 | 0.66595 |
Table of results:
|
u_explicit |
f(u_explicit) |
|
0 |
0 |
|
0.25 |
1.578125 |
|
0.365108 |
2.708836 |
|
0.45625 |
3.7934542 |
|
0.53415 |
4.8568485 |
|
0.603255 |
5.9078111 |
|
0.66595 |
6.9505464 |
Implicit Method: In this method, we use Newton-rhapson method to achieve equilibrium of the forces such that, f_internal = F_external
Similar to process followed in implicit method:
let u0 = 0 be the initial displacement of the rod.
f(u) =u3+9u2+ 4u –given
delta_F = 1 – External force step.
f'(u) = k(u) = 3u2+18u+4
k(u) being the stiffness of the rod
f(0) = 0, k(0) = 4
u(0)=0
delta_u1 = 1/4 = 0.25
u1 = 0.25
f_int = f(u1) = 1.578125
Residual = f_int - F_ext = 0.578125
Tolerance check = Residual < 0.01, False.
As residual >0.01, we perform Newton-Rhapson iterations to reach this tolerance by making changes to f(u) and hence 'u'.
Using Taylor series for the residuals function R(u) = f_int - F_ext for current displacement term 'u_n',
R_n+1 = R_n(u_n) + d/du(R(u_n))*delta(u_n+1) = R_n + k(u_n)*delta(u_n+1)
putting R_n+1 = 0, neglecting the higher order terms, and solving for delta(u_n+1),
delta(u_n+1) = - R_n*1/[k(u_n)]
This term can be used as the 'correction factor' to get a more appropriate figure for 'u_n'. We can repeat this for a 2nd iteration if the corrected displacement doesn't fullfill the tolerance condition over and over.
1st iteration:
u1_0 = 0.25
k(u1_0) = 8.6875
Residual = f_int - F_ext = R_n = 0.578125
u1_1 = u1_0 + delta(u1_0)
delta(u1_0) = - R_0*1/[k(u1_0)]
delta(u1_0) = - (0.578125 * 1/8.6875) = - 0.066546763
Therefore, u1_1 = 0.25 + (- 0.066546763) = 0.183453237
f(u1_1) = 1.042883, R_1 = 1.042883 – 1 = 0.042883 > 0.01. Hence, we need to iterate once more.
2nd iteration:
u1_2 = 0.183453237, k(u1_2) = 7.403124
u1_2 = u1_0 + delta(u_2)
delta(u1_2) = u1_1 - R_1*1/[k(u1_1)] = 0.183453237 – (0.042883 /7.403124)
delta(u1_2) = - 0.072339304
Hence,
u1_2 = 0.177661, f(u1_2) = 1.00032, R_2 = 1.00032 – 1 = 0.000320254 < 0.01
f(u1_2) = 1.00032
k(u1_2) = 7.403124
u1_0 = 0.25
delta(u1_2) = - 0.072339
u1_2 = 0.177661
R1_2 = 0.000320254
This procedure is then repeated for several steps.
| 1 | f(u)_n-1 | k(u)_n-1 | u_n-1 | delta_u(n) | u_n | Residual = f_int-F_ext | Tolerence check | ||
| u1 = ? | 0 | 4 | 0 | 0.25 | 0.25 | 0.578125 | FALSE | ||
| u1= u0+delta_u1 | |||||||||
| delta_u = delta_F / k(u) | |||||||||
| delta_F | 1 | ||||||||
| f_int = | f(u_n)= | 1.578125 | |||||||
| 1st iteration: | |||||||||
| u1_0 = | 0.25 | k(u1_0) = | 8.6875 | ||||||
| u1_1 =? | |||||||||
| u1_1 = u1_0 + delta(u1_0) | |||||||||
| delta(u_1) = - R_0*1/[k(u1_0)] | |||||||||
| delta(u_1) = | -0.066546763 | u1_1 = | 0.183453 | ||||||
| Now: | 1.042883 | 0.042882898 | FALSE | ||||||
| 2nd iteration: | |||||||||
| u1_1 = | 0.183453237 | k(u1_1) = | 7.403124 | ||||||
| u1_2 = u1_0 + delta(u_2) | |||||||||
| delta(u_2) = u1_1 - R_1*1/[k(u1_1)] | |||||||||
| delta(u_2) = | -0.072339304 | u1_2 = | 0.177661 | ||||||
| Now: | 1.00032 | 0.000320254 | TRUE | ||||||
| Hence, | |||||||||
| 1.00032 | 7.403124 | 0.25 | -0.072339 | 0.177661 | 0.000320254 | TRUE | |||
| 2 | u2 = ? | 1.00032 | 7.292583 | 0.177661 | 0.137126 | 0.314786 | 0.18215154 | FALSE | |
| u2= u1+delta_u2 | |||||||||
| delta_u2 = delta_F / k(u) | |||||||||
| delta_F | 1 | ||||||||
| f_int = | f(u_n)= | 2.182152 | |||||||
| 1st iteration: | |||||||||
| u2_0 = | 0.31478633 | k(u2_0) = | 9.963425 | ||||||
| u2_1 =? | |||||||||
| u2_1 = u2_0 + delta(u2_0) | |||||||||
| delta(u2_0) = - R2_0*1/[k(u2_0)] | |||||||||
| delta(u2_0) = | -0.01828202 | u2_1 = | 0.296504 | ||||||
| Now: | 2.003318 | 0.003317615 | TRUE | ||||||
| Hence, | |||||||||
| 2.003318 | 9.963425 | 0.177661 | -0.018282 | 0.296504 | 0.003317615 | TRUE | |||
| 3 | u3 = ? | 2.003318 | 9.600822 | 0.296504 | 0.104158 | 0.400662 | 0.111737314 | FALSE | |
| u3= u2+delta_u3 | |||||||||
| delta_u3 = delta_F / k(u) | |||||||||
| delta_F | 1 | ||||||||
| f_int = | f(u_n)= | 3.111737 | |||||||
| 1st iteration: | |||||||||
| u3_0 = | 0.400662058 | k(u3_0) = | 11.69351 | ||||||
| u3_1 =? | |||||||||
| u3_1 = u3_0 + delta(u3_0) | |||||||||
| delta(u3_0) = - R3_0*1/[k(u3_0)] | |||||||||
| delta(u3_0) = | -0.0095555 | u3_1 = | 0.391107 | ||||||
| Now: | 3.000931 | 0.000930646 | TRUE | ||||||
| Hence, | 3.000931 | 11.69351 | 0.296504 | -0.009556 | 0.391107 | 0.000930646 | TRUE | ||
| 4 | u4 = ? | 3.000931 | 11.49881 | 0.391107 | 0.086966 | 0.478072 | 0.078529187 | FALSE | |
| u4= u3+delta_u4 | |||||||||
| delta_u4 = delta_F / k(u) | |||||||||
| delta_F | 1 | ||||||||
| f_int = | f(u_n)= | 4.078529 | |||||||
| 1st iteration: | |||||||||
| u4_0 = | 0.478072071 | k(u4_0) = | 13.29096 | ||||||
| u4_1 =? | |||||||||
| u4_1 = u4_0 + delta(u4_0) | |||||||||
| delta(u4_0) = - R4_0*1/[k(u4_0)] | |||||||||
| delta(u4_0) = | -0.005908468 | u4_1 = | 0.472164 | ||||||
| Now: | 4.000364 | 0.000364052 | TRUE | ||||||
| Hence, | 4.000364 | 13.29096 | 0.391107 | -0.005908 | 0.472164 | 0.000364052 | TRUE | ||
| 5 | u5 = ? | 4.000364 | 13.16776 | 0.472164 | 0.075943 | 0.548107 | 0.060877573 | FALSE | |
| u5= u4+delta_u5 | |||||||||
| delta_u5 = delta_F / k(u) | |||||||||
| delta_F | 1 | ||||||||
| f_int = | f(u_n)= | 5.060878 | |||||||
| 1st iteration: | |||||||||
| u5_0 = | 0.548106662 | k(u5_0) = | 14.76718 | ||||||
| u5_1 =? | |||||||||
| u5_1 = u5_0 + delta(u5_0) | |||||||||
| delta(u5_0) = - R5_0*1/[k(u5_0)] | |||||||||
| delta(u5_0) = | -0.004122491 | u5_1 = | 0.543984 | ||||||
| Now: | 5.000181 | 0.000180829 | TRUE | ||||||
| Hence, | 5.000181 | 14.76718 | 0.472164 | -0.004122 | 0.543984 | 0.000180829 | TRUE | ||
| 6 | u6 = ? | 0 | 14.67947 | 0.543984 | 0.068122 | 0.612107 | 0.049836172 | FALSE | |
| u6= u5+delta_u6 | |||||||||
| delta_u6 = delta_F / k(u) | |||||||||
| delta_F | 1 | ||||||||
| f_int = | f(u_n)= | 6.049836 | |||||||
| 1st iteration: | |||||||||
| u6_0 = | 0.612106515 | k(u6_0) = | 16.14194 | ||||||
| u6_1 =? | |||||||||
| u6_1 = u6_0 + delta(u6_0) | |||||||||
| delta(u6_0) = - R5_0*1/[k(u6_0)] | |||||||||
| delta(u6_0) = | -0.003087372 | u6_1 = | 0.609019 | ||||||
| Now: | 6.000103 | 0.000103261 | TRUE | ||||||
| Hence, | 6.000103 | 16.14194 | 0.543984 | -0.003087 | 0.609019 | 0.000103261 | TRUE | ||
Table of the implicit method results:
|
u_implicit |
f(u_implicit) |
|
0 |
0 |
|
0.177661 |
1.00032 |
|
0.296504 |
2.003318 |
|
0.391107 |
3.000931 |
|
0.472164 |
4.000364 |
|
0.543984 |
5.000181 |
|
0.609019 |
6.000103 |
Plot of the acquired results:
Using the equation for internal force in the rod f(u), and inputting various values for ‘u’ 0 to 6.6, the exact solution table is acquired and used to compare the results.
Conclusions from the plot:
- The equation for internal force generated in a rod is nonlinear (Cubic) in nature and hence the exponential plot. The analysis is done using an incremental external force.
- Due to the small sample space, and relative simplicity of the problem, it is NOT VERY evident but the explicit analysis seems to drift a bit from the exact solution. Implicit analysis with Newton-Raphson iterations helps reaching equilibrium between internal and external forces and hence follows close to the exact solution. Even though the difference is very minute, a more complicated problem of this type would favour implicit method.
- The concepts and methods used in this problem are similar to one's used in a nonlinear finite element analysis. Therefore, the given example is a useful learning tool for nonlinear problems.