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}~.}