Use Maple to reduce your manual labor.
The problem becomes easier to solve if we consider numerical values of the parameters. Let the local nodes numbers of element 5 {\displaystyle 5} be 1 {\displaystyle 1} for node 5 {\displaystyle 5} , and 2 {\displaystyle 2} for node 6 {\displaystyle 6} .
Let us assume that the beam is divided into six equal sectors. Then,
θ = π 4 ; θ 1 = π 3 ; θ 2 = π 6 . {\displaystyle \theta ={\cfrac {\pi }{4}}~;~~\theta _{1}={\cfrac {\pi }{3}}~;~~\theta _{2}={\cfrac {\pi }{6}}~.} Let r 1 = 1 {\displaystyle r_{1}=1} and r 2 = 1.2 {\displaystyle r_{2}=1.2} . Since the blue point is midway between the two, r = 1.1 {\displaystyle r=1.1} .
Also, let the components of the stiffness matrix of the composite be
C 11 = 10 ; C 33 = 100 ; C 12 = 6 ; C 13 = 60 ; C 44 = 30 . {\displaystyle C_{11}=10;~~C_{33}=100;~~C_{12}=6;~~C_{13}=60;~~C_{44}=30~.} Let the velocities for nodes 1 {\displaystyle 1} and 2 {\displaystyle 2} of the element be
v 1 = [ v 1 x v 1 y ω 1 ] = [ 1 2 1 ] ; v 2 = [ v 2 x v 2 y ω 2 ] = [ 2 1 1 ] {\displaystyle \mathbf {v} _{1}={\begin{bmatrix}v_{1}^{x}\\v_{1}^{y}\\\omega _{1}\end{bmatrix}}={\begin{bmatrix}1\\2\\1\end{bmatrix}}~;~~\mathbf {v} _{2}={\begin{bmatrix}v_{2}^{x}\\v_{2}^{y}\\\omega _{2}\end{bmatrix}}={\begin{bmatrix}2\\1\\1\end{bmatrix}}} The x {\displaystyle x} and y {\displaystyle y} coordinates of the master and slave nodes are
[ x 1 y 1 x 2 y 2 ] = [ r cos θ 1 r sin θ 1 r cos θ 2 r sin θ 2 ] {\displaystyle {\begin{bmatrix}x_{1}\\y_{1}\\x_{2}\\y_{2}\end{bmatrix}}={\begin{bmatrix}r\cos \theta _{1}\\r\sin \theta _{1}\\r\cos \theta _{2}\\r\sin \theta _{2}\end{bmatrix}}}
[ x 1 − y 1 − x 2 − y 2 − ] = [ r 1 cos θ 1 r 1 sin θ 1 r 1 cos θ 2 r 1 sin θ 2 ] {\displaystyle {\begin{bmatrix}x_{1-}\\y_{1-}\\x_{2-}\\y_{2-}\end{bmatrix}}={\begin{bmatrix}r_{1}\cos \theta _{1}\\r_{1}\sin \theta _{1}\\r_{1}\cos \theta _{2}\\r_{1}\sin \theta _{2}\end{bmatrix}}}
[ x 1 + y 1 + x 2 + y 2 + ] = [ r 2 cos θ 1 r 2 sin θ 1 r 2 cos θ 2 r 2 sin θ 2 ] {\displaystyle {\begin{bmatrix}x_{1+}\\y_{1+}\\x_{2+}\\y_{2+}\end{bmatrix}}={\begin{bmatrix}r_{2}\cos \theta _{1}\\r_{2}\sin \theta _{1}\\r_{2}\cos \theta _{2}\\r_{2}\sin \theta _{2}\end{bmatrix}}} Since there are two master nodes in an element, the parent element is a four-noded quad with shape functions
N 1 − ( ξ , η ) = 1 4 ( 1 − ξ ) ( 1 − η ) N 2 − ( ξ , η ) = 1 4 ( 1 + ξ ) ( 1 − η ) N 1 + ( ξ , η ) = 1 4 ( 1 − ξ ) ( 1 + η ) N 2 + ( ξ , η ) = 1 4 ( 1 + ξ ) ( 1 + η ) . {\displaystyle {\begin{aligned}N_{1-}(\xi ,\eta )&={\cfrac {1}{4}}(1-\xi )(1-\eta )&N_{2-}(\xi ,\eta )&={\cfrac {1}{4}}(1+\xi )(1-\eta )\\N_{1+}(\xi ,\eta )&={\cfrac {1}{4}}(1-\xi )(1+\eta )&N_{2+}(\xi ,\eta )&={\cfrac {1}{4}}(1+\xi )(1+\eta )~.\end{aligned}}} The isoparametric map is
x ( ξ , η ) = x 1 − N 1 − ( ξ , η ) + x 2 − N 2 − ( ξ , η ) + x 1 + N 1 + ( ξ , η ) + x 2 + N 2 + ( ξ , η ) y ( ξ , η ) = y 1 − N 1 − ( ξ , η ) + y 2 − N 2 − ( ξ , η ) + y 1 + N 1 + ( ξ , η ) + y 2 + N 2 + ( ξ , η ) . {\displaystyle {\begin{aligned}x(\xi ,\eta )&=x_{1-}~N_{1^{-}}(\xi ,\eta )+x_{2-}~N_{2^{-}}(\xi ,\eta )+x_{1+}~N_{1^{+}}(\xi ,\eta )+x_{2+}~N_{2^{+}}(\xi ,\eta )\\y(\xi ,\eta )&=y_{1-}~N_{1^{-}}(\xi ,\eta )+y_{2-}~N_{2^{-}}(\xi ,\eta )+y_{1+}~N_{1^{+}}(\xi ,\eta )+y_{2+}~N_{2^{+}}(\xi ,\eta )~.\end{aligned}}} Therefore, the derivatives with respect to ξ {\displaystyle \xi } are
x , ξ = ∂ x ∂ ξ y , ξ = ∂ y ∂ ξ {\displaystyle {\begin{aligned}x_{,\xi }={\frac {\partial x}{\partial \xi }}\\y_{,\xi }={\frac {\partial y}{\partial \xi }}\end{aligned}}} Since the blue point is at the center of the element, the values of ξ {\displaystyle \xi } and η {\displaystyle \eta } at that point are zero.
Therefore,
x , ξ = x , ξ e x − y , ξ e y , ‖ x , ξ ‖ = x , ξ 2 + y , ξ 2 . {\displaystyle \mathbf {x} _{,\xi }=x_{,\xi }\mathbf {e} _{x}-y_{,\xi }\mathbf {e} _{y},~~\lVert \mathbf {x} _{,\xi }\rVert _{}={\sqrt {x_{,\xi }^{2}+y_{,\xi }^{2}}}~.} The local laminar basis vector e ^ x {\displaystyle {\widehat {\mathbf {e} }}_{x}} is given by
e ^ x = x , ξ ‖ x , ξ ‖ {\displaystyle {\widehat {\mathbf {e} }}_{x}={\cfrac {\mathbf {x} _{,\xi }}{\lVert \mathbf {x} _{,\xi }\rVert _{}}}} The laminar basis vector e ^ y {\displaystyle {\widehat {\mathbf {e} }}_{y}} is given by
e ^ y = e z × e x ^ {\displaystyle {\widehat {\mathbf {e} }}_{y}=\mathbf {e} _{z}\times {\widehat {\mathbf {e} _{x}}}} To compute the velocity gradient, we have to find the velocities at the slave nodes using the relation
[ v i − x v i − y v i + x v i + y ] = [ 1 0 − ( y i − − y i ) 0 1 ( x i − − x i ) 1 0 − ( y i + − y i ) 0 1 ( x i + − x i ) ] [ v i x v i y ω i ] {\displaystyle {\begin{bmatrix}v_{i-}^{x}\\v_{i-}^{y}\\v_{i+}^{x}\\v_{i+}^{y}\end{bmatrix}}={\begin{bmatrix}1&0&-(y_{i-}-y_{i})\\0&1&(x_{i-}-x_{i})\\1&0&-(y_{i+}-y_{i})\\0&1&(x_{i+}-x_{i})\end{bmatrix}}{\begin{bmatrix}v_{i}^{x}\\v_{i}^{y}\\\omega _{i}\end{bmatrix}}} For master node 1 of the element (global node 5), the velocities of the slave nodes are
[ v 1 − x v 1 − y v 1 + x v 1 + y ] {\displaystyle {\begin{bmatrix}v_{1-}^{x}\\v_{1-}^{y}\\v_{1+}^{x}\\v_{1+}^{y}\end{bmatrix}}} For master node 2 of the element (global node 6), the velocities of the slave nodes are
[ v 2 − x v 2 − y v 2 + x v 2 + y ] {\displaystyle {\begin{bmatrix}v_{2-}^{x}\\v_{2-}^{y}\\v_{2+}^{x}\\v_{2+}^{y}\end{bmatrix}}} The interpolated velocity within the element is given by
v ( ξ , t ) = D D t [ x ( ξ , t ) ] = ∑ i − = 1 2 ∂ ∂ t [ x i − ( t ) ] N i − ( ξ , η ) + ∑ i + = 1 2 ∂ ∂ t [ x i + ( t ) ] N i + ( ξ , η ) = ∑ i − = 1 2 v i − ( t ) N i − ( ξ , η ) + ∑ i + = 1 2 v i + ( t ) N i + ( ξ , η ) . {\displaystyle {\begin{aligned}\mathbf {v} ({\boldsymbol {\xi }},t)&={\cfrac {D}{Dt}}[\mathbf {x} ({\boldsymbol {\xi }},t)]\\&=\sum _{i-=1}^{2}{\frac {\partial }{\partial t}}[\mathbf {x} _{i-}(t)]~N_{i^{-}}(\xi ,\eta )+\sum _{i+=1}^{2}{\frac {\partial }{\partial t}}[\mathbf {x} _{i+}(t)]~N_{i^{+}}(\xi ,\eta )\\&=\sum _{i-=1}^{2}\mathbf {v} _{i-}(t)~N_{i-}(\xi ,\eta )+\sum _{i+=1}^{2}\mathbf {v} _{i+}(t)~N_{i+}(\xi ,\eta )~.\end{aligned}}} The velocity gradient is given by
∇ v = v i , j = ∂ v i ∂ x j . {\displaystyle {\boldsymbol {\nabla }}\mathbf {v} =v_{i,j}={\frac {\partial v_{i}}{\partial x_{j}}}~.} The velocities are given in terms of the parent element coordinates (ξ , η {\displaystyle \xi ,\eta } ). We have to convert them to the (x , y {\displaystyle x,y} ) system in order to compute the velocity gradient. To do that we recall that
∂ v x ∂ ξ = ∂ v x ∂ x ∂ x ∂ ξ + ∂ v x ∂ y ∂ y ∂ ξ ∂ v x ∂ η = ∂ v x ∂ x ∂ x ∂ η + ∂ v x ∂ y ∂ y ∂ η ∂ v y ∂ ξ = ∂ v y ∂ x ∂ x ∂ ξ + ∂ v y ∂ y ∂ y ∂ ξ ∂ v y ∂ η = ∂ v y ∂ x ∂ x ∂ η + ∂ v y ∂ y ∂ y ∂ η . {\displaystyle {\begin{aligned}{\frac {\partial v_{x}}{\partial \xi }}&={\frac {\partial v_{x}}{\partial x}}{\frac {\partial x}{\partial \xi }}+{\frac {\partial v_{x}}{\partial y}}{\frac {\partial y}{\partial \xi }}\\{\frac {\partial v_{x}}{\partial \eta }}&={\frac {\partial v_{x}}{\partial x}}{\frac {\partial x}{\partial \eta }}+{\frac {\partial v_{x}}{\partial y}}{\frac {\partial y}{\partial \eta }}\\{\frac {\partial v_{y}}{\partial \xi }}&={\frac {\partial v_{y}}{\partial x}}{\frac {\partial x}{\partial \xi }}+{\frac {\partial v_{y}}{\partial y}}{\frac {\partial y}{\partial \xi }}\\{\frac {\partial v_{y}}{\partial \eta }}&={\frac {\partial v_{y}}{\partial x}}{\frac {\partial x}{\partial \eta }}+{\frac {\partial v_{y}}{\partial y}}{\frac {\partial y}{\partial \eta }}~.\end{aligned}}} In matrix form
[ ∂ v x ∂ ξ ∂ v x ∂ η ] = [ ∂ x ∂ ξ ∂ y ∂ ξ ∂ x ∂ η ∂ y ∂ η ] [ ∂ v x ∂ x ∂ v x ∂ y ] = J [ ∂ v x ∂ x ∂ v x ∂ y ] {\displaystyle {\begin{bmatrix}{\frac {\partial v_{x}}{\partial \xi }}\\{\frac {\partial v_{x}}{\partial \eta }}\end{bmatrix}}={\begin{bmatrix}{\frac {\partial x}{\partial \xi }}&{\frac {\partial y}{\partial \xi }}\\{\frac {\partial x}{\partial \eta }}&{\frac {\partial y}{\partial \eta }}\\\end{bmatrix}}{\begin{bmatrix}{\frac {\partial v_{x}}{\partial x}}\\{\frac {\partial v_{x}}{\partial y}}\end{bmatrix}}=\mathbf {J} {\begin{bmatrix}{\frac {\partial v_{x}}{\partial x}}\\{\frac {\partial v_{x}}{\partial y}}\end{bmatrix}}} and
[ ∂ v y ∂ ξ ∂ v y ∂ η ] = [ ∂ x ∂ ξ ∂ y ∂ ξ ∂ x ∂ η ∂ y ∂ η ] [ ∂ v y ∂ x ∂ v y ∂ y ] = J [ ∂ v y ∂ x ∂ v y ∂ y ] . {\displaystyle {\begin{bmatrix}{\frac {\partial v_{y}}{\partial \xi }}\\{\frac {\partial v_{y}}{\partial \eta }}\end{bmatrix}}={\begin{bmatrix}{\frac {\partial x}{\partial \xi }}&{\frac {\partial y}{\partial \xi }}\\{\frac {\partial x}{\partial \eta }}&{\frac {\partial y}{\partial \eta }}\\\end{bmatrix}}{\begin{bmatrix}{\frac {\partial v_{y}}{\partial x}}\\{\frac {\partial v_{y}}{\partial y}}\end{bmatrix}}=\mathbf {J} {\begin{bmatrix}{\frac {\partial v_{y}}{\partial x}}\\{\frac {\partial v_{y}}{\partial y}}\end{bmatrix}}~.} Therefore,
[ ∂ v x ∂ x ∂ v x ∂ y ] = J − 1 [ ∂ v x ∂ ξ ∂ v x ∂ η ] and [ ∂ v y ∂ x ∂ v y ∂ y ] = J − 1 [ ∂ v y ∂ ξ ∂ v y ∂ η ] . {\displaystyle {\begin{bmatrix}{\frac {\partial v_{x}}{\partial x}}\\{\frac {\partial v_{x}}{\partial y}}\end{bmatrix}}=\mathbf {J} ^{-1}{\begin{bmatrix}{\frac {\partial v_{x}}{\partial \xi }}\\{\frac {\partial v_{x}}{\partial \eta }}\end{bmatrix}}\qquad {\text{and}}\qquad {\begin{bmatrix}{\frac {\partial v_{y}}{\partial x}}\\{\frac {\partial v_{y}}{\partial y}}\end{bmatrix}}=\mathbf {J} ^{-1}{\begin{bmatrix}{\frac {\partial v_{y}}{\partial \xi }}\\{\frac {\partial v_{y}}{\partial \eta }}\end{bmatrix}}~.} The rate of deformation is given by
D = 1 2 [ ∇ v + ( ∇ v ) T ] . {\displaystyle {\boldsymbol {D}}={\frac {1}{2}}[{\boldsymbol {\nabla }}\mathbf {v} +({\boldsymbol {\nabla }}\mathbf {v} )^{T}]~.} The global base vectors are
e x = [ 1 0 ] ; e y = [ 0 1 ] . {\displaystyle \mathbf {e} _{x}={\begin{bmatrix}1\\0\end{bmatrix}};~\qquad \mathbf {e} _{y}={\begin{bmatrix}0\\1\end{bmatrix}}~.} Therefore, the rotation matrix is
R lam = [ e x ∙ e x ^ e x ∙ e y ^ e y ∙ e x ^ e y ∙ e y ^ ] {\displaystyle \mathbf {R} _{\text{lam}}={\begin{bmatrix}\mathbf {e} _{x}\bullet {\widehat {\mathbf {e} _{x}}}&\mathbf {e} _{x}\bullet {\widehat {\mathbf {e} _{y}}}\\\mathbf {e} _{y}\bullet {\widehat {\mathbf {e} _{x}}}&\mathbf {e} _{y}\bullet {\widehat {\mathbf {e} _{y}}}\end{bmatrix}}} Therefore, the components of the rate of deformation tensor with respect to the laminar coordinate system are
D lam = R lam T D R lam {\displaystyle \mathbf {D} _{\text{lam}}=\mathbf {R} _{\text{lam}}^{T}\mathbf {D} \mathbf {R} _{\text{lam}}} The rate constitutive relation of the material is given by
D D t [ σ 11 σ 22 σ 33 σ 23 σ 31 σ 12 ] = [ C 11 C 12 C 13 0 0 0 C 12 C 11 C 13 0 0 0 C 13 C 13 C 33 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 44 0 0 0 0 0 0 ( C 11 − C 12 ) / 2 ] [ D 11 D 22 D 33 D 23 D 31 D 12 ] . {\displaystyle {\cfrac {D}{Dt}}{\begin{bmatrix}\sigma _{11}\\\sigma _{22}\\\sigma _{33}\\\sigma _{23}\\\sigma _{31}\\\sigma _{12}\end{bmatrix}}={\begin{bmatrix}C_{11}&C_{12}&C_{13}&0&0&0\\C_{12}&C_{11}&C_{13}&0&0&0\\C_{13}&C_{13}&C_{33}&0&0&0\\0&0&0&C_{44}&0&0\\0&0&0&0&C_{44}&0\\0&0&0&0&0&(C_{11}-C_{12})/2\end{bmatrix}}{\begin{bmatrix}D_{11}\\D_{22}\\D_{33}\\D_{23}\\D_{31}\\D_{12}\end{bmatrix}}~.} Since the problem is a 2-D one, the reduced constitutive equation is
D D t [ σ 11 σ 33 σ 31 ] = [ C 11 C 13 0 C 13 C 33 0 0 0 C 44 ] [ D 11 D 33 D 31 ] . {\displaystyle {\cfrac {D}{Dt}}{\begin{bmatrix}\sigma _{11}\\\sigma _{33}\\\sigma _{31}\end{bmatrix}}={\begin{bmatrix}C_{11}&C_{13}&0\\C_{13}&C_{33}&0\\0&0&C_{44}\end{bmatrix}}{\begin{bmatrix}D_{11}\\D_{33}\\D_{31}\end{bmatrix}}~.} The laminar x {\displaystyle x} -direction maps to the composite 3 {\displaystyle 3} -direction and the laminar y {\displaystyle y} -directions maps to the composite 1 {\displaystyle 1} -direction. Hence the constitutive equation can be written as
D D t [ σ y y σ x x σ x y ] = [ C 11 C 13 0 C 13 C 33 0 0 0 C 44 ] [ D y y D x x D x y ] . {\displaystyle {\cfrac {D}{Dt}}{\begin{bmatrix}\sigma _{yy}\\\sigma _{xx}\\\sigma _{xy}\end{bmatrix}}={\begin{bmatrix}C_{11}&C_{13}&0\\C_{13}&C_{33}&0\\0&0&C_{44}\end{bmatrix}}{\begin{bmatrix}D_{yy}\\D_{xx}\\D_{xy}\end{bmatrix}}~.} Rearranging,
D D t [ σ x x σ y y σ x y ] = [ C 33 C 13 0 C 13 C 11 0 0 0 C 44 ] [ D x x D y y D x y ] . {\displaystyle {\cfrac {D}{Dt}}{\begin{bmatrix}\sigma _{xx}\\\sigma _{yy}\\\sigma _{xy}\end{bmatrix}}={\begin{bmatrix}C_{33}&C_{13}&0\\C_{13}&C_{11}&0\\0&0&C_{44}\end{bmatrix}}{\begin{bmatrix}D_{xx}\\D_{yy}\\D_{xy}\end{bmatrix}}~.} The plane stress condition requires that σ y y = 0 {\displaystyle \sigma _{yy}=0} in the laminar coordinate system. We assume that the rate of σ y y {\displaystyle \sigma _{yy}} is also zero. In that case, we get
0 = D D t ( σ y y ) = C 13 D x x + C 11 D y y {\displaystyle 0={\cfrac {D}{Dt}}(\sigma _{yy})=C_{13}~D_{xx}+C_{11}~D{yy}} or,
D y y = − C 13 C 11 D x x . {\displaystyle D_{yy}=-{\cfrac {C_{13}}{C_{11}}}~D_{xx}~.} To get the global stress rate and rate of deformation, we have to rotate the components to the global basis using
D D t ( σ ) = R lam D D t ( σ lam ) R lam T ; D = R lam D lam R lam T . {\displaystyle {\cfrac {D}{Dt}}({\boldsymbol {\sigma }})=\mathbf {R} _{\text{lam}}{\cfrac {D}{Dt}}({\boldsymbol {\sigma }}_{\text{lam}})\mathbf {R} _{\text{lam}}^{T}~;\qquad \mathbf {D} =\mathbf {R} _{\text{lam}}\mathbf {D} _{\text{lam}}\mathbf {R} _{\text{lam}}^{T}~.}