4. STEPS OF NEEDLE GRIPPER DESIGN FOR DEMOULDING FLEXIBLE
4.2 Determining the pull-out force for one needle and its stress zeroing radius
4.2.2 Solving the initial value problems for the differential equations in Chapter 4.2.1 and optimising the
4.2.2.2 Numerical solution for the initial values
To obtain a particular solution for a second order ODE two initial conditions have to be given and with them the Cauchy problem have to be solved. These values can be the value of the function at a given point and the initial slope of the integral curve at the same point. To determine these values physical considerations will be used. In the function
r r (r)
) (
Q = ρ and
0 0
0)or (0) ) (
0 (
Q r r
r ρ
= ρ (107)
as equation (51) in Chapter 4.2.1.2 only R0 is known that is the radius of the inserted needle.
The initial value of ρ0 is not known at this moment, the condition for the calculations assumes that ρ(0) < r0. At the end of this Chapter a method for determining the value of ρ0 theoretically will be shown.
For determining the second constant for the initial value problem the following assumption will be used. The effect of the needle is considered local, this way the disturbance of the enlargement of ρ0 to R0 only acts at the vicinity of the needle. This effect cannot be detected at large distances and at infinity it disappears. The formulation of this assumption is:
1
) 1 R(
lim ) q(
Q(r) lim lim
r r
r = =
∞
→
∞
→
∞
→ r
r (108)
this is the second initial condition. The problem with these conditions that the two values are given at different points (r=0 and r=∞) of the manifold. Because equation (97) can only be solved numerically these conditions cannot be implied, because the numerical values cannot be determined at infinity. The problem has to be transformed to an initial value problem, where the two conditions are set at the same point of the manifold. In the new set of problems both values will be given at r=0 the way that the original conditions are satisfied:
r W Q(0) ρ(0)
0
=
= and Q′(0)=K (109)
where:
dr dQ(r) Q
K= 0′ = , the value of the derivative of Q(r) at 0. The 0 point will be set at the boundary of the inserted 1 mm diameter needle as in the calculation, this way r0=0.5 and it is measured from the axis of the needle.
r(0.5) W ρ(0.5)
Q(0.5)= = and Q′(0.5)=K (110)
The value of K cannot be chosen arbitrarily, because condition (110) has to be in correlation with (107) and (108). That is W and K has to be chosen the way that also satisfies condition (108). For different needles different W and K has to be determined. The problem is that K cannot be determined directly from a physical point of view, only indirectly from (108). To determine K an iterative program method will be used.
As mentioned earlier the numerical methods have problems when dealing with infinity. To handle the problem a very large number (T) shall be chosen (e.g.: T= 5·105 which is >>R0) and the behaviour of the integral curve will be analysed at that point of the manifold. The rou-tine that calculates the approximate value of K is achieved by an interval splitting method2 . In this method two different values (a1, a2) of K will be chosen arbitrarily (from the previous Chapter we know that the value ranges somewhere between 0.2 and 30, but for the first run we can choose -100 and 100 to be surely in the range). A third value is determined as the mean value of the two a3= (a1+a2)/2. Equation (97) is solved independently for the three val-ues for K (and of course for W). The numerical method produces the routine for the calcula-tions of Q(r) at discreet points of the manifold.
Let us determine Q(T) which can be considered as the value of Q at infinity. The next step is to subtract one from the value of Q(T), this is
sol1= Q(T)1-1, (111)
this is done for the other two solutions as well (sol2, sol3) and analysed whether they are negative or positive. These values then can be compared with each other. The final solution would be that solx=0, but the equations will be solved for a given error.
For the cycle a small number will be chosen which will set the size of the error (this can be e.g.:10-8). Two of the three solutions will be multiplied together and analysed whether they are larger or smaller than zero.
If sol2·sol3>0 (112)
then the solution can definitely be found in the other interval.
If sol1·sol3<0 (113)
then the solution can be found inside this interval. (If the solution is not negative in either case, than different ai-s have to be chosen). After sol1·sol3<0 is achieved then this interval shall be split into two, by a1 staying the same and a2 becoming a2new and a3 becoming a3=(a1+a2new)/2. This way the interval is split into two and the cycle can run again, until the
2 Originally programmed by Dr. Laczik Bálint
differences become less than 10-8. With the gained values K will become a3, and with this the initial condition is solved. With the operators of Maple the values of Q(r) can be calculated for the arbitrary values of r (as long as r ≥ R0).
This routine converges fast, so if no singular solution is present between the intervals the so-lution will be found in 30-40 steps.
The start values used in the program are a1=0.2, a2=100. The calculation takes a long time, because the ODE have to be solved numerically, then the values have to be substituted into Q(T) and the results have to be compared (according to equations: (111), (112), (113) in every cycle. (For a Pentium 4 (2.4 GHz) computer it takes approximately 10-30 seconds of computation time).
The numerical solution is given in Maple by the following procedures:
:=
sol1 proc (x_rkf45) ... end proc
:=
sol1
, ,
r = (proc ( ) ... end procr ) Q( )r = (proc ( ) ... end procr ) = d
d
rQ( )r (proc ( ) ... end procr )
An example of the procedure list for the value T is the following:
=
( )
sol1500000
, ,
=
( )
r 500000 500000 Q( )r(500000) = 1.00000000059982908 ( ) =
d d
rQ( )r 500000 -0.306523319533648560 10-14
To acquire the data from the list, operators have to be used. OP[1] is the operator for the ra-dial coordinate. OP[2] is the function Q(r) and OP[3] is its derivate. For the substitution of Q(T) the operator OP[2](T) is needed. When the cycle stops K is determined to the given pre-cision (in this case 10-8). With the acquired initial values the ODE can be solved and the solu-tion will satisfy the original condisolu-tion given by equasolu-tion (110). Equations (98)-(101) can be solved with the same procedure. For the equations obtained from different elastic potentials (specialised version of the general Blatz-Ko potential) the constant K will be different. The difference is due to the fact that these equations include terms of the function Q(r) of different degree and were determined from different invariants. Some values of K for the general and the two special cases obtained from the routine for a 1 mm diameter needle, taking ρ0=0.45 mm is given here
:=
KGenBK 4.826931164 stepnumber := 29 :=
KSpecBK1 4.327344839 stepnumber := 30 :=
KSpecBK2 12.11344370 stepnumber := 27
Values of K for the three cases including the step number of the iterative calculation
The integral curve (solutions) for the three cases are shown in the next figure (Fig. 36)
Q(r)
r
Fig. 36: The numerical solutions of the equations of equilibrium of the three elastic potential for a 1 mm diameter needle, with determined initial values
The curves show the same characteristics that were determined in the previous analysis, this way they can be considered as the solution of the equations. These statements are also valid for the results obtained from the other two material laws. The equations of equilibrium of these cases are similar to the general one with the exception that they contain terms of lower degree. The characteristics of the solution functions are the same as in the general case, with the differences that the value of K of the initial value for Q’ is different in every case.
Notes
The calculations showed that the integration constants are not independent of each other in each case. They both have to be determined simultaneously. In some cases when the initial condition for Q(0) is not prescribed properly (e.g.: too large number is given for it, meaning a very large initial deformation) the other K cannot be determined. It would lead to a value of the initial slope too large that the even the mathematical software could not handle.