 So, welcome to lecture 39. Here in this lecture we will see how do we use finite element method to calculate forces in electrical machines and equipment. So, what we will do is we will see two types of four calculations. One type we will see in L39 and other calculation we will see in L40. So, the first in this lecture we will see Maxwell stress tensor best calculation. Now, let us see some basics. So, mainly two types of force calculations are required. One is total force for rigidly moving system that means a body is rigid and it does not get locally deformed and body as a whole moves. So, that is called as rigidly moving system. Second is local force calculation wherein you have local deformations and then you need to calculate forces in a different way and you know at local level and how to do that we will see. Now, these you know two types of forces can be obtained using the following methods. So, total force the first type can be obtained by using either Lorentz force density that is J cross B, Maxwell stress tensor, MST, virtual work that is VW or Coulomb virtual work CVW methods whereas local force is obtained by J cross B, virtual work and CVW methods. MST cannot be used for local force because in MST based calculation we integrate you know force along some contour. So, obviously you cannot you are not then calculating local forces there. So, now you know secondly forces on current carrying conductors they can be calculated either by using J cross B method or MST or virtual work or Coulomb virtual work method. Whereas you know saliency or reluctance forces are calculated by MST virtual work and CVW method. J cross B method cannot be used here because J is unknown in the structure experiencing the force. MST, J cross B and CVW methods require only one FEM solution whereas virtual work method requires two FEM solutions. Now, one may wonder if virtual work the CVW is as the name suggests is an extension of virtual work then how come this requires two whereas CVW requires one FEM solution. See in virtual work technique we require you know two FEM solutions because we want to calculate force as change in energy divided by displacement when the body under consideration moves from one position to another hence we require energy values at both positions requiring two FEM solutions. Whereas in case of CVW that differentiation DW by DX X being the displacement DW by DX is taken care in the FEM formulation itself and that is why we require only one FEM solution. Of course we are not going to see this Coulomb virtual work method in this lecture anyone in the lecture but you know those who are interested they can see literature to understand more on Coulomb virtual work method. Now, let us go into details of Maxwell stress tensor method. Now, to calculate the force on the body what we have to do we have to enclose it by a surface which is a hypothetical surface and then integrate Maxwell stress tensor on this surface. Now, what is this MST we will see in the coming slides and in the 2D formulation the surface integral reduces to a counter integral. Obviously, from 3D to 2D the complexity reduces from surface integration to counter integration. Accuracy of the MST method depends on the choice of the integrating surface. So, one has to be careful and if there are ion air interfaces then we have to use very small size elements to reduce errors. So, this method can be used for both ferromagnetic materials and current carrying conductors as was explained on the first slide. So, now some theory about virtual work method although the application we are going to see in the next lecture that is L40. So, in this virtual work method the movable part of the system is virtually displaced in simulation and change in energy is calculated. If the magnetic energy Wf is used for calculation flux linkage lambda is kept constant in the first approach whereas in the second approach if magnetic co-energy Wco is used for calculation then current I is kept constant. So, this is second approach. So, both the approaches we will discuss in the next lecture. So, this is what is you know summarized these two statements are summarized by this you know statements here mathematical you know statements. So, f is dwf by dx wf by dx at flux constant condition can give you force or you can get the force value by Wco by Wx with current held constant. This method is less sensitive to choice of mesh and is fast converging as compared to the MST method and that is why it is many times preferred for force calculation. Now, let us understand what is this Maxwell stress tensor. Now, we know J cross B is the force density force per unit volume and J is del cross H and we are of course this neglecting here displacement current density. So, del cross H is equal to J and now H is replaced by B upon mu naught. So, fv force per unit volume is given by del cross B upon mu naught cross B wherein J is replaced by this expression. So, now first we let us calculate del cross B upon mu naught which is given by this just 1 over mu naught is taken out and then this is curl of B and when you expand curl of B you will get this you know expression. And now if you use this expression here these components x ax hat ay hat and az hat components in x, y and z directions these 3 components if you use now here in this cross product because this is what finally we want to calculate estimate. So, del cross B upon mu naught cross B will be given by these 3 components written here and then B x B y B z corresponding to this B. Now, you know I am going to skip some of the intermediate steps but for those who are interested in knowing all steps they can refer this book wherein detailed derivation is given. When you expand it then the x component you will get as this. So, when you further actually you know manipulate that previous expression then you get this expression and noting that divergence B is equal to 0 as always and mod of B square is B x square plus B y square plus B z square right. We can rewrite this whole thing as this. Now, this has been sort of rewritten in matrix multiplication form a row vector is multiplied with a column vector. Now, actually speaking this is really you know what you see in the bracket is basically del operator written in row matrix form. So, actually what is implicit here is A x hat here, A y hat here and A z hat here and similarly you have A x hat here, A y hat and A z hat. So, this is what you know but when we write it in matrix multiplication form we actually do not write A x, A y and A z hats del dot this vector gives you this. Similarly, F y you can get it as this and F z you get it as this. So, now you can write complete force. Now, force has three components F x, F y and F z. So, as I mentioned earlier this when you write it in this kind of matrix row matrix row vector form then this is nothing but actually F x, A x hat plus F y, A y hat plus F z, A z hat. So, this is nothing but force right but since we are writing in matrix form we can drop those unit vectors. So, now if we combine all three expressions that we obtained in the previous slide then we can actually combine those individual three column matrices and then we get this 3 by 3 matrix. This I have already explained this is del operator and now we can then write this F x, F y, F z as divergence of tensor t where t is this entire matrix. Now, here divergence of this column which is a vector is the scalar F x, divergence of this is F y, divergence of this is F z right. So, remember again just I am repeating this is in fact this is this multiplied by A x hat this is this into A y hat and A z hat and when you take divergence then you will get you know F x as was you know obtained earlier. So, now this capital F force vector is volume integral F v dv right and then F v we can write this as divergence of tensor dv and then we basically you know invoke divergence theorem and then write this as close surface integral t dot ds and where in df is equal to t dot ds. And in case of 2d this df will become t dot A n hat dl where ds is dl you know ds vector is ds as a magnitude into A n hat the corresponding vector representing ds. And then area ds magnitude is nothing but dl into 1 since we are doing 2d approximation you know it will become you know just dl. So, here the divergence t therefore is a vector F x, F y, F z 3 components are there. Now, let us understand you know little bit more about this you know tangential and normal unit vectors which will be required you know in further calculations for example, here you require unit normal vector. Now, how do we basically you know express it mathematically we will see now. Now, let us take you know an elemental edge in 2d a plane reduces to an edge. Now, let us take this this you know surface and then this actually this is a contour. So, now we are going from 3d to 2d and then we will see how do we express A n hat and A t hat tangential and normal vectors unit normal vectors. So, A tangential unit vector can be written as L x A x hat plus L y A y hat wherein since this is a unit vector square root of L x square plus L y square is equal to 1 right L x and L y individually can be positive or negative right. So, it depends on how this A tangential is oriented. So, A normal hat will be then given by minus L y A x hat plus L x A y hat. Now, if you take the dot product of A t hat and A n hat then you will find that the dot product is 0 because this will give minus L x L y plus L x L y so dot product is 0 as it should be because they are orthogonal vectors right. So, now how do we find out these you know values of L x and L y, L x and L y for each of these vectors. For the surface under consideration which is an agent. Then we have to worry about how do we express A t and A n in terms of you know chosen coordinate system. So, now let us see as an example this you know A t hat along this edge and we are considering this corresponding surface. So, A n hat for the entire this surface will be in this direction. So, now if we bring this A n hat at this point the point under consideration on this edge remember A n hat is you know at every point on this surface A n hat is normal and it is in the same direction. So, this A n hat as a representation is shown at the center of this center of this face but I can always bring this A n hat at the point under consideration on this edge. Now actually if you bring it here and now if you actually you know make a plane pass through you know A t and A at hat A t hat and A n hat and we call that as you know x y plane. Then you will get A t hat and A n hat as shown here and then you can see here for this A t hat both L x and L y will be negative because L y is x is negative and y also is negative. So, both these components will be negative and A n hat L x is positive and L y is negative. And then that ensures when these both are both these are negative and here you know one is positive and one is negative it ensures that their dot product is 0. So, now I am just writing these expressions for the two unit vectors again here and now since we are considering 2D approximation that tensor matrix reduces to only 2 by 2 and then you know this is dL is this A n hat is this minus L y and L x right and again we are not since we are writing it is a column matrix we are not writing unit vectors. Similarly, there are unit vectors here also A x hat, A y hat, A x hat, A y hat but we are not writing those. So, now if we expand this you will get dFx and dFy as these two expressions a simple just a product of these two you will get these two expressions. So, going further now let us see how do we calculate force between two circular conductors first we will see by J cross B approach and then we will also see how to calculate force using MST approach. Now these B y and B expressions are very well known to you we have seen these many number of times. So, B y is given by this and B x is given by this exploration here on the right hand side and then the B net is square root of B x square plus B y square and for this we are talking you know for one element we always run for loop and then for each element we calculate some performance parameter same thing is being done here. Now what we do we find out you know whether the element under consideration in the for loop whether it is you know lying in this first conductor. Now what we are doing these are the two circular conductors here 4 and 5 and then you know we have enclosed them into you know rectangular boundaries which will be required basically for you know the MST best calculation force calculation. So, now if the you know the element is in sub domain 4 then f x and f y are given by these two you know expressions how do we get them f is nothing but j cross B is the force per unit volume. So, that you multiply by volume area of the element into 1 because it is 2 dimensional approximation with z dimension taken as 1 meter. So, j cross B is this remembering that j only has z component. So, if you take B as only having B x and B y components and j with only j z component you can easily you know verify that j cross B will come as this. So, now here f x the x component will be j j B y minus j B y into area of the element. Similarly, f y will be j B x. So, j B x into area of the element. So, you calculate you know these forces for you know this conductor and call them as f x 1 and f y 1 for each element and actually you know go on you know go to the next element and keep doing this and store all the force values in you know these column vectors x and y components. Similar thing you do for second circular conductor or remember this is a circular but this is a cylindrical conductor. So, z direction is not shown here it is a 2 dimensional approximation. So, then you calculate and store values of forces on all elements within these finite elements within this conductor phi in these 2 column vectors x and y components of forces and then you basically sum all the forces. So, this sum of f x 1 will give you the total x component of x component of force acting on conductor 4 the resultant y component of force acting on conductor 4 similarly for conductor phi. So, by this we have got f x and f y acting on individual these 2 conductors. So, now let us see how do we calculate force between the same 2 circular conductors by using the MST approach. So, what do we do we are here we are just initializing by these 2 commands. Then what we do is these 2 conductors and as I as I as I mentioned previously these 2 conductors are enclosed in these 2 fictitious counters over which we will do integration of Maxwell stress tensor. And now what we do we run from run a loop from i equal to 1 to number of edges. Now what are these n underscore edges? Now there are these 4 main edges here 9, 10, 11, 12. Now if you take age 9 then age 9 will have number of sub edges because now you know you are here you will have one triangular element like this, other triangular element like this, third triangular element like this and so on. So, this will form these segments will form sides of triangular elements like this. So, you can call these as sub edges on this main age 9. So, when we say i goes from 1 to n underscore edges we are basically going over all these sub edges. When we are running this for loop if the sub edge under consideration if it is lying either on 9, 10, 11 or 12th main edge then what you do in nodes 1 column vector you take using the E matrix you take the coordinates which are there in second and third row of the E matrix corresponding to the age i because each age i will have 2 nodes and these are these nodes. So, you have got by using this command node numbers of end points of each sub edge and then using these 2 commands you basically get the x and y coordinates of those end points of the sub edge under consideration. And then you take average of the x and y coordinates of the end points and then you take the you know by using this command you find out length of that sub edge. Then you find out in which triangular element means midpoints of the sub edge under consideration lie having got that triangular element number in which that you know midpoint lies then we can get the corresponding b x and b y values corresponding to that triangle and then we can calculate net value of b which is square root of b x square plus b y square. So, now let us you know see how we can calculate forces now depending upon you know sub edge where actually you know it is lying whether it is part of main age 9 or 10 or 11 or 12 we will execute corresponding sets of commands. Now, if for example for age 9 lx is 1 and ly is 0 why because we have seen here a hat tangential is given by lx ax hat plus ly a y hat and this relation holds. So, now here for any sub edge on this age 9 which is along x direction lx will be 1 and ly will be 0 if lx is 1 and ly is 0 then only this component will be there for any age along this any sub edge along this 9. So, that is the reason that here it is this is 1. So, it is by bx upon mu 0 into dl this is dl divided by mu 0 which I already explained and then similarly the y component y component will be again for this sub edge it is 1 and 0 so ly is 0. So, this term goes to 0 in case of dfy and this term remains. So, it is by square minus 0.5 b net square into dl divided by mu 0 and anyway lx is 1. So, this is how you can calculate x and y components for any sub edge which is part of main age 9. Similarly, then you can calculate for you know all other sub edges if they are lying on either 10, 11 or 12th main edge and remember for age 10 which is vertically down lx will be 0, ly will be minus 1 because it is in negative y direction similarly for 11 and 12. Now the FEM solution we will see now the solution by both these methods Lorentz force equation and Maxwell stress tensor and also we will actually compare with analytical formula. So, Lorentz force equation j cross b using that method which was already described you get force values as this whereas Maxwell stress tensor method gives you forces on conductor 1 and conductor 2 as this and these values are quite close to each other confirming that our MST approach is you know validated and for these two conductors we are assuming that currents are in opposite directions the field will be as shown like this and since the currents are in opposite directions these conductors will repel each other and that is the reason we are getting forces in opposite directions. So, conductor 1 force is in negative x direction and for conductor 2 the force is in positive x direction. The analytical formula is mu 0 i1 i2 upon 2 pi d into the length in z direction so this gives the force between the two parallel conductors. Now here since we are calculating you know per unit depth because it is two-dimensional approximation so the length in z direction is 1 so this actually here into 1 so that is that is why it is not considered. So, now actually if you substitute the values then you get this as 1.03 Newton per meter. Now here diameter of two conductors is taken as 10 mm and current density assumed is 5 ampere per mm square so distance between two conductors is 30 mm. So, i1 and i2 will be j into area so pi d square by 4 where d is 0.01 meter so you get 392.7 amperes that where you substituted here for both this i1 and i2 and the distance as 0.03 you get this 1.03 which again is quite close to these values. So, all the three you know methods analytical formula base calculation, MSD base calculation and Lorez force equation base calculation all the three are giving values in the same range. So, this verifies our you know calculation methods. Finally, we actually also calculate force on a core or plunger by using MSD method. Now here this is one core C shape core and this is again a core piece which you can call it as a plunger and when there is a coil here when you pass current to this coil you get flux established like this and now this you know this core C shape core will attract this piece of core because you know if it is north pole here there will be south pole here and then these two opposite poles will attract each other same thing will happen here. So, now here the MSD in MSD base calculation what we do we enclose this piece by using by you know a fictitious contour and we integrate MSD over this contour. So, the dimensions are given here air gap length is 40 mm. So, this is 40 mu R is mu R for this and this is taken as 5000 and the other dimensions are given in this figure right and coil dimensions are also given. So, when we actually you know calculate force using the MSD method which is which was described in this lecture we get force as 2.71 into 10 raise to 4 Newton per meter again remember this is per meter depth because we are doing 2D approximation. So, we do calculations for unit depth in Z direction. So, we basically we will calculate the same force value for the last case study using virtual work method in the next lecture. Thank you.