Solving a simple problem using both Explicit and Implicit Analysis
OBJECTIVE: To solve the given problem using both explicit and implicit methods and compare using the Force v/s Displacement plots.
PROBLEM: Consider the case of a simple bar in tension as shown in the Figure. Suppose that the force in the bar is a nonlinear function of the displacement, `F(u) =(u)^3 + 9(u)^2 + 4u`.
From the force function, the stiffness of the bar is `k(u) =d/(du)F(u)=3u^2 + 18u + 4`
Clearly it is possible to plot a graph of force, f, versus displacement, u, by using the force equation. However, if the nonlinear force versus displacement relationship is not known, but the stiffness is known, it is possible to construct the graph of force versus displacement numerically. One method is to use an incremental explicit load control scheme wherein, equal load increments (or steps) of dF = 1 unit each are used to load the tension bar. Similarly, an implicit method can also be used to solve the problem by using Newton-Raphson iterations to establish convergence after each step. Solve the problem using both methods and plot the Force v/s displacement graph for both methods.
THE EXPLICIT AND IMPLICIT SCHEMES OF ANALYSIS
The explicit and implicit methods are simply two different approaches used to solve ODEs and PDEs used to depict physical processes as computer simulations. While the Explicit scheme uses the value from the present time step to calculate the value in the next timestep, the implicit scheme calculates the value in the next time step by using the present values and expressing the source term as a function of the next time step. This creates a need for Matrix inversion in the case of the implicit method. Having said all this, we must note that these schemes are just two different ways of solving the same problem. A fundamental difference between the two schemes is that the Explicit scheme assumes convergence for all time steps and thus doesn't check for convergence after each timestep. On the other hand, the implicit scheme checks for convergence after every timestep before proceeding to the next. This makes the implicit scheme stable but demanding in terms of computations as most of the times, convergence is reached only after several iterations. As a result, the time taken for computing for a single time step is bigger. The explicit scheme can be considered as more direct as there are no convergence checks done in between but this generates stability problems. To avoid this, the timesteps have to be kept extremely small and in most cases, the user doesn't have much control over the selection of this extremely small timestep. This is not the case for the implicit method. Since the convergence check is performed after every time step, the user has control over how big each timestep should be, in the implicit method.
USING THE EXPLICIT SCHEME
Force function, `F(u) =(u)^3 + 9(u)^2 + 4u`
Stiffness, `k(u) =d/(du)F(u)=3u^2 + 18u + 4`
Displacement = Force/Stiffness`=F/k`
Force increment`= dF = 1`
Displacement increment, `du=(dF)/(k)`
Displacement, ` u_(i+1) = u_i+du_i`
Step 1:
`dF=1`
`k(u_0)= 3(u_0)^2 + 18u_0 + 4 = 3*0^2+18*0+4=4`
`du_1= (dF)/(k(u_0))= 1/4=0.25`
`u_1 = u_0 + du_1=0+0.25=0.25`
`F(u_1) =(u_1)^3 + 9(u_1)^2 + 4u_1=0.25^3+9*0.25^2+4*0.25=1.578`
`F_(ext)=dF=1`
`u_(exact) =0.177` { by solving the equation, `F_(ext) =(u_(exact))^3 + 9(u_(exact))^2 + 4u_(exact)`}
Step 2:
`dF=1`
`k(u_1)= 3(u_1)^2 + 18u_1 + 4 = 3*0.25^2+18*0.25+4=8.688`
`du_2= (dF)/(k(u_1))= 1/8.688=0.115`
`u_2 = u_1 + du_2=0.25+0.115=0.365`
`F(u_2) =(u_2)^3 + 9(u_2)^2 + 4u_2=0.365^3+9*0.365^2+4*0.365=2.709`
`F_(ext)=dF+dF=1+1=2`
`u_(exact) =0.296` { by solving the equation, `F_(ext) =(u_(exact))^3 + 9(u_(exact))^2 + 4u_(exact)`}
Step 3:
`dF=1`
`k(u_2)= 3(u_2)^2 + 18u_2 + 4 = 3*0.365^2+18*0.365+4=10.972`
`du_3= (dF)/(k(u_2))= 1/10.972=0.091`
`u_3 = u_2 + du_3=0.365+0.091=0.456`
`F(u_3) =(u_3)^3 + 9(u_3)^2 + 4u_3=0.456^3+9*0.456^2+4*0.456=3.793`
`F_(ext)=dF+dF+dF=1+1+1=3`
`u_(exact) =0.391` { by solving the equation, `F_(ext) =(u_(exact))^3 + 9(u_(exact))^2 + 4u_(exact)`}
| Step i | ui | dF | k(ui) | dui= dF/k(ui) | ui+1 | Fext | Fint | R= Fint - Fext | uexact |
| 0 | 0 | 0 | 4 | 0 | 0.25 | 0 | 0 | 0 | 0 |
| 1 | 0.25 | 1 | 8.688 | 0.25 | 0.365 | 1 | 1.578 | 0.578 | 0.177 |
| 2 | 0.365 | 1 | 10.972 | 0.115 | 0.456 | 2 | 2.709 | 0.709 | 0.296 |
| 3 | 0.456 | 1 | 12.837 | 0.091 | 0.534 | 3 | 3.793 | 0.793 | 0.391 |
USING IMPLICIT SCHEME
Unlike the Explicit scheme which assumes the Equilibrium to exist at all times and directly moves forward to calculate the variables, the implicit scheme checks convergence after each time step and then only moves into the calculation of the variables. The convergence is brought about by several Newton-Raphson iterations. For convergence, the `F_(ext)` should be equal to `F` or the Residual, `R = F-F_(ext) = 0`. After each timestep, the residual need to be brought within the tolerance value(here, we are assuming the tolerance as 0.01). This is done by using the Newton-Raphson iteration to calculate the correction values until the required tolerance is reached.
Expanding the residual, R as a Taylor series about the current displacement `u^j` gives:
`R ^(j+1) = R ^j (u^ j ) + (dR(u ^j ) )/(du) δu^(j+1) + ... = R ^j + k(u ^j )δu^(j+1) + ...`
By neglecting higher-order terms in the series, setting `R^(j+1) = 0` and solving for `δu^(j+1)`, the following correction is obtained :
`δu^(j+1) = −[k(u^ j )]^(−1)R ^j`
Notice that the iteration variable is `j` and it may take several iterations for the norm of the residual, `||R||`, to be reduced below the chosen tolerance criteria(0.01 in this case). The implicit analysis proceeds as follows:
Force function, `F(u) =(u)^3 + 9(u)^2 + 4u`
Stiffness, `k(u) =d/(du)F(u)=3u^2 + 18u + 4`
Displacement = Force/Stiffness`=F/k`
Force increment`= dF = 1`
Displacement increment, `du=(dF)/(k)`
Displacement, ` u_(i+1) = u_i+du_i`
Step 1:
`dF=1`
`k(u_0)= 3(u_0)^2 + 18u_0 + 4 = 3*0^2+18*0+4=4`
`du_1= (dF)/(k(u_0))= 1/4=0.25`
`u_1 = u_0 + du_1=0+0.25=0.25`
`F(u_1) =(u_1)^3 + 9(u_1)^2 + 4u_1=0.25^3+9*0.25^2+4*0.25=1.578`
`F_(ext)=dF=1`
`R = F(u_1)-F_(ext) = 1.578 - 1.000 = 0.578`
`||R|| > 0.01`. Therefore, Newton-Raphson iterations have to be done to obtain a correction factor that will bring the `R` value within our tolerance limit.
| No. | ITERATIONS | |
|---|---|---|
| 1. | `u_1^((0)) = u_1= 0.25` | |
|
Correction factor for `u_1= delta u_1^((1)) = −[k(u_1^ ((0)) )]^(−1)R^((0))` `delta u_1^((1)) = −[k(u_1^ ((0)) )]^(−1)R^((0)) = -[8.688]^(-1)*0.578 = -0.06653` |
`k(u_1^ ((0)) ) = k(u_1 ) =3*0.25^2+18*0.25+4= 8.688` `R^((0)) = 0.578` |
|
| Corrected `u_1 = u_1^((1)) = u_1 + delta u_1^((1)) = 0.25 + (-0.06653) = 0.183` | ||
|
Updated Residual, `R^((1)) = F^((1))-F_(ext) = 1.039-1 = 0.039` |
`F^((1)) = (u_1^((1)))^3 + 9(u_1^((1)))^2 + 4u_1^((1))` `=0.183^3+9*0.183^2+4*0.183= 1.039` |
|
| `||R||>0.01`. Proceed to the next iteration as `||R||` is not within the tolerance limit. | ||
| 2. |
Correction factor for `u_1= delta u_1^((2)) = delta u_1^((1))−[k(u_1^ ((1)) )]^(−1)R^((1))` `delta u_1^((2)) = delta u_1^((1))−[k(u_1^ ((1)) )]^(−1)R^((1)) =-0.06653 -[7.394]^(-1)*0.039 = -0.0718` |
`k(u_1^ ((1)) ) = 3*(0.183)^2 + 18*(0.183) + 4 =7.394` `R^((1)) = 0.039` |
| Corrected `u_1 = u_1^((2)) = u_1 + delta u_1^((2)) = 0.25 + (-0.0718) = 0.1782` | ||
| Updated Residual, `R^((2)) = F^((2))-F_(ext) = 1.004-1 = 0.004` |
`F^((2)) = (u_1^((2)))^3 + 9(u_1^((2)))^2 + 4u_1^((2))` `=0.1782^3+9*0.1782^2+4*0.1782= 1.004` |
|
| `||R||<0.01`. `||R||` has reached within the tolerance limit. Thus, no need for further iterations | ||
Therefore, `F(u_1) = F(u_1^((2)))= 1.004` and the final corrected value of `u_1= u_1^((2)) = 0.1782`
Step 2:
`dF=1`
`k(u_1)= 3(u_1)^2 + 18u_1 + 4 = 3*0.1782^2+18*0.1782+4= 7.302`
`du_2= (dF)/(k(u_1))= 1/7.302=0.137`
`u_2 = u_1 + du_2=0.1782+0.137=0.315`
`F(u_2) =(u_2)^3 + 9(u_2)^2 + 4u_2=0.315^3+9*0.315^2+4*0.315=2.184`
`F_(ext)=dF+dF=1+1 = 2`
`R = F(u_2)-F_(ext) = 2.184 - 2.000 = 0.184`
`||R|| > 0.01`. Therefore, Newton-Raphson iterations have to be done to obtain a correction factor that will bring the `R` value within our tolerance limit.
| No. | ITERATIONS | |
|---|---|---|
| 1. | `u_2^((0)) = u_2` | |
|
Correction factor for `u_2= delta u_2^((1)) = −[k(u_2^ ((0)) )]^(−1)R^((0))` `delta u_2^((1)) = −[k(u_2^ ((0)) )]^(−1)R^((0)) = -[9.967]^(-1)*0.184 = -0.0185` |
`k(u_2^ ((0)) ) = k(u_2 ) = 9.967` `R^((0)) = 0.184` |
|
| Corrected `u_2 = u_2^((1)) = u_2 + delta u_2^((1)) = 0.315 + (-0.0185) = 0.296` | ||
|
Updated Residual, `R^((1)) = F^((1))-F_(ext) = 1.998-2 = -0.002` |
`F^((1)) = (u_2^((1)))^3 + 9(u_2^((1)))^2 + 4u_2^((1))` `=0.296^3+9*0.296^2+4*0.296= 1.998` |
|
| `||R||<0.01`. `||R||` has reached within the tolerance limit. Thus, no need for further iterations | ||
Therefore, `F(u_2) = F(u_2^((1)))= 1.998` and the final corrected value of `u_2= u_2^((1)) = 0.296`
Step 3:
`dF=1`
`k(u_2)= 3(u_2)^2 + 18u_2 + 4 = 3*0.296^2+18*0.296+4= 9.59`
`du_3= (dF)/(k(u_2))= 1/9.59=0.104`
`u_3 = u_2 + du_3=0.296+0.104=0.4`
`F(u_3) =(u_3)^3 + 9(u_3)^2 + 4u_3=0.4^3+9*0.4^2+4*0.4=3.104`
`F_(ext)=dF+dF+dF=1+1+1 = 3`
`R = F(u_3)-F_(ext) = 3.104 - 3.000 = 0.104`
`||R|| > 0.01`. Therefore, Newton-Raphson iterations have to be done to obtain a correction factor that will bring the `R` value within our tolerance limit.
| No. | ITERATIONS | |
|---|---|---|
| 1. | `u_3^((0)) = u_3 = 0.4` | |
|
Correction factor for `u_3= delta u_3^((1)) = −[k(u_3^ ((0)) )]^(−1)R^((0))` `delta u_3^((1)) = −[k(u_3^ ((0)) )]^(−1)R^((0)) = -[11.68]^(-1)*0.104 = -0.0089` |
`k(u_3^ ((0)) ) = k(u_3 ) =3*0.4^2+18*0.4+4= 11.68` `R^((0)) = 0.104` |
|
| Corrected `u_3 = u_3^((1)) = u_3 + delta u_3^((1)) = 0.4 + (-0.0089) = 0.3911` | ||
|
Updated Residual, `R^((1)) = F^((1))-F_(ext) = 3.0008-3 = 0.0008` |
`F^((1)) = (u_3^((1)))^3 + 9(u_3^((1)))^2 + 4u_3^((1))` `=0.3911^3+9*0.3911^2+4*0.3911= 3.0008` |
|
| `||R||<0.01`. `||R||` has reached within the tolerance limit. Thus, no need for further iterations | ||
Therefore, `F(u_3) = F(u_3^((1)))= 3.0008` and the final corrected value of `u_3= u_3^((1)) = 0.3911`
| Step i | Displacement | dF | k(ui) | dui= dF/k(ui) | ui+1 | Fext | Force | Residual | uexact | |||
| ui | Corrected Ui | Fint | Corrected Fint | R= Fint - Fext | Corrected R | |||||||
| 0 | 0.000 | 0.000 | 0 | 4.000 | 0.000 | 0.250 | 0 | 0.0000 | 0.0000 | 0.000 | 0.0000 | 0 |
| 1 | 0.250 | 0.178 | 1 | 8.688 | 0.115 | 0.365 | 1 | 1.5781 | 1.0043 | 0.578 | 0.0043 | 0.177 |
| 2 | 0.365 | 0.296 | 1 | 10.972 | 0.091 | 0.456 | 2 | 2.7088 | 1.9985 | 0.709 | -0.0015 | 0.296 |
| 3 | 0.456 | 0.391 | 1 | 12.837 | 0.078 | 0.534 | 3 | 3.7935 | 3.0009 | 0.793 | 0.0009 | 0.391 |
RESULT and DISCUSSION:
- The results from both the Explicit and the implicit method has been listed alongside the exact solution in the table below.
Steps Force Displacement Explicit Implicit Exact Step 0 0 0 0 0 Step 1 1 0.25 0.178 0.177 Step 2 2 0.365 0.296 0.296 Step 3 3 0.456 0.391 0.391 - A graph has also been plotted to visualize the variation in displacements calculated by the two methods in comparison with the exact displacements. The displacements have been plotted in the X-axis and the Y-axis shows the force applied in incremental steps of 1 unit.
- It can be observed that the explicit solutions are deviating as more and more load is being applied. On the other hand, the implicit solution and exact solution are overlapping each other perfectly. The Newton-Raphson iterations have corrected the incremental steps so that they land on the exact solution according to the specified tolerance.
- The table below contains the explicit results when the load increment per step is reduced to 0.01 units instead of 1 unit, i.e, the total load of 3 units is applied in 24 steps(each with dF= 0.01 unit) instead of just 3 steps. The Force vs displacement plot for this case is as shown below
Step i ui dF k(ui) dui= dF/k(ui) ui+1 Fext Fint R= Fint - Fext 0 0 0 4.000 0.003 0.003 0 0 0 1 0.003 0.01 4.045 0.005 0.007 0.01 0.010056 5.63E-05 2 0.007 0.02 4.134 0.007 0.015 0.03 0.030277 0.000277 3 0.015 0.03 4.265 0.009 0.024 0.06 0.060752 0.000752 4 0.024 0.04 4.435 0.011 0.035 0.1 0.101548 0.001548 5 0.035 0.05 4.640 0.013 0.048 0.15 0.152703 0.002703 6 0.048 0.06 4.876 0.014 0.063 0.21 0.214227 0.004227 7 0.063 0.07 5.139 0.016 0.078 0.28 0.286115 0.006115 8 0.078 0.08 5.426 0.017 0.095 0.36 0.368345 0.008345 9 0.095 0.09 5.733 0.017 0.112 0.45 0.46089 0.01089 10 0.112 0.1 6.058 0.018 0.130 0.55 0.56372 0.01372 11 0.130 0.11 6.398 0.019 0.149 0.66 0.676805 0.016805 12 0.149 0.12 6.751 0.019 0.168 0.78 0.800115 0.020115 13 0.168 0.13 7.116 0.020 0.188 0.91 0.933625 0.023625 14 0.188 0.14 7.491 0.020 0.208 1.05 1.077311 0.027311 15 0.208 0.15 7.876 0.020 0.228 1.2 1.231154 0.031154 16 0.228 0.16 8.268 0.021 0.249 1.36 1.395134 0.035134 17 0.249 0.17 8.668 0.021 0.270 1.53 1.569237 0.039237 18 0.270 0.18 9.074 0.021 0.291 1.71 1.75345 0.04345 19 0.291 0.19 9.486 0.021 0.312 1.9 1.94776 0.04776 20 0.312 0.2 9.903 0.021 0.333 2.1 2.152158 0.052158 21 0.333 0.21 10.326 0.021 0.354 2.31 2.366635 0.056635 22 0.354 0.22 10.754 0.021 0.376 2.53 2.591183 0.061183 23 0.376 0.23 11.185 0.021 0.397 2.76 2.825797 0.065797 24 0.397 0.24 11.621 0.022 0.419 3 3.070469 0.070469 - The graph shows that the new Explicit plot (Explicit2) lies very close to the exact solution plot. The application of the same load in a greater number of steps, by reducing the load increment per step to 0.01 units has reduced the deviation of the explicit results. This depicts the fact that the timesteps have to be kept extremely small for explicit analysis to achieve convergence. This limitation is due to the absence of any convergence check after each time step.
CONCLUSION:
- The problem was solved using both methods.
- The results were analyzed using the Force vs displacement plots.
- The explicit solutions were found to have deviated from the exact solution whereas the implicit solutions coincided with the exact solutions.
- Furthermore, the problem was solved again using explicit method, only this time, the load increment per step was reduced to 0.01. The new explicit results are closer to the exact solution than before, i.e, the deviation has reduced significantly, thus confirming the need to keep the timestep as extremely small for explicit schemes.